Quantifying the Stellar Halo's Response to the LMC's Infall with Spherical Harmonics

A fundamental challenge facing the study of the Milky Way (MW) galaxy is that most of its mass is in dark matter (DM). Because we cannot directly observe the MW's DM halo, we must use tracer populations, such as halo stars, globular clusters, and MW satellites, to study the MW's DM halo indirectly. In particular, the velocity distributions of tracer populations can be used to derive estimates of the MW's mass, utilizing methods derived from Jeans (1915) modeling (e.g., Dehnen et al. 2006; Watkins et al. 2009, 2019; Gnedin et al. 2010; Deason et al. 2012; Eadie et al. 2017; Sohn et al. 2018; Wegg et al. 2019). Prior to the era of Gaia, underpinning essentially all studies of the global kinematic structure of the MW stellar halo, as well as estimates of 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, we can confront the ways in which these assumptions are violated and use this information to understand our Galaxy on a deeper level.

One major source of disequilibrium in the MW is its most massive satellite, the Large Magellanic Cloud (LMC). The classical picture of the LMC is as a relatively low mass ($\sim {10}^{10}{M}_{\odot }$) satellite orbiting the MW on a T ∼ 2 Gyr orbit (e.g., Avner & King 1967; Hunter & Toomre 1969; Murai & Fujimoto 1980; Lin & Lynden-Bell 1982; Lin et al. 1995; Bekki & Chiba 2005; Mastropietro et al. 2005; Connors et al. 2006). However, proper motion (PM) measurements of the LMC using the Hubble Space Telescope (HST) revealed that the total velocity of the LMC is much higher than previously measured (v ∼ 320 $\mathrm{km}\,{{\rm{s}}}^{-1}$; 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 backward orbital integrations (Besla et al. 2007; Kallivayalil et al. 2013) and statistical predictions from cosmological simulations (Boylan-Kolchin et al. 2011; Busha et al. 2011; González et al. 2013; Patel et al. 2017). In addition, there is mounting evidence that the LMC is more massive than previously thought ($\sim {10}^{11}{M}_{\odot }$), including arguments based on models of the Magellanic system (e.g., Besla et al. 2010, 2012; Pardy et al. 2018); abundance matching (e.g., Behroozi et al. 2010; Guo et al. 2010; Moster et al. 2010, 2013); the measured rotation curve of the LMC (van der Marel & Kallivayalil 2014); the presence of satellites around the LMC, including the Small Magellanic Cloud (e.g., Kallivayalil et al. 2018; Erkal & Belokurov 2020; Patel et al. 2020; Pardy et al. 2020); the timing argument (Peñarrubia et al. 2016); and perturbations in the Orphan Stream (Erkal et al. 2019; Koposov et al. 2019).

While concerns about the LMC's influence on dynamics in the MW were raised as early as Avner & King (1967), the revised picture of the LMC as a massive ($\sim {10}^{11}{M}_{\odot }$) satellite approaching the MW for the first time is increasingly worrisome for estimates of the MW gravitational potential that neglect the LMC, as the LMC mass is a significant fraction of the MW halo mass. Several studies of simulations of the LMC's infall demonstrate that a massive LMC invalidates the assumption of an inertial Galactocentric reference frame, as the center of mass (COM) can be substantially displaced (by as much as 30 kpc; Gómez et al. 2015) from the center of the Galaxy, resulting in net motion of the halo with respect to the MW disk. This net COM motion is predicted to be ∼40 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (GC19, Erkal et al. 2019; Petersen & Peñarrubia 2020); Gómez et al. (2015) find that the net motion could be as high as 75 $\mathrm{km}\,{{\rm{s}}}^{-1}$. Erkal et al. (2020) show that ignoring the influence of the LMC leads to systematic overestimates (as high as 50%) of the MW mass when using equilibrium models. Therefore, when considering the motions of tracer populations that we use to study the MW DM halo, we must account for the influence of the LMC in our models.

In addition to perturbations as a result of the COM motion of the LMC, MW halo tracers are also predicted to be perturbed by the LMC-induced DM wake (Garavito-Camargo et al. 2019, hereafter GC19). In the ΛCDM cosmological paradigm, host halos are predicted to respond to the infall of satellites; this response can be thought of as a gravitational or density wake. One component of this wake arises owing to local interactions of particles with the satellite. The satellite transfers kinetic energy to nearby resonant particles, creating an overdensity trailing in its orbit and causing an effective drag force on the satellite (i.e., dynamical friction; e.g., Chandrasekhar 1943; White 1983; Tremaine & Weinberg 1984). GC19 refer to this component of the wake as the Transient response, as it is expected to weaken over time. In addition to the Transient response, there is also a global response in the DM halo, resulting in large-scale over- and underdensities in the DM halo (e.g., Weinberg 1989) that can potentially even excite structure in the disk (Weinberg 1998; Weinberg & Blitz 2006; Laporte et al. 2018a). GC19 refer to this component of the wake as the Collective response. For the benefit of the reader, we have included these definitions in Table 1. Using detailed N-body simulations, GC19 demonstrated that the density wake induced by the infall of the LMC gives rise to distinct, correlated kinematic patterns in the MW stellar halo.

Table 1.  Useful Definitions

	LMC-Induced Wake Components from GC19
Transient response	The overdensity trailing the LMC in its orbit, which arises due to local scattering. This component can be
	thought of in terms of classical dynamical friction.
Collective response	Refers to both the overdensity in the north (which arises owing 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 the
	MW halo to the LMC's infall.
	Spherical Harmonics
θ	Colatitudinal angle, measured downward from the z-axis; $\cos \theta =z/R$
ϕ	Azimuthal angle, angle in the x–y plane measured from the x-axis; $\tan \phi =y/x$
ℓ	Order of spherical harmonic ${Y}_{{\ell }m}$
m	Degree of spherical harmonic ${Y}_{{\ell }m}$
${a}_{{\ell }m}$	Spherical harmonic coefficient for mode ℓ m
${\varphi }_{{\ell }m}$	Phase of spherical harmonic coefficient ${a}_{{\ell }m}$
Cℓ	Average power in mode of order ${\ell };$ total power is $(2{\ell }+1)\times {C}_{{\ell }}$
Zonal spherical harmonics	Spherical harmonics with ℓ = 0; rotational symmetry about z-axis
Sectorial spherical harmonics	Spherical harmonics with ${\ell }=| m|$

Note. We use the function arctan2 in NumPy (Oliphant 2006) to compute azimuthal angle ϕ and phase φ, such that both angles take values over the range [−π, π].

Download table as:  ASCIITypeset image

GC19 explored these wake signatures in the context of two MW halo models, one isotropic halo and one radially varying, radially anisotropic halo. They find that while there are similarities in the wake morphology for the two models, there are also key differences: the Transient response is much stronger for the model with the radially anisotropic halo, whereas the Collective response is stronger for the isotropic halo. Therefore, understanding how the velocity anisotropy $\beta =1-{\sigma }_{T}^{2}/{\sigma }_{R}^{2}$ behaves in the MW is important for the predicted morphology of the LMC-induced DM wake. Until relatively recently, our knowledge of the motions of halo tracers was limited to one component of motion, the line-of-sight (LOS) velocity; given this major observational constraint, it was necessary to make assumptions about the tangential motions of stars, and isotropy (β = 0) was the most common assumption. However, simulations predict that β should become increasingly radially biased as a function of radius (see, e.g., Rashkov et al. 2013; Loebman et al. 2018), and β in the solar neighborhood is radially biased (β ∼ 0.5 − 0.7; Smith et al. 2009; Bond et al. 2010).

We now have full phase-space information for distant halo tracers, from the Gaia mission and HST PM studies, and can measure β outside the solar neighborhood directly. Estimates of β outside of the solar neighborhood have generally found radially biased β, using GCs (Sohn et al. 2018; Watkins et al. 2019) and halo stars (Bird et al. 2019; Cunningham et al. 2019b; Lancaster et al. 2019b).

Detecting the halo response to the LMC-induced DM wake would be an exciting advancement in testing our assumptions about the properties of DM, as well as providing key constraints on the potential of the MW and the mass and orbital history of the LMC. However, the GC19 simulations give predictions for the response in the context of smooth MW DM and stellar halos. In reality, the MW stellar halo contains 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 recent review) and stellar clouds (e.g., Newberg et al. 2002; Rocha-Pinto et al. 2004; Jurić et al. 2008; Li et al. 2016). In addition, using a sample of MW halo main-sequence turnoff stars from the HALO7D survey (Cunningham et al. 2019a), Cunningham et al. (2019b) observed that the estimated parameters of the velocity ellipsoid (i.e., $\langle {v}_{\phi }\rangle ,{\sigma }_{\phi },{\sigma }_{R},{\sigma }_{\theta }$) were different in the different survey fields; these differing estimates could be interpreted as evidence that the halo is not phase mixed over the survey range ($\langle r\rangle =23$ kpc). They also showed maps of the halo velocity anisotropy β in two halos from the Latte suite of FIRE-2 simulations (introduced in Wetzel et al. 2016), finding that the anisotropy can vary over the range $[-1,1]$ across the sky. Some of the variation in the β estimates appeared to correlate with stellar overdensities in the halos, indicating that galactic substructure is at least in part responsible for the different velocity distributions. While some substructure in the halo can be clearly identified as overdensities in phase space and removed from analysis, the presence of velocity substructure in the halo could complicate attempts to detect signatures of the LMC-induced DM wake. For example, Belokurov et al. (2019) recently argued that the Pisces Overdensity (Sesar et al. 2007; Watkins et al. 2009; Nie et al. 2015) might be stars in the wake trailing the LMC in its orbit, because of their net negative radial velocities. However, it remains difficult to conclusively argue this scenario given that these stars could also be in Galactic substructure (or, perhaps, stars that are in substructure and have been perturbed by the DM wake).

In summary, there is substantial observational evidence that MW stellar halo is in disequilibrium, on average radially anisotropic, and rife with un-phase-mixed substructure, in clear violation of the three central assumptions of equilibrium models. The velocity field of the halo contains information about the potential of the MW, the dwarf galaxies that were consumed as the MW assembled its mass, and the properties of its current most massive perturber, the LMC. However, separating out the different origins of the features in the MW halo velocity field remains a formidable challenge. One way forward is to consider the spatial scale of perturbations: we expect substructure to cause velocity variation on relatively small spatial scales, as opposed to the large-scale perturbations from the LMC-induced DM wake. Therefore, to disentangle these effects, we seek a quantitative description of the kinematic structure of the halo that incorporates variation on different spatial scales. Spherical harmonic expansion is a natural tool to address this problem.

While ideally we would embark on a full basis function expansion (BFE) of the phase-space structure of the halo, in this work we focus on the spherical harmonic expansion of the three components of motion in spherical coordinates, over different distance ranges in the halo (as a complement to this work, N. Garavito-Camargo et al. (2020, in preparation), will present full BFEs of the spatial distributions for these simulations). GC19 explored the density structure and the properties of the velocity dispersions in addition to the mean velocities; we choose to focus on the mean velocities here because of the challenges of estimating densities (i.e., deeply understanding completeness and survey selection functions) and the fact that estimates of the mean of a distribution require fewer tracers than dispersion estimates.

This paper is organized as follows. In Section 2, we present a brief overview of spherical harmonic expansion and define the notation used in the rest of the paper. We then show the results of using spherical harmonic expansion on the velocity fields from the GC19 simulations of the LMC's infall into the MW in Section 3. We demonstrate how the different components of the wake are described in terms of spherical harmonics. In Section 4, we investigate how Galactic substructure might complicate our ability to measure perturbations to the velocity field as a result of the LMC-induced DM wake, by studying the Bullock & Johnston (2005) purely accreted stellar halos. In Section 5, we use two models of the Sagittarius stream to estimate how the MW's most massive stream might influence the angular power spectrum of the MW halo velocity field and interfere with signatures from the wake. We summarize our conclusions in Section 6.

We seek to describe the variation on different spatial scales in halo velocity fields by using spherical harmonic expansion. In this section, we define the notation we use throughout the paper for spherical harmonics. As a reference, we have also included many of these definitions in Table 1. Laplace's spherical harmonics of order ℓ and degree m are defined as

$$\begin{eqnarray}&&{Y}_{{\ell }}^{m}(\theta ,\phi )=\sqrt{\displaystyle \frac{2{\ell }+1}{4\pi }\displaystyle \frac{({\ell }-m)!}{(l+m)!}}{P}_{{\ell }}^{m}(\cos \theta ){e}^{{im}\phi },\end{eqnarray} \tag{ 1 }$$

where θ is the colatitudinal angle (i.e., the polar angle measured downward from the north pole), ϕ is the azimuthal angle (i.e., the angle in the x − y plane measured from the x-axis), and ${P}_{{\ell }}^{m}$ are the associated Legendre polynomials.

Spherical harmonics comprise an orthogonal basis for any function f(θ, ϕ) defined on the surface of the sphere:

$$\begin{eqnarray}&&f(\theta ,\phi )=\displaystyle \sum _{{\ell }=0}^{\infty }\displaystyle \sum _{m=-{\ell }}^{m={\ell }}{a}_{{\ell }m}{Y}_{{\ell }}^{m}(\theta ,\phi ),\end{eqnarray} \tag{ 2 }$$

where ${a}_{{\ell }m}$ are the spherical harmonic coefficients, given by

$$\begin{eqnarray}&&{a}_{{\ell }m}={\int }_{{\rm{\Omega }}}f(\theta ,\phi ){Y}_{{\ell }}^{m* }(\theta ,\phi )d{\rm{\Omega }}.\end{eqnarray} \tag{ 3 }$$

When the spherical harmonics are complex, the coefficients are also complex; we define the phase ${\varphi }_{{\ell }m}$ of a spherical harmonic coefficient as

$$\begin{eqnarray}&&{\varphi }_{{\ell }m}={\tan }^{-1}\left(\displaystyle \frac{\mathrm{Im}[{a}_{{\ell }m}]}{\mathrm{Re}[{a}_{{\ell }m}]}\right),\end{eqnarray} \tag{ 4 }$$

where we use the function arctan2 implemented in NumPy (Oliphant 2006) to compute the inverse tangent, taking into account the quadrant in which ${a}_{{\ell }m}$ lies in the complex plane.

The angular power spectrum C_ℓ can be computed from the spherical harmonic coefficients ${a}_{{\ell }m}$:

$$\begin{eqnarray}&&{C}_{{\ell }}=\displaystyle \frac{1}{2{\ell }+1}\displaystyle \sum _{m}| {a}_{{\ell }m}{| }^{2}.\end{eqnarray} \tag{ 5 }$$

The total power in a given order ℓ is thus $(2{\ell }+1)\times {C}_{{\ell }}$, as there are $(2{\ell }+1)$ total m values for a given ℓ value. Therefore, in this paper, power spectra will always have the quantity $(2{\ell }+1)\times {C}_{{\ell }}$ (in units of ($\mathrm{km}\,{{\rm{s}}}^{-1}$)²) plotted on the y-axis, as we are expanding the velocity field.

In this work we use the Python package healpy (Zonca et al. 2019),⁸ based on the Healpix scheme (Górski et al. 2005) to perform all of our analysis relating to spherical harmonics. All maps are made using the healpy plotting routine mollview, power spectra and spherical harmonic coefficients are computed using the function anafast, and synthetic maps are generated using synfast.

While healpy works with the spherical harmonics in complex form, the real spherical harmonics are defined as

$$\begin{eqnarray}{Y}_{{\ell }m}(\theta ,\phi )=\left\{\begin{array}{ll}{P}_{{\ell }m}(\cos \theta )\cos (m\phi ) & m\geqslant 0\\ {P}_{{\ell }| m| }(\cos \theta )\sin (| m| \phi ) & m\lt 0\end{array}\right..\end{eqnarray} \tag{ 6 }$$

For illustrative purposes, we show the real spherical harmonics from ℓ = 0 to ℓ = 4 in Figure 1. For a given ${Y}_{{\ell }m}$, the degree m corresponds to the number of waves along a line of constant latitude. The order ℓ, in conjunction with degree m, determines how many times zero is crossed along a line of constant longitude: there are ${\ell }-| m|$ zero crossings along a meridian. In the case of ${\ell }=| m|$, there are no zero crossings along the meridian (outer column in each row of Figure 1), and a total of ℓ complete waves along the equator; these modes are referred to as the sectorial spherical harmonics. When m = 0, there are ℓ zero crossings along the meridian (central column of Figure 1), and no change in amplitude with longitude; these are known as the zonal spherical harmonics and have symmetry about the z-axis. Modes with $m\lt 0$ are of sine type, while modes with m > 0 are of cosine type; as demonstrated by Figure 1, for the real spherical harmonics, changing the sign of m results in a 90° rotation about the z-axis.

Zoom In Zoom Out Reset image size

Spherical harmonic expansion and angular power spectra have many applications in astrophysics, most famously in studies of the cosmic microwave background (e.g., Planck Collaboration et al. 2018, 2019, and references therein). While spherical harmonics are commonly used to describe the angular dependence in full BFEs of the potential of DM halos (e.g., Hernquist & Ostriker 1992; Weinberg 1996, 1999; Lowing et al. 2011), they have not generally been used to describe velocity fields. In the following sections, we discuss spherical harmonic expansion of the velocity fields of several different types of simulations.

In this section, we perform spherical harmonic expansion of the velocity fields on the high-resolution, N-body simulations of the LMC's infall into the MW from GC19. The kinematic patterns (for mean velocities in all components, as well as the densities and velocity dispersions) are discussed in detail in GC19. Here, we discuss how spherical harmonic expansion of the mean velocities can be used to characterize the MW halo response to the LMC's infall.

This section is organized as follows. We summarize key GC19 simulation properties in Section 3.1. In Section 3.2, we discuss the spherical harmonic expansion of the velocity maps at 45 kpc in detail. In Section 3.3, we discuss the radial evolution of the power spectrum, for both the isotropic and radially anisotropic MW models. The dependence of the power spectra on the LMC mass is explored in Section 3.4.

3.1.  GC19 Simulation Details

3.2. The Velocity Field near ${R}_{\mathrm{LMC}}$

We first discuss the velocity field from the GC19 fiducial LMC model (${M}_{\mathrm{LMC}}=1.8\times {10}^{11}{M}_{\odot }$) with the radially anisotropic 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 expect the velocity maps over this radial range to be sensitive to the Transient response, given that the LMC passed through very recently, and the Transient response arises owing to local scattering. The mean velocity maps in three components of spherical motion (with respect to the Galactic center) are shown in the top panels of Figure 2.

Zoom In Zoom Out Reset image size

At z = 0 in these simulations, the LMC is located at 50 kpc, and its angular position indicated by the star symbol in the top panels of Figure 2. The color of the star symbol indicates the sign of the motion of the LMC in each component of motion: ${({v}_{R},{v}_{\theta },{v}_{\phi })}_{\mathrm{LMC}}=(99,-345,-46)\,\mathrm{km}\,{{\rm{s}}}^{-1}$. At 45 kpc, the motions of the stars in the Transient response near the LMC's position trace the COM motion of the LMC. The net motion outward in v_R and the converging motions v_θ and v_ϕ near the LMC's position result from stars in the Transient response, accelerating toward the LMC. In addition, stars at 45 kpc are being accelerated toward the overdensity in 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 ${{\ell }}_{\max }=5$ expansion of each component of velocity. Because the kinematic variation induced by the LMC occurs on large scales, the dominant features in the velocity maps are effectively captured by a low-order spherical harmonic expansion. The bottom panels of Figure 2 show the magnitudes of the spherical harmonic expansion coefficients ($| {a}_{{\ell }m}|$), color-coded by m value. In the radial velocity map (left panels), the dominant term is the ${\ell }=2,m=\pm 1$ mode; this mode captures the outward radial motion in the upper left quadrant and the lower right quadrant (near the LMC), as well as the inward radial motion in the lower left quadrant and upper right quadrant. In the v_θ maps, the monopole term (ℓ = m = 0) is dominant, reflecting the net upward motion of the halo with respect to the disk. In the v_ϕ maps, the ${\ell }=2,m=\pm 2$ mode dominates; the sectorial ℓ = 2 mode captures the converging motions of stars in the Transient response near the location of the LMC.⁹ The fact that there is no power at ℓ = 0 in v_ϕ is indicative of the fact that the GC19 MW models have no net rotation. If the MW does have any net rotation (which has been observed, but not with high statistical significance; Deason et al. 2017), this would result in power at ℓ = 0 in v_ϕ, which would not interfere with the predicted wake signature.

We note that we have not included error bars on the spherical harmonic coefficients in Figure 2, nor do we include error bars on our estimates of the power spectra in subsequent figures. This is because the simulations are very high resolution, so the statistical errors on the coefficients are very small; the dominant sources of uncertainty here are in the models, not in the noise in the simulations.

3.3. Evolution with Distance

Figure 3 shows the velocity maps in the three components of spherical motion for this simulation at 45, 70, and 100 kpc, in the Galactocentric frame. Radial velocity maps are shown in the left panels, polar motion v_θ is plotted in the middle panels, and azimuthal motion is shown in the right panels. As a result of the Collective response, there is a radial velocity dipole that increases in strength as a function of 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 terms of v_z: in these simulations, while the net motion in the plane of the disk is fairly stable ($\langle {v}_{x}\rangle =\langle {v}_{y}\rangle \sim 0$ $\mathrm{km}\,{{\rm{s}}}^{-1}$  over all radii), there is net upward motion in the halo ($\langle {v}_{z}\rangle \gt 0$), which increases as a function of distance from the MW's disk. Erkal et al. (2020) show that the MW globular clusters and dwarf satellites also show net motion in v_z (and no net motion in v_x, v_y); however, they note the caveat that these tracers may not be phase mixed in the MW potential. We note that the velocity shifts in these simulations result from two sources: the DM overdensity in the north (which gets stronger with distance, because the LMC spends more time at larger radii), and the net acceleration of the MW disk toward the LMC (e.g., Petersen & Peñarrubia 2020). Disentangling the relative contributions to the overall velocity shift from these two sources is beyond the scope of this work.

Zoom In Zoom Out Reset image size

The power spectra for these maps are shown by the bold lines in Figure 4. The kinematic patterns are concisely summarized by the angular power spectra. The radial velocity dipole that increases in magnitude with Galactocentric distance is reflected by the increasing power in the ℓ = 1 modes; the increasing mean polar velocity as a function of distance is captured by the increasing power in ℓ = 0 (i.e., the monopole). The power in v_ϕ is strongest at 45 kpc, where the stars in the Transient response are closest to the present-day position of the LMC and are accelerated by its COM motion.

Zoom In Zoom Out Reset image size

The faded lines in Figure 4 are the resulting power spectra for the GC19 simulation using an isotropic MW halo. As discussed in GC19, the Collective response is stronger in the isotropic simulations. This is reflected in the resulting power spectra. In the radial velocity power spectrum, at large radii, the magnitude of the velocity dipole is larger for the isotropic halo, as is the magnitude of the ${v}_{\theta }$ monopole, reflecting the stronger 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}_{\mathrm{LMC},\mathrm{vir}}=0.8,1.0,1.8,2.5\times {10}^{11}{M}_{\odot }$. Figure 5 shows the resulting power spectra for the mean velocities in the three spherical components of motion at 45, 70, and 100 kpc (top, middle, and bottom panels, respectively). The different line styles in each figure represent the different LMC masses.

Zoom In Zoom Out Reset image size

While the shapes of these power spectra at a given distance for a given component of motion are overall very similar to one another (highlighted by the bottom panel, which shows the power spectra on a logarithmic scale), the total power at a given ℓ value clearly trends with LMC mass. In GC19, by quantifying the "strength" of the wake as the magnitude of the density fluctuations, they find that the strength of the wake is comparable at 45 kpc for all LMC masses (see their Figure 25; though this is for the isotropic MW model, which has a very weak Transient response). Based on the top panels of Figure 5, we can see that the power at different ℓ values increases strongly with LMC mass even at 45 kpc. We emphasize that the y-axis labels are different in each panel, to emphasize the effect of changing the LMC mass; we note that the peak of the power spectrum is largest in v_R at all distances, and the v_ϕ power spectrum has the least amount of power at all distances.

The sensitivity of these signals to the LMC mass shown in Figure 5 emphasizes the significance of the paradigm shift from a ${10}^{10}{M}_{\odot }$ mass LMC to a favored $\sim {10}^{11}{M}_{\odot }$ mass LMC. If the LMC were only ${10}^{10}{M}_{\odot }$, this would be a factor of eight less massive than the least massive LMC simulated by GC19; based on the power spectra in Figure 5, we can see that the signatures of the DM wake in this scenario would be very weak. Because the LMC is favored to be ∼10% of the MW mass, as opposed to ∼1% (like the other MW satellites), we cannot ignore its gravitational influence on MW halo tracers.

3.5. Summary

The GC19 simulations model the MW DM halo (and, as result, their stellar halo) as smooth; however, we know that the MW stellar halo is structured. In this section, we investigate how Galactic substructure might complicate our ability to characterize the LMC-induced DM wake using the spherical harmonic expansion of the velocity field, using the Bullock & Johnston (2005) suite of simulations of purely accreted stellar halos (hereafter BJ05).

The publicly available BJ05 simulations are N-body simulations of accreted dwarf galaxies onto an MW-like parent galaxy. The full suite of simulations consists of 1515 individual accretion events, with a variety of masses, orbital parameters, and accretion times, that together make up the 11 traditional BJ05 halos. Each disrupted satellite is modeled with 10⁵ DM particles. The parent galaxy is represented by a time-evolving potential with disk, halo, and bulge components. Johnston et al. (2008) explore in depth how the observable properties of the substructure in these simulations are related to the properties of their satellite progenitors. Here, we explore the links between a galaxy's accretion history, the resulting velocity maps, and the angular power spectra of the velocity field. We note that because there are no DM particles in the parent galaxy halo, there are no density wakes induced in the BJ05 simulations. In this section we are only concerned with spatially varying mean velocities arising from debris from accreted dwarfs.

Specifically, we consider the six halos with "artificially constructed" accretion histories discussed in Johnston et al. (2008). While the standard 11 BJ05 halos have accretion histories from merger trees constructed in the ΛCDM cosmological context, the artificially constructed halos contain debris from accretion events selected from the full BJ05 library that have the desired properties. While all six halos end up with a total luminosity $L\sim {10}^{9}{L}_{\odot }$, they assemble their halos very differently. The six artificial halos are as follows:

• 1.  
  Circular (radial) orbits: halo assembled only from accretion events with ${J}_{\mathrm{sat}}/{J}_{\mathrm{circ}}\gt 0.75$ (${J}_{\mathrm{sat}}/{J}_{\mathrm{circ}}\lt 0.2$).
• 2.  
  Recent (early) accretion: halo assembled only from accretion events with t_acc < 8 Gyr (t_acc > 11 Gyr).
• 3.  
  High L_sat (low L_sat): halo assembled only from accretion events with ${L}_{\mathrm{sat}}\gt {10}^{7}{L}_{\odot }$ (${L}_{\mathrm{sat}}\lt {10}^{7}{L}_{\odot }$).

Figure 6 shows velocity maps of the six artificially constructed BJ05 halos for stars in the distance range of 30–50 kpc (with respect to the Galactic center), projected in Mollweide coordinates. As in previous figures, left panels are radial velocity (v_R) maps, middle panels are polar velocity (v_θ) maps, and right panels are azimuthal velocity (v_ϕ) maps. We emphasize that the color bars for these maps range from $[-200,200]\,\mathrm{km}\,{{\rm{s}}}^{-1}$, a much larger range than shown for the GC19 simulations; the amplitudes of the velocity fluctuations in these maps are much greater than those due to the LMC-induced DM wake as seen in GC19. As a result, the signatures from substructure are difficult to compare with the wake signatures using the maps alone; the angular power spectra corresponding to these maps, along with the GC19 power spectra, are plotted in Figure 7. We also note that due to the resolution of the BJ05 simulations, there are pixels that contain no star particles; these pixels are assigned to have $\langle v\rangle =0\,\mathrm{km}\,{{\rm{s}}}^{-1}$, consistent with the assumptions of an equilibrium model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a result of the different (and extreme) accretion histories experienced by these halos, their velocity maps look very different. Unsurprisingly, the halo that experienced only early accretion has no features in any components of motion in its velocity maps (bottom panels of Figure 6), given that all of its accreted material has had sufficient time to become phase mixed. In contrast, the halo that accreted massive, high-luminosity satellites has large-scale features in all components of motion (second row in Figure 6). It is also worth noting that this halo accreted ∼35% of its mass within the past 8 Gyr (see Figure 7 of Johnston et al. 2008); it has therefore experienced recent accretion in addition to massive accretion. The halo that experienced mostly circular accretion events (top row of Figure 6) has low levels of variation in its radial velocity map, with many thin features of nearly constant (and approximately 0 $\mathrm{km}\,{{\rm{s}}}^{-1}$) radial velocity, corresponding to streams. These streams have more energy in tangential than radial motion: they appear as bands of nearly constant but high velocity in the v_θ and v_ϕ maps. The halo that experienced only recent accretion (third row of maps in Figure 6) contains both debris from massive accretion events (as seen by large patches of stars at common velocities) and 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 features in any components of motion; however, it is important to keep in mind that these halos also have had relatively quiescent recent accretion histories. Both halos had assembled $\sim 95 \%$ of their mass 8 Gyr ago (see again Figure 7 of Johnston et al. 2008); therefore, any features in these maps are due to relatively recent, low-mass events. A radial stream appears 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 velocity maps are plotted in Figure 7. The halos with the most power at all ℓ values are the halo built from recent accretion and the halo built from high-luminosity satellites. The halo built from circular accretion events also has high power in v_θ and v_ϕ. The thin streams in the low-L_sat halo also result in substantial power (> 100 $\mathrm{km}\,{{\rm{s}}}^{-1}$) over many ℓ in the three components of motion, though not as much power as the recently accreted, high-L_sat, and circular halos.

The halos that experienced only radial accretion events and only early accretion events both have nearly featureless maps in v_θ and v_ϕ; while with few to no features in radial velocity, we note that the v_R maps appear noticeably noisier than the other components of motion. This is a result of the fact that these halos have radially biased velocity anisotropy: their radial velocity dispersions are much greater than their tangential velocity dispersions. The resulting radial velocity maps have greater fluctuations from pixel to pixel than their tangential velocity counterparts. In addition, because there are many fewer particles in these simulations than in GC19 (even though we are looking at a larger radial range in BJ05), the pixel-to-pixel variation is higher for the BJ05 maps than the GC19 maps. This results in some power in the radial velocity power spectrum, albeit with less power than the other halos, and with hardly any power in the tangential velocity components.

The purple shaded region in Figure 7 shows the range of power at 45 kpc for the GC19 simulations, from the lowest-mass to the highest-mass LMC. The power from the halo formed through recent accretion and the high-L_sat halo is much greater than the power from the LMC-induced DM wake at nearly all ℓ in all components of motion. Even the circular and low-L_sat halos can have comparable signals to the wake in the tangential components of motion. However, the shape of the power spectrum is very different for Galactic substructure than for the LMC-induced DM wake. The power spectra are characterized by a sawtooth pattern, with peaks at odd ℓ values for v_R and v_θ and peaks at even ℓ values for v_ϕ. This sawtooth pattern is indicative of the fact that spherical harmonics are not the ideal basis for the velocities of stars in substructure; we explore why the power spectra have these features in the Appendix.

Based on the power spectra plotted in Figure 7, if the MW stellar halo is dominated by debris from recent, massive accretion events, the signal from the LMC-induced DM wake in the velocity field would be overwhelmed by the Galactic substructure. This signal is from debris that has not yet phase mixed; based on the velocity maps, we can see that this substructure is clearly visible as overdensities in phase space, not just in velocity space. Therefore, one could take advantage of the fact that many of these features would be clearly identifiable observationally as overdensities and could be removed from the analysis relatively easily.

In addition, while it is likely that debris from an early, massive accretion event dominates the inner halo (i.e., the Gaia Sausage/Gaia-Enceladus, ∼10 Gyr ago; e.g., Belokurov et al. 2018; Helmi et al. 2018), the current consensus is that the MW has had a fairly quiescent recent accretion history. This consensus has emerged based on numerous studies, including studies of the structure and kinematics of stars of the Galactic disk plane (e.g., Gilmore et al. 2002; Hammer et al. 2007; Ruchti et al. 2015), the steep stellar density profile beyond ∼25 kpc in the halo (e.g., Deason et al. 2013; Pillepich et al. 2014; Deason et al. 2018; Fukushima et al. 2018, 2019; Thomas et al. 2018; Starkenburg et al. 2019), and the amount of substructure in the halo relative to predictions from simulations (e.g., Lancaster et al. 2019a). An alternative scenario is one in which the MW has experienced more recent low-luminosity or radial accretion events with debris that is harder to find observationally; for example, Donlon et al. (2019) suggest that a recent (∼2 Gyr), radial merger (with $M\sim {10}^{9}{M}_{\odot }$) that mixes efficiently could explain the Virgo Overdensity (Vivas et al. 2001). Regardless, the MW's accretion history is not believed to be dominated by recent massive accretion.

The major exceptions to the picture of the MW as having a quiescent (massive) recent merger history are the relatively recent accretion of Sagittarius (Sgr; ∼6 Gyr ago) and the LMC (∼2 Gyr ago). In the following section, we explore how the presence of debris from Sgr might impact the spherical harmonic expansion of the halo velocity field.

In Section 3, we discussed in detail the spherical harmonic expansion of the kinematic variation that arises because of the infall of the LMC. In Section 4, we saw that recent, massive accretion events can cause kinematic variation on large scales, which in turn can result in substantial power over many ℓ values in the power spectrum. In the MW, the debris from the most recent, massive accretion event (aside from the LMC) is found in the Sgr stream (Ibata et al. 2001). The progenitor of Sgr was a relatively massive, luminous satellite (L ∼ 10⁸, $M\gt {10}^{9}{M}_{\odot };$ e.g., Peñarrubia et al. 2010; Niederste-Ostholt et al. 2012; Deason et al. 2019. Gibbons et al. (2017) suggest an even higher mass M > 6 × 10¹⁰). Debris from the Sgr accretion event is found all over the sky over Galactocentric distances ranging from ∼15 to ∼130 kpc (e.g., Majewski et al. 2003; Belokurov et al. 2014; Hernitschek et al. 2017; Sesar et al. 2017). In this section, we investigate how the presence of Sgr stars might obscure the signal from the LMC-induced DM wake in the velocity field.

We consider two different models for the Sgr stream. The first Sgr model we consider is a fit of the Sagittarius stream in the presence of the LMC, which we dub the Erkal model. This model uses the same stream-fitting machinery as in Erkal et al. (2019), which accounts for the reflex motion of the MW due to the LMC. This technique rapidly generates streams using the modified Lagrange Cloud stripping technique from Gibbons et al. (2014). For this model, we fit the radial velocity and distances from Belokurov et al. (2014) and on-sky positions from Belokurov et al. (2006) and Koposov et al. (2012) for the bright stream.

Motivated by the results of Law & Majewski (2010), we model Sgr as a $2.5\times {10}^{8}{M}_{\odot }$ Plummer sphere with a scale radius of 0.85 kpc. The progenitor is rewound for 5 Gyr in the combined presence of the MW and LMC and then disrupted to the present to form the Sgr stream. For the MW potential, we take the triaxial NFW generalization from Bowden et al. (2013), which allows for different inner 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 necessarily aligned with the galactic Cartesian coordinates. We also include a similar disk and bulge to the MWPotential2014 from Bovy (2015): a Miyamoto–Nagai disk (Miyamoto & Nagai 1975) with a mass of $6.8\times {10}^{10}{M}_{\odot }$ with a scale radius of 3 kpc and a scale height of 0.28 kpc, and a Hernquist bulge (Hernquist 1990) with a mass of $5\times {10}^{9}{M}_{\odot }$ and a scale radius of 0.5 kpc. We use the dynamical friction prescription of Jethwa et al. (2016) for the dynamical friction both from the MW on the LMC and from the MW on Sgr. The distance, radial velocity, and PMs of the Sgr are left as free parameters with priors set by observations (McConnachie 2012). For the LMC, we give it a fixed position and velocity based on its mean observed distance (Pietrzyński et al. 2013), radial velocity (van der Marel et al. 2002), and PM (Kallivayalil et al. 2013). We model the LMC as a Hernquist profile (Hernquist 1990) with a scale radius of 25 kpc and a free mass with a uniform prior from $0\,\mathrm{to}\,3\times {10}^{11}{M}_{\odot }$. In order to account for the fact that Sgr was initially more massive and had a substantial DM component that would have experienced more dynamical friction, we have an additional free parameter that increases the mass of Sgr by λ_DF when computing its dynamical friction. This has a uniform prior between 0 and 20. Thus, altogether we have 15 free parameters: the mass and scale radius of the NFW profile (${M}_{\mathrm{NFW}},{r}_{s\mathrm{NFW}}$), an inner and outer minor- and intermediate-axis flattening (${q}_{0},{p}_{0},{q}_{\infty },{p}_{\infty }$), three angles to describe the rotation of the triaxial halo, the mass of the LMC (M_LMC), the mass multiplier λ_DF, and finally the PMs, radial velocity, and distance of the Sgr progenitor at the present day.

We use the MCMC package from Foreman-Mackey et al. (2013) to estimate the parameters of Sgr. We use 100 walkers for 2000 steps with a 1000-step burn-in. The best-fit parameters require an LMC mass of $2.0\times {10}^{11}{M}_{\odot }$, flattenings of ${q}_{0}=0.68,{p}_{0}=0.87,{q}_{\infty }=0.81,{p}_{\infty }=0.94$. The mass multiplier λ_DF = 9.5, suggesting that fitting the Sgr stream requires more dynamical friction than the low mass we have assumed would provide. The best-fit MW mass is $6.76\times {10}^{11}{M}_{\odot }$ with a scale radius of 15.3 kpc. Although this MW mass is relatively modest, we note that the scale radius is also quite small. Despite the flexibility of this model, we note that it does not perfectly match the distance, but most importantly for this work, it gives a good match to the radial velocity across the sky. A comparison of the model with radial velocities from Belokurov et al. (2014) is shown in Figure 8. The positions of the stars for this Sgr model in the x–z plane, color-coded by heliocentric LOS velocity, are shown in the top left panel of Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The second Sgr model we consider is the publicly available model from Dierickx & Loeb (2017, hereafter DL17).¹⁰ The x–z positions of stars from this model, again color-coded by heliocentric LOS velocity, are shown in the right panel of Figure 9. In DL17, they first utilize a semianalytic approach to derive initial conditions for the Sgr progenitor, by integrating the equations of motion forward in time over 7–8 Gyr and comparing the resulting position and velocity vector to the observed properties of the Sgr remnant. They assume virial masses ${M}_{\mathrm{Sgr}}={10}^{10}{M}_{\odot }$ and ${M}_{\mathrm{MW}}={10}^{12}{M}_{\odot }$. To then model the disruption of Sgr, they run an N-body simulation using the derived initial conditions from the semianalytic approach, modeling both a live MW and Sgr. While this model does reproduce many of the features of the Sgr stream, including the positions of stars observed in the Two Micron All Sky Survey (Majewski et al. 2004) and the Sloan Digital Sky Survey (SDSS; Belokurov et al. 2014) and the large apocentric distances observed in Sesar et al. (2017), we emphasize that the N-body simulation is not tuned to fit the observations of the Sgr stream. As a result, certain properties of the stream (e.g., the LOS velocities along the leading arm) are not well matched by the data.

To investigate how stars from Sgr might impact the power spectrum of the MW halo's velocity field, we overlay the two models for the Sgr stream onto the fiducial anisotropic GC19 simulation (with ${M}_{\mathrm{LMC}}=1.8\times {10}^{11}{M}_{\odot }$). To combine the two independent simulations, we assign the total Sgr stellar mass to be 10% of the total stellar mass of the GC19 halo. This ratio is consistent with current estimates of the total stellar mass of the MW ($\sim {10}^{9}{M}_{\odot };$ e.g., Deason et al. 2019; Mackereth & Bovy 2020) and Sgr (${M}_{\mathrm{Sgr},* }\sim {10}^{8}{M}_{\odot };$ e.g., Niederste-Ostholt et al. 2012; Deason et al. 2019).

The resulting velocity maps at 45 kpc of the two Sgr models overlaid on the GC19 simulations are shown in the middle panels (for the Erkal model) and bottom panels (for the DL17 model) of Figure 9. The corresponding power spectra for the velocity maps in Figure 9 are shown by the thick solid lines in Figure 10. We only show the full power spectra at 45 kpc, as Sgr does not substantially contribute to the overall power spectrum at larger radii (with the exception of the DL17 model at 70 kpc; this results from stars accelerating toward and away from the stream apocenters, which are at larger distances for the DL17 model than for the Erkal model). From left to right, power spectra for v_R, v_θ, v_ϕ are plotted; the thick solid lines are the power spectra from the combination of the GC19 simulation with the Sgr models (top panels show the results when using the Erkal model; bottom panels show the results from the DL17 model) at 45 kpc. The dashed line shows the power spectrum from the GC19 simulation (excluding Sgr). Dotted–dashed lines show the difference between the halo including Sgr and excluding Sgr (i.e., the contribution of Sgr to the overall power spectrum), at 45 kpc (purple), 70 kpc (orange), and 100 kpc (blue), computed in 5 kpc shells.

Zoom In Zoom Out Reset image size

The dotted–dashed lines in Figure 10, representing the contribution to the power spectrum due to Sgr, have a similar morphology to the power spectra discussed in Section 4 and in the Appendix: they are characterized by a sawtooth pattern, with peaks at odd ℓ values in v_R and v_θ. Figure 10 shows that the two Sgr models result in different signatures. For the Erkal model, including Sgr increases the peak at ℓ = 1 in v_R, while the v_ϕ power spectrum is mostly unaffected. The DL17 model hardly affects the low ℓ power in v_R, while the power in v_ϕ is slightly enhanced. Both models result in higher power in v_θ at all ℓ; at 45 kpc, v_θ is the only component of motion for which the signal from Sgr is comparable to the signal of the LMC-induced DM wake.

While including Sgr does increase the overall power in orders ℓ that contain signatures of the LMC-induced DM wake, the features sensitive to the Collective response (the ℓ = 1 peak in v_R, as well as the monopole ℓ = 0 peak in v_θ) are stronger at all distances than the power due to both Sgr models alone (the dotted–dashed lines in Figure 10; though we note that the DL17 model does substantially increase the power of the monopole in v_θ). In addition, the signatures from the Collective 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 the DL17 model in v_θ). While the signal from Sgr in these key ℓ values does not overwhelm the signal from the LMC, it does contribute to the overall power in the orders sensitive to the wake signatures at 45 kpc. Sgr also contributes power in specific l, m modes that are sensitive to the LMC-induced DM wake (e.g., l = 1, m = 0 in v_R); however, the phases of the different coefficients are not strongly affected by Sgr in the low ℓ values. Modeling the influence of the inclusion of Sgr will be essential in quantifying the strengths of the LMC-induced wake components observationally at smaller Galactocentric distances.

In summary, while Sgr stars will affect the power spectra of the MW halo velocity field at smaller Galactocentric distances, the signatures from the LMC-induced DM wake as predicted by GC19 should still be distinguishable from the Sgr signatures. The primary signature of the Transient response (power in ℓ = 2 in v_ϕ) is largely unaffected by the inclusion of Sgr for both models. In v_R and v_ϕ, the dominant features in the power spectra that arise owing to the LMC-induced DM wake are stronger than the features arising owing to the inclusion of Sgr, and the power spectra have different morphologies. Because Sgr stars are likely to contribute power at ℓ = 1 in v_R and ℓ = 0 in v_θ, modeling Sgr will be important in characterizing the Collective response at closer distances in the halo (∼<50 kpc).

In addition, Sgr has been extensively studied, and our knowledge of its velocity structure is better known than ever before with the release of Gaia DR2 (e.g., Antoja et al. 2020; Ibata et al. 2020; Ramos et al. 2020). Like the substructure studied in BJ05, the debris from Sgr is also well known to be overdense on the sky, and many halo studies remove stars believed to be associated with the stream. Therefore, given that we have shown that the wake signatures should still be identifiable even if all Sgr stars are included in the analysis, the prospects for observing these signatures only improve if debris from Sgr is subtracted, even if this subtraction is imperfect.

However, we note that we have only considered the effects of including Sgr stars in the analysis of the stellar halo, and not the effects the infall of the Sgr dwarf may have had on the MW disk and halo over the past ∼6–8 Gyr. Depending on the assumed mass of Sgr, it may have resulted in a substantial shift in the MW barycenter and also caused a wake in the MW DM halo (Laporte et al. 2018b). GC19 compare the magnitude of the density perturbations that arise owing to Sgr and the LMC using the Laporte et al. (2018b) simulations and find that the contribution from Sgr is negligible (see GC19's Figure 26); however, they did not discuss the perturbations to the velocity field. We leave the exploration of these effects to future work.

In this paper, we use spherical harmonic expansion to describe the perturbed velocity fields of the MW as a result of the LMC's infall using the simulations from GC19. We explore the ways in which Galactic substructure might obscure the signatures from the wake in the power spectrum, using the BJ05 simulations and two models for the Sgr stream. We summarize our primary findings as follows:

• 1.  
  We study the perturbation to the velocity field caused by the LMC-induced DM wake using the simulations from GC19. We find that low-order spherical harmonic expansion of the velocity field in these simulations usefully captures the salient features of the LMC-induced DM wake. We found that increasing power with Galactocentric radius in ℓ = 1 in v_R and ℓ = 0 in v_θ are signatures of the Collective response. At 45 kpc (near the LMC), powers in ℓ = 2 in v_R and v_ϕ are signatures of the Transient response. We find that the amplitudes of the power spectra scale with LMC mass.
• 2.  
  We investigate how Galactic substructure might affect the angular power spectrum of the MW's velocity field, using the BJ05 simulations with artificially constructed accretion histories. We find that massive, recent accretion causes large-scale, high-amplitude fluctuations in the velocity field. Velocity substructure arising as a result of debris from recent, massive satellites creates much more power in the power spectrum of the velocity field than the perturbation due to the LMC. However, the MW is not believed to have experienced much recent, massive accretion, with the exceptions of Sgr and the LMC itself.
• 3.  
  Given that Sgr is the most recent, massive accretion event experienced by the MW (with the exception of the LMC), we investigate how Sgr stars could impact measurements of the overall MW power spectra and our ability to measure the signatures associated with the LMC-induced DM wake. The power spectrum on large scales (i.e., low ℓ) remains generally dominated by the signatures from the LMC-induced DM wake; this result complements the GC19 findings that the amplitude of the density wake induced by the LMC's infall is much greater than the density wake induced by Sgr. In addition, overall power due to Sgr decreases as a function of distance, in contrast to the Collective response signatures. However, including Sgr stars does increase the power in modes that are sensitive to the Collective response, especially at 45 kpc. Care should therefore be taken to model the impact of Sgr stars on the power spectrum in studies attempting to use this method for detecting and quantifying the Collective response.

Based on our findings, performing spherical harmonic expansion in the MW velocity field could be a method for identifying and characterizing the LMC-induced DM wake, which would in turn provide constraints on the mass and orbital history of the LMC. There remain technical challenges associated with implementing this technique with observational data (e.g., incorporating measurement uncertainties, limited sky coverage of spectroscopic programs, combining data from different surveys); we leave a detailed exploration of how to estimate the spherical harmonic expansion coefficients from realistic data to future work. However, the future is bright given the upcoming observational programs that will map our Galaxy's phase-space structure. The first two Gaia data releases have already transformed our understanding of the MW's kinematic structure; with future releases from the Gaia mission in conjunction with the Gaia spectroscopic follow-up programs (e.g., DESI, 4MOST, WEAVE), as well as future astrometric (e.g., Rubin Observatory Legacy Survey of Space and Time, WFIRST) and spectroscopic programs (e.g., SDSS-V MW Mapper), the halo velocity field will be better known than ever before.

While GC19 only include the MW and the LMC, they reveal the complex behavior that arises in the phase-space structure of the MW halo owing to the infall of the LMC. However, many complications remain to be addressed. While we explored how Galactic substructure might obscure the signals from the wake, the assumption of smoothness of the MW halo in GC19, as well as their mapping of the stellar halo onto the DM halo, remain important to keep in mind. While much of the stellar halo is phase mixed in the inner regions of the MW, at larger distances (e.g., r > 50 kpc) we may need to rely on debris that is not phase mixed in order to see wake signatures. In addition, the stellar halo may not be in equilibrium with the DM halo for several reasons. First of all, simulations show that accretion is not the only mechanism by which stellar halos form: simulated halos show substantial fractions of stars that formed in situ (e.g., Zolotov et al. 2009, 2010). Yu et al. (2020) show that in situ stars (formed in outflows of the host galaxy) can be ejected to large distances in the halo and can compose 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 release to argue that the MW halo has an in situ component. Second, because 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 dwarf galaxies have been observed to vary over orders of magnitude (e.g., McConnachie 2012): the relative amounts of mass accreted in stars and DM are therefore not constant. Finally, a significant fraction of the DM accretion is "smooth" (e.g., Angulo & White 2010; Genel et al. 2010), in contrast to the relatively "lumpy" accretion that builds up the stellar halo. While methods have been developed using cosmological simulations to determine DM halo distributions based on stellar halo distributions (Necib et al. 2019), these complex effects were not included in the creation of stellar halos from DM halos in GC19. In addition, while GC19 varies the mass of the LMC, they assume a single mass for the MW and do not simulate a range of orbital histories.

Many of these complications can be addressed by studying DM wakes more generally in cosmological simulations. These simulations contain both in situ and accreted halo components and have experienced a variety of accretion histories. By applying this technique to encounters between MW-like galaxies and massive satellites in cosmological simulations, we can identify wake signatures for a range of mass ratios and orbital histories. In the present era of wide-field 3D kinematic data sets coupled with detailed high-resolution simulations, we can move beyond equilibrium models and develop new methods for characterizing the complex processes that have shaped our Galaxy's formation.

We thank the anonymous referee for the helpful comments and suggestions. E.C.C. and R.L. are supported by Flatiron Research Fellowships at the Flatiron Institute. The Flatiron Institute is supported by the Simons Foundation. A.D. is supported by a Royal Society University Research Fellowship. K.V.J.'s contributions were supported by NSF grant AST-1715582. N.G.-C. and G.B. acknowledge support from HST grant AR 15004 and NASA ATP grant 17-ATP17-0006. C.F.P.L.'s contributions are supported by world premier international research center initiative (WPI), MEXT, Japan. This project was developed in part at the 2019 Santa Barbara Gaia Sprint, hosted by the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara. This research was supported in part at KITP by the Heising-Simons Foundation and the National Science Foundation under grant No. NSF PHY-1748958. E.C.C. would like to thank Andrew Wetzel, David Hogg, Adrian Price-Whelan, and David Spergel for helpful scientific discussions.

Software: Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), Healpy (Zonca et al. 2019), IPython (Perez & Granger 2007), 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).

In Section 4, we see that the power spectra for all of the BJ05 halos with debris that is not yet phase mixed are characterized by a sawtooth pattern, with peaks at odd ℓ values for v_R and v_θ and peaks at even ℓ values for v_ϕ. To understand why these features in 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 A1. The top panels of Figure A1 show intensity maps, while the bottom panels show the corresponding power spectra. While a single spot gives rise to a rather smooth power spectrum, with power in even as well as odd ℓ values, two spots of opposing sign (located 180° apart from one another) result in a power spectrum with peaks at odd ℓ values. As seen in the right panels, banding can also give rise to the sawtooth pattern: banding at constant intensity yields peaks at even ℓ values, whereas a band weighted by a dipole creates peaks at odd ℓ values.

Zoom In Zoom Out Reset image size

To understand why the halos have these features in their velocity maps and power spectra, we explore the properties of two specific accretion events. First, we consider BJ05 satellite 1066 (${M}_{\mathrm{Sat}}=9.1\times {10}^{10}\,{M}_{\odot },{t}_{\mathrm{Sat}}=4.58\,\mathrm{Gyr},{J}_{\mathrm{Sat}}/{J}_{\mathrm{circ}}=0.68$). This halo is accreted by both the recently accreted BJ05 halo and the high-L_sat halo. The top panels of Figure A2 show the positions of all stars from this satellite in Cartesian coordinates, color-coded by their radial velocities. Dashed lines indicate the radial range of 30–50 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 located between 30 and 50 kpc from the center, we see that our cross section does not include the apocenters of these shells (where stars have v_R ∼ 0 $\mathrm{km}\,{{\rm{s}}}^{-1}$), but rather stars that are moving quickly either toward an apocenter (yellower spots) or away from apocenter back toward the center of the galaxy (bluer points). In the radial velocity map, the patches of stars with opposite velocity are separated by approximately 180° from each other in this projection. This accretion event has a much narrower range of velocities in v_θ and v_ϕ; the v_θ map has a fairly smooth continuum, going from positive v_θ as the stars moves toward and away from the first apocenter to negative v_θ as the stars pass through the second apocenter. This results in banding weighted by a dipole, a version of the map shown in the rightmost panel of Figure A1. In v_ϕ, the velocities are always negative; the velocity map appears as a relatively constant band.

Zoom In Zoom Out Reset image size

The power spectra for these velocity maps are shown as the purple curves in Figure A3. This radial accretion event has the most power in the radial velocity power spectrum; the power spectra peak at odd ℓ values for v_R and v_θ, while the v_ϕ power spectrum is dominated by peaks at even ℓ values.

Zoom In Zoom Out Reset image size

As a different example, we can consider the debris from Satellite 1288, shown in the bottom panels of Figure A2. This relatively low mass (${M}_{\mathrm{Sat}}=5.54\times {10}^{10}\,{M}_{\odot }$), recently accreted (${t}_{\mathrm{Sat}}=6.55\,\mathrm{Gyr})$, circular (${J}_{\mathrm{Sat}}/{J}_{\mathrm{circ}}=0.99$) accretion event leaves a w-wrapped kinematically cold stream that lies approximately in the (x, y) plane. When we take the cross section of stars from 30 to 50 kpc, we capture many of the stream stars, which form nearly continuous bands in the Mollweide projection. Because the stream is not perfectly circular, at opposite points in the orbit, the radial and polar velocities will be of similar magnitude but of opposite sign, while the tangential velocity is approximately constant along the stream.

The power spectra are shown as the green curves in Figure A3. 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 ℓ values, while the v_ϕ power spectrum is peaked at odd ℓ values.

The fact that power spectra have these sawtooth patterns is indicative that spherical harmonic expansion of the different velocity components 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 address this; we leave this to future work.