Quantifying the impact of the Large Magellanic Cloud on the structure of the Milky Way’s dark matter halo using Basis Function Expansions

The suite of N–body simulations used in this work were presented in G19

The suite consists of eight high-resolution simulations, with dark matter particle mass $m_{p}=1.5\times 10^{4}M_{\odot}$¹¹1Note, that this value was incorrectly quoted as $4\times 10^{4}$ M_⊙ in Table 1 of G19, run with Gadget-3 (Springel et al., 2008). We use two MW models, each with the same Hernquist DM halo density profile with virial mass $M_{vir}=1.2\times 10^{12}$M_⊙, an exponential disk of mass $M_{d}=5.78\times 10^{10}$M_⊙, and a bulge of mass $M_{b}=0.7\times 10^{10}$M_⊙. We adopt two different with different halo distribution functions, where one is isotropic and the other is radially biased. These are represented by anisotropy parameters, $\beta(r)=0$ and $\beta(r)=-0.5-0.2\frac{d\rm{ln}\rho(r)}{d\rm{ln}r}$, respectively.

For the LMC, we have four models with different halo masses but with the same enclosed mass within 9 kpc, fixed by the rotation curve of the LMC (van der Marel & Kallivayalil, 2014). This observational constraint implies that the concentration of the halo is set for each LMC mass model and that the inner mass profile is similar for each LMC mass model. This has important consequences to the halo response, as we will illustrate.

Both the MW and the LMC halos were initialized with GalIC (Yurin & Springel, 2014) with a Hernquist DM halo, but matched to an NFW halo in the inner parts following van der Marel & Kallivayalil (2014). Table 1 summarizes the main properties of these simulations. We take simulation #7 as our fiducial simulation since the MW halo model (radially biased) and LMC mass are consistent with current accepted values (e.g., Moster et al., 2013; Cunningham et al., 2019).

The DM halos of the MW and the LMC were evolved separately in isolation for 2 Gyrs in to order to guarantee energetic equilibrium in the halos. The LMC was placed at the virial radius of the MW halo (ICs values are reported in Table 8 of G19). The simulations were run for $\sim$2 Gyr, and follow the evolution of the LMC on its first infall into the MW. At the present time the position and velocity vectors of the simulated LMC are within 2$\sigma$ of the observed values reported in Kallivayalil et al. (2013).

To analyze the present-day snapshot from our N–body simulations we use BFEs. Several BFE expansions have been developed in the last three decades. For DM halos, Clutton-Brock (1973) built a BFE whose zeroth-order basis is the Plummer profile (Plummer, 1911). Similarly, Hernquist & Ostriker (1992) built a BFE based on the Hernquist profile (Hernquist, 1990); given its accuracy, this expansion has been widely used in the literature (e.g Johnston et al., 2001, 2002b, 2002a; Sanders et al., 2020). Recently, Lilley et al. (2018b) presented an analytic flexible BFE expansion based on the NFW profile (Navarro et al., 1997), and presented a more general

family of double-power law expansions parametrized with two parameters in (Lilley et al., 2018a). Another option is to solve the Poisson equation numerically. which is a particular form of the Sturm-Liouville equation (Weinberg, 1999). For example, in Dai et al. (2018), the NFW basis was solved numerically for the Eris cosmologically simulated halo (Guedes et al., 2011).

In our case, the Hernquist BFE is a natural choice since our simulations were initialized with Hernquist DM halos. Furthermore, (Sanders et al., 2020) demonstrate a higher performance and accuracy of the Hernquist BFE with respect to the other expansions.

In this section we summarize the main equations of the Hernquist expansion. For completeness, the derivation can be found in Appendix A. For a detailed and comprehensive derivation of this expansion we refer the reader to sections 2.2 and 3.1 in Hernquist & Ostriker (1992) or section 2.1 in Lowing et al. (2011).

The density and potential for a system described by particles is represented by the following expansions:

Where the coefficients $S_{nlm}$ and $T_{nlm}$ are defined as:

The expressions for $\rho_{nl}$ and $\Phi_{nl}$ can be found in Appendix A. The normalization factor $I_{nl}$ is defined as:

And $K_{nl}$ is defined as:

Python public implementations of the Hernquist & Ostriker (1992) BFE can be found in the galactic dynamics Python-libraries Gala (Price-Whelan, 2017; Price-Whelan et al., 2017) ²²2http://gala.adrian.pw/en/latest/examples/Arbitrary-density-SCF.html, and Galpy (Bovy, 2015). An example of how this BFE decomposes the density and potential for a prolate, oblate, and triaxial halo can be found in Appendix D. In this paper we used a modified version of the Gala library, which implements a new noise reduction method developed in this work (Appendix B) and is also adjusted for reading N–body snapshots.

A challenge with BFEs is choosing the order of the expansion, i.e. $n_{max}$ and $l_{max}$. This choice is not obvious and depends on the number of particles and morphology of the system. Weinberg (1996) developed a method to truncate the expansion by studying the noise generated by the discrete nature of the simulation. In Appendix B, we discuss the nature of the noise in BFE and how the noise is reduced by truncating the expansion.

Here we utilize BFEs for MW particles alone to describe the MW’s DM halo response due to the passage of the LMC. Panel $c)$ in Figure 1 shows the BFE for the present-day MW’s DM halo without the LMC particles. Panel $f)$ shows the density contrast for the MW’s DM halo, computed as:

where $\rho_{000}$ corresponds to the density computed with only the first term in the expansion.

Panel $f)$ clearly illustrates the halo response.

A collisionless self-gravitating system responds to a perturbation through excitation of specific modes (e.g: Kalnajs, 1977; Weinberg, 1993). For the collisionless Boltzmann equation, the spectrum of solutions is represented by discrete and continuous modes. Overall, the discrete modes are sustained by resonances (e.g: Weinberg, 1994) and damp slowly over time (the Collective Response) while the continuous modes (Vandervoort, 2003) are non-resonant and damp quickly (the dynamical friction DM Wake). In the case of a satellite sinking in a larger host, the host response is commonly known as the ”wake” which encompasses both discrete and continuous modes. However, given their markedly distinct and striking morphologies we make a distinction between the two as follows:

LMC’s Dynamical Friction Wake: The wake is the overdensity of DM particles trailing the LMC as illustrated in Figure 2. The LMC excites modes that will phase-mix on time scales shorter than an orbital time. These modes can be thought of as wave-packets that are excited by the perturbation. Classically, the dynamical friction wake is described as the local scattering of DM halo particles as a massive perturber travels through the medium of a larger host. In the case of a homogeneous and infinite background medium, this leads to the classic dynamical friction equation derived in (Chandrasekhar, 1943). Mulder (1983) and Weinberg (1986) showed that dynamical friction was the result of the wake formation. Here we define the LMC’s dynamical friction wake as the non-resonant continuous modal response of the halo to the passage of the LMC. As such, the dynamical friction wake will damp quickly after a few crossing times.

Collective Response:

The Collective Response is shown in the upper left of Figure 2. It is a large overdensity in the halo, mainly in the Northern hemisphere. It has been recognized that satellite-halo interactions lead to a barycenter excursion (Weinberg, 1989; Choi, 2007). The amplitude of the barycenter excursion (or displacement) depends on the mass and orbit of the satellite. For a galaxy, the barycenter displacement is complex because the system is not a point mass and therefore the response is described by the excitation of its $l=1$ weakly damped mode (Weinberg, 1989). In contrast to the dynamical friction wake, the Collective Response is the result of the excitation of discrete modes (e.g. $l=1,m=1$), which are sustained by resonances and self-gravity. These perturbations will damp slowly over time-scales longer than the orbital time. Here we define the Collective Response as the discrete modal response of the halo to the passage of the LMC predominantly seen in the ($l=1$ mode). As such, the Collective Response will persist over many crossing times.

In the remainder of this section, we will study the combined effect of both the dynamical friction wake and the Collective Response using the MW BFE (no LMC particles).

A 3D animated rendering of the density contrast (equation 6), illustrating the halo response to the LMC’s passage, can be found here https://vimeo.com/546207117. The present-day frame from the animation for the density contrast (equation 6) is shown in Figure 2.

We do not include an expansion for the baryonic components of the MW (disk and bulge), although they are included in the N-body simulation as live components (see Section 2.1). Consequently, the MW’s DM halo and the LMC do feel the gravitational potential of the MW disk and bulge.

We use BFEs to compute the gravitational potential $\Phi$ of the LMC using all the LMC’s DM particles. The expansion is centered on the LMC, where the cusp of the LMC was computed using the shrinking sphere algorithm described in Power et al. (2003). We compute the expansion up to order $n_{max}=l_{max}=20$. The gravitational potential is computed using only the coefficients with $\Gamma>$5. We keep the halo scale length fixed during the iterative procedure since the inner region of the LMC’s halo does not change significantly during its orbit. We test this by fitting the scale length in each iteration to the LMC’s DM particle distribution.

We define particles as bound to the LMC if their kinetic energy (KE) is less than the gravitational potential energy ($\Phi$). This is similar to the iterative procedure used in Johnston et al. (1999). In practice we use the energy per particle mass and hence a particle is bound if: $\dfrac{1}{2}v_{k}^{2}<\Phi_{BFE}$, where $\Phi_{BFE}$ is computed using equation A2 at the location of each particle, $k$. The velocity $v_{k}$ is computed with respect to the LMC’s cusp.

We iterate the computation of the bound particles for every particle until we reach 1% convergence in the number of bound particles. The mean number of iterations required is $\sim$ 5.

Panel $d)$ in Figure 1 shows the BFE reconstruction of the density field of the bound LMC particles at the present time for the fiducial model. The LMC DM halo exhibits an S-shape produced by the tidal field of the MW. This is most pronounced in the galactocentric y-z plane (coincident with the orbital plane of the LMC).

The total bound mass of the LMC at the present-day ($M_{\rm{LMC}}$) depends on both the orbit of the LMC and on the DM halo profile of the LMC. Figure 3 illustrates the bound mass of the LMC at the present day, $M_{\rm{LMC}}$, as a function of the LMC’s initial virial infall mass ($M_{\rm{LMC},vir}$), for the four LMC models used in this study (all Hernquist halos).

Bound $M_{\rm{LMC}}$ ranges from $\sim 4\times 10^{10}$ (50% of the infall mass of the lowest mass LMC model) up to $\sim 7\times 10^{10}M_{\odot}$ (30% of the infall mass of the highest mass LMC model). The most massive LMC models lose a larger fraction of mass compared to the less massive models. This is expected because the concentration

of these halo models are different. Each LMC’s DM halo is initially represented by a Hernquist profile with a scale radius chosen to reproduce the observed rotation curve of the LMC (see G19). As such, low mass models have a higher halo concentration compared to high mass models. For these low concentration models, the DM particles in the outskirts of the LMC’s DM halo are less bound to the LMC than in the higher concentration models (low mass models). As a result, the ratio of bound to infall virial mass is lower in the higher mass models (see lower panel of Figure 3).

We illustrate the extent of the bound LMC’s DM halo in Figure 4 for all of our four models, in the Galactocentric $y-z$ plane. As a reference, the MW’s disk particles are also plotted. The present-day bound LMC’s DM halo is very extended, where the edge of the halo scale length (white circle) can be as close as $\sim$30 kpc from the Galactic center.

Over the last 2 Gyrs, the LMC has lost more than 50% of its pre-infall mass (see lower panel in Figure 3). This DM debris extends over large distances throughout the DM halo. The LMC’s DM debris can impact the shape of the MW’s halo (see Section 4.4) and can extend the reach of direct detection experiments to lower WIMP mass (Besla et al., 2019).

The mass of the LMC’s DM debris is found using the same iterative method described in section 3.2, but now requiring that the kinetic energy exceeds the potential energy, thus defining unbound LMC particles.

We utilize one high-order ($n_{max}=l_{max}=30$) expansion to capture all LMC particles, i.e. both bound and unbound particles. The BFE for the bound LMC only requires 20 terms. However, including the highly distorted debris, 1831 terms and a signal to noise threshold of 2 (see Appendix B) were needed to reproduce the density field computed from the LMC debris particles (top vs. bottom panels in Figure 5).

Through comparison with the BFE for the bound particles (previous section), we can identify the properties of the LMC’s DM debris. The debris is found to be a much more extended asymmetric structure, as illustrated in panel $e)$ in Figure 1 for the fiducial model (relative to panel $d$).

In Figure 5, we illustrate the density contrast of the LMC debris with respect to the MW (panel $c$ of Figure 1). This is computed using equation 6, where the denominator is now the MW density (not just the monopole). We isolate the structure of the LMC debris by masking high density regions in the LMC particle distribution (density contrast $>$1.5), which correspond to the bound LMC. Two components comprise the LMC’s DM debris: trailing and leading,

1) Trailing component: LMC DM particles that extend from the LMC out to $-200$ kpc in $\hat{z}$ and beyond 200 kpc in $\hat{y}$ (Galactocentric). The DM density enhancement of the trailing debris relative to the MW (which includes the halo response, i.e. the wake) is $\sim$100% for the fiducial model. Meaning that in the regions of the halo traced by the debris, the mass in the debris equals the local DM mass of the MW in that same volume.

This high ratio is aided by the fact that the dynamical friction wake causes a significant decrease in the MW’s DM content in that region, by $\sim 30\%-40\%$ (panel $f$ in Figure 24).

2) Leading component: LMC DM particles in the north galactic hemisphere, extending from the LMC out to 150 kpc in $\hat{z}$ and $\hat{y}$. The leading component of the LMC debris overlaps with the Collective Response (panel $d$ vs. $e$ in Figure 24). The shape of the leading component is asymmetric, resembling the shell structure typically formed through mergers of satellites on radial orbits (Piran & Villumsen, 1987; Heisler & White, 1990). The dynamics of tidal debris from satellite galaxies has been extensively studied (e.g Choi et al., 2009; Sanderson & Bertschinger, 2010; Amorisco, 2015; Hendel & Johnston, 2015; Drakos et al., 2020). In general the resulting morphology of tidal debris is governed by the satellite’s orbit, internal rotation, mass, and internal structure. However, for massive satellites, the self gravity of the debris, the satellite’s gravitational influence on the debris, and the host halo’s response to the satellite, are also important. All of these processes are captured in our simulations.

In order to fully constrain and predict the morphology of the LMC’s debris, simulations that cover a more complete parameter space are needed. For example, varying the density profile of the LMC and including perturbations and mass loss from the SMC. We conclude that, in our simulations, the mass of the LMC’s DM debris is non-negligible and can dominate the density field of the DM around the MW in some regions of the southern hemisphere.

The distribution of the debris is highly asymmetric and evolves in time - it cannot be well-represented by the infall of a static LMC halo. We expect that the distribution of the LMC DM debris can affect the dynamics of stellar streams, GCs, and satellite galaxies; this will be studied in subsequent papers.

With the three optimal BFEs computed for each component of the system (MW, LMC, LMC+Debris) we can now quantify the density, potential, acceleration and shape of the combined MW–LMC system in a compact and accurate manner. Specifically, we combine the BFE for the MW (section 3.1) with the BFE for the LMC+Debris (section 3.3). The combined expansion is seen in panel $a$ in Figure 24).

For the fiducial simulation, the combined expansion consists of 236 coefficients for the MW and 1831 for the LMC+debris (see Table 2 and previous section). Each combined expansion contains the information of the entire 100 million particle simulation.

Figure 6,

illustrates the projected density (left panel), acceleration (middle panel), and potential (right panel) fields for the combined MW–LMC system computed with the BFE. Each row shows a different size-scale of the system, from a radius of 75 kpc (top panel) to 150 kpc (bottom panel). Note that, to create these combined potentials, the center of each of the expansions is different. The expansion of the MW

is centered on the MW COM, while the expansion of the LMC+Debris is centered on the LMC’s COM.

The inner region ($<$35 kpc) of the MW’s DM halo is governed by the presence of the disk and the bulge. The contribution from the MW disk to plotted the density, potential and acceleration fields, was added by fitting the simulated MW disk from the N-body simulations with a Miyamoto-Nagai profile (Miyamoto & Nagai, 1975). The bulge was also added analytically following the Hernquist profile used in the initial conditions from G19. While it is possible to use a BFE for the bulge and the disk, respectively (Petersen et al., 2019), in this study we are primarily concerned with the perturbations at large distances, where the deviations from the analytic representations of the disk and bulge are negligible.

An interesting feature observed in the simulated acceleration field is the local minimum value between the MW disk and the LMC. This minimum is located at $\vec{r}=(-2.5\hat{x},-31\hat{y},-20\hat{z})$ kpc (see white cross in bottom right panel of Figure 6), which is along the separation vector between the LMC and MW. The position corresponds to ($l$=280, $b=-$32) in galactic coordinates. The location of this minimum does not vary among our 8 simulations since the enclosed mass within 50 kpc of the LMC is similar in all models. The combined acceleration field provides a snapshot of the non-axisymmetric accelerations experienced by objects in orbit about the MW at the present-day.

Outside 35 kpc, the perturbations from the LMC become apparent. Formed by the combined contribution of the MW’s DM halo response, the bound LMC and the LMC’s DM debris, the perturbation is anisotropic. The DM distribution outside 35 kpc is elongated in the $\hat{y}$ direction, which is in the direction of the DM dynamical friction wake and the DM debris. Hence, the DM halo shape/distribution contains information about the location and amplitude of the DM dynamical friction wake.

Ultimately, our goal is to quantify the perturbation caused by the LMC in the MW’s DM halo using the BFEs, we will discuss this in Section 4. We will discuss how these results change as a function of the MW anisotropy profile and LMC mass in section 5.2.

One of the advantages of using BFEs is the ability to decompose the halo response into modes whose amplitudes are the coefficients of the expansion. The contribution to the total gravitational potential energy $U_{nlm}$ of each mode can be computed with the corresponding coefficients. We quantify the halo response to the passage of the LMC in terms of the energy in each expansion term. The total gravitational potential $U$ for the basis expansion is defined:

Note that equation 8 has units of energy, the coefficients $S$ and $T$ are unit-less, but the units comes from the normalization factor $I_{nl}$ (see Eq. A7). The coefficients used in the energy calculation are already smoothed, as described in Section B. The most energetic coefficients represent the modes that contribute the most to the total gravitational potential energy of the system.

Figure 7 shows the gravitational potential energy, $U_{nlm}$, of all coefficients decomposed in $n$ and $l$, with $m$ indicated per column. These coefficients were computed for our fiducial simulation #7. The main characteristics of the MW BFEs are summarized as follows:

1. 1.

   Globally, the most energetic term is the monopole, which contributes $99$% of the total gravitational potential energy in both cases. This term is the spherically symmetric Hernquist DM halo.

2. 2.

   The radial modes, $n>0$, contribute up to order $\sim$ 20, for the $m=l=0$ case.

3. 3.

   The angular, $l$, modes appear up to order $l=14$, but most of the energy in those terms is contained within order $l<6$. The angular, $m$, modes exhibit contributions up to $m=8$ and $m=12$; however, most of the energy is contained in the $m=0,1,2$ modes (shown in Figure 7).

Although globally the halo energetics are reasonably described by the original, unperturbed Hernquist DM halo potential, the existence of energetic angular $l,m$ modes prove the existence of asymmetries induced by the LMC’s passage within the halo. Furthermore, both the halo Collective Response and the dynamical friction wake are localized perturbations that will have a non-negligible impact on the density (G19) and shape of the halo (see Section 4.4).

Studying the effect of the LMC’s DM debris is also crucial to interpret the shape of the MW’s halo. Since we have separate expansions for the LMC+DM debris and the MW, we cannot self-consistently quantify the increase in energy in the MW BFE from the LMC’s debris. However, we will study the effects of the debris on the density field of the MW-LMC system in section 4.4.

In this section we qualitatively describe the terms that most contribute to the halo response (no LMC BFE). Figure 8 shows the projected density field of the halo response (Collective Response and dynamical friction wake), dissected using different sets of coefficients, ranked by their order. Each column shows the reconstruction up to a different maximal order in $l$ and $m$, from $l_{max}=m_{max}=0$ (left column) up to $l_{max}=m_{max}=4$ (right column). In each row the maximal order of the radial $n$ term is increased. This means that the upper left panel only contains two coefficients, $n=0$ and $n=1$, while the bottom right panel has the most terms, up to $n_{max}=20$ and $l_{max}=4$. Below we list the the salient features for each columns.

The left column illustrates spherically symmetric density fields. As the order in $n$ increases, each term contributes to a perturbation at a characteristic radius. The combination of these $n$-modes creates a more complex perturbation that sets the radial scale of the halo response. For example, in the lower left panel, the

overdensities appear mainly beyond 50 kpc, which coincides with where the dynamical friction wake starts to pick up.

The middle and right columns: As we increase the order in $l$ and $m$, we see deviations from spherically symmetric perturbations. For $l_{max}=m_{max}=1$ we clearly see a dipole in the density field, corresponding to the Collective Response. The amplitude of the $l_{max}=2$ mode is responsible for applying torques to the MW’s disk (e.g Weinberg, 1998b; Gómez et al., 2013, 2017; Laporte et al., 2018b, a). The bottom row illustrates that the radial structure of the dipole is defined by the $n$ terms. The level of asymmetry needed to reveal the dynamical friction wake is reached at $l_{max}=m_{max}$=4.

Distinguishing between the Collective Response and the dynamical friction wake is not straightforward, but this exercise illustrates the order of the modes that contribute to these structures. Dissecting all the dynamical features in the simulations would require additional analysis of the coefficients and its time evolution.

In our simulations we start with an idealized spherical halo. However, cosmologically, halos are not expected to be spherical (e.g White & Silk, 1979; Zemp et al., 2011; Vera-Ciro et al., 2011; Chua et al., 2019; Emami et al., 2020).

In this section, we use the BFE to analyze the shape of the MW–LMC halo shape and compare it with MW–like halos proposed in the literature. More details about the BFE of triaxial halos can be found in Appendix D.

We compute BFEs for simulations of oblate, prolate, triaxial halos whose main properties are summarized in table 3. The triaxial halo model is the mean result for MW–like halos found in the Illustris cosmological simulation (Chua et al., 2019). In addition, we include the shape measurement derived by Law & Majewski (2010) from fitting the Sgr stream These halos are characterized by the three axes describing the shape of ellipsoids: $a$ is the major, $b$ the intermediate, and $c$ the minor axis. We then utilize axis ratios to define the ellipsoid: $s_{\rho}=c/a$, the minor to major axis ratio; and $q_{\rho}=b/a$, the intermediate to major axis ratio. The triaxiality of the halo can be defined in terms of $s$ and $q$ as:

If $T_{\rm{qs}}\geq 0.67$, the halo is considered to be prolate. If $0.67\geq T_{\rm{qs}}\geq 0.33$ the halo is triaxial, and if $T_{\rm{qs}}\leq 0.33$ the halo is oblate.

Figure 9, summarizes our findings for the BFE for our simulated oblate, prolate and triaxial halos. Each square in the $l-m$ grid shows the energy corresponding to the sum of all the $n$-modes with the same $l$ and $m$. That is, $U_{lm}=\sum_{n}U_{nlm}$.

Prolate halos only have $m=0$ modes. Corresponding to the ‘zonal’ spherical harmonics that do not depend on longitude and whose lobes are perpendicular to the plane of the MW’s disk. Oblate halos on the other hand, have a major contribution from the $|m|=l$ modes. These are represented by the “sectoral” spherical harmonics, whose lobes are parallel to the MW’s disk. The BFE for the triaxial halo is a mixture of that of the prolate and oblate halos.

The sign of the coefficients can also be used as an indicator of the halo shape. Oblate halos have negative values for the $l=2^{i}_{odd}$ modes. See Appendix  D for further details. Note that, in all cases, as the axis ratios increase, higher-order terms will be needed to properly characterize the structure. The main characteristic of these idealized halos is that the odd $l$ and $m$ modes do not appear in the expansion, as they are not radially symmetric. As such, the existence of odd modes would signify divergence from these standard halos.

With this intuition in mind, we can now move forward and analyze how a more complicated, asymmetric structure, such as the DM dynamical friction wake produced by the LMC, will manifest in the BFE space.

In Figure 10, we plot the energy contribution of each $l,m$ mode, summed over all $n$-modes, for the MW BFE (section 3.1).

The MW halo responds to the presence of the LMC, but the LMC BFE is not included in this analysis.

Comparing Figure 9 with Figure 10, we find that the most energetic coefficients in the perturbed MW halo are distinct from those observed in idealized or cosmologically motivated halos.

In the MW BFE, the monopole term is the dominant mode, followed by the $l=1$ mode, as also found by Tamfal et al. (2020) and in the kinematics of the halo particle by Cunningham et al. (2020). This corresponds to the dipole response, which is induced by the barycenter motion (see 4.4). In other words, we find strong contributions from radially asymmetric modes (odd $l$-modes), which are completely absent in the idealized halos shown in Figure 9. The LMC’s impact on the MW halo is not expected to be distorted in a manner consistent with a prolate, oblate or triaxial DM distribution.

In Figure 11, we plot the energy as a function of $l$ and $m$ modes, respectively, to better illustrate the difference between our simulated MW and idealized halos.

Overall, the energy of the even $l$ and $m$ terms is higher for the chosen oblate, prolate, and triaxial halos than in the MW, although these idealized halos were chosen to be somewhat extreme.

There is pronounced difference between the MW BFE and the LM+10 halo, particularly between $l=2-6$. This implies that the MW halo response to the LMC alone is not sufficient to generate the deformations measured by (Law & Majewski, 2010). Also, it is still unknown if the response to the passage of the LMC of an initial triaxial MW can explain the resulting halo shape of Law & Majewski (2010). In addition, the direct torque from the LMC on the Sgr. stream (Vera-Ciro & Helmi, 2013; Gómez et al., 2015) has to be taken into account to properly interpret the measured shape of Law & Majewski (2010), as demonstrated in Vasiliev et al. (2020). Note that we have not included the LMC BFE in this analysis.

Evidence of the direct torque from the LMC in stellar streams has been detected in the Orphan stream (Erkal et al., 2019b), where a highly flattened, prolate

MW DM halo was preferred in order to reproduce the morphology of the stream using N-body simulations. However, our analysis indicates that the halo response cannot be captured by such static, axi-symmetric models. This reinforces the point that efforts to measure halo shapes using streams must to take into account both the direct torque from the LMC as well as the halo response, as recently shown in Vasiliev et al. (2020) for the Sag. Stream.

Interestingly, the average halo from Illustris measured by Chua et al. (2019) (green stars) is triaxial with power in the $l=4,6$ terms consistent with that of our simulated MW. This indicates that it may be challenging to disentangle the impact of the LMC from the structure induced by the cosmological assembly history of the MW. The quadrupole terms ($l$ and $m$) may serve as a discriminant of the halo triaxiality and effect from the LMC. A boost in the quadrupole relative to the higher order terms may signal underlying triaxiality. The slope of the energy as a function of $l$ and $m$ is correlated with the degree of triaxiality ($T$), where the slope flattens as the axis ratios decrease. Also note that the power increases for more highly-triaxial halos, as shown by the black stars.

This analysis illustrates that it is not straightforward to disentangle the halo response due to a passage of a massive satellite from the halo’s intrinsic/initial shape in simulations. A full diagnostic will require a proper assessment of the covariances of the coefficients to quantify the halo response similar to Ghil et al. (2002); Darling & Widrow (2019).

We now turn our attention towards understanding the origin of the MW’s shape in our simulations. The present-day shape of the density field of the MW’s DM halo is sensitive to the recent infall of the LMC. The following physical processes take place over the last 2 Gyrs:

1) Density asymmetries occur owing to the presence of the LMC (bound DM halo and unbound DM debris) and the dynamical friction wake;

2) Barycenter motion: Direct torques from the LMC and dynamical friction wake move the barycenter of the MW–LMC system. Since the dynamical times of the outer halo of the MW are longer than those of the inner halo, the inner halo will respond faster to the presence of the LMC. Consequently, the COM of the inner halo shifts with respect to the outer halo as a function of time. As a result, the barycenter of the MW–LMC system not only changes over time, but also as a function of Galactocentric distance. Because of the barycenter displacement, the density in the northern hemisphere is larger than that in the south, at a given Galactocentric radius ($>$ 30 kpc, measured from the disk COM). The barycenter motion is a manifestation of the Collective Response.

As a consequence of both processes, observers in the MW’s disk should observe an apparent reflex motion, as has been recently measured (Petersen & Peñarrubia, 2021; Erkal et al., 2020).

In this section, we quantify the contribution of the density asymmetries and the barycenter motion to the shape of the MW’s DM halo density distribution at the present day.

As stated earlier, the COM of the MW’s disk has moved as a function of time owing to the infall of the LMC (e.g Gómez et al., 2015; Petersen & Peñarrubia, 2020), but most of this motion has occurred recently. This is illustrated in the top panel of Figure 12, which shows the change in COM position of the MW’s disk with respect to its present day location, as a function of lookback time. The bottom panel shows the cumulative COM motion as a function of lookback time. From both of these plots we conclude that the COM of the MW’s disk moves more than 50% of the total displacement in the last 0.5 Gyrs, which corresponds to a distance of 20 kpc for the fiducial LMC mass model (LMC3), in agreement with previous studies (Gómez et al., 2015; Erkal et al., 2019b; Petersen & Peñarrubia, 2020). The rapid motion of the disk COM over such a short time-scale implies that the impact on specific objects in orbit about the MW disk will depend on their radius and orbital period.

The entire halo of the MW will not respond as a rigid body to the COM motion of the disk. Instead, the motion of the halo will vary as a function of radius. The amplitude of the corresponding barycenter motion of the MW’s halo depends on the infall mass of the LMC, its orbit, and MW’s halo response.

The inner regions of the halo will move with the disk, but the outer regions will lag, and hence are displaced from the disk COM. Our goal is to distinguish the radius of this inner halo and characterize the relative COM motion of the outer halo.

Focusing on the present day halo, we quantify the inner halo as those regions that follows the COM displacement of the disk. We compute the COM of isodensity contours of the halo using the BFE, as a function of distance. We further distinguish the effects of the LMC and the LMC’s DM debris in these calculations by comparing results using the BFE with or without those particles (see section 3).

The left panel of Figure 13 shows the resulting COM position of isodensity contours within the present day halo of the MW, as a function of distance, $r_{ell}$, which is the ellipsoid radius:

In the right panel of Figure 13, we provide a visualisation of isodensity contours in the combined MW–LMC system.

The left panel of Figure 13 illustrates that the COM of the MW halo is not coincident with that of the MW disk at all radii. Within $r_{ell}\sim$30 kpc, the MW halo moves coherently with the MW’s disk COM (pink regions in both left and right panels).

However, beyond $r_{ell}\sim$30 kpc, isodensity contours are distorted by the LMC (right panel). Consequently, the halo COM position begins to diverge from that of the inner halo. These outer radii are affected by different factors. When MW particles alone are considered (black dashed lines), the offset reflects the halo response, including the barycenter motion. The inclusion of bound LMC particles increases the displacement from $r_{ell}=$ 30-40 kpc and 50-90 kpc (orange shaded region). The inclusion of unbound debris particles from the LMC augments the COM displacement at $r_{ell}$ larger than 90 kpc (blue shaded region). The orange line denotes the COM displacement for the entire system combined. Although we have computed the COM within isodensity shells, the measured displacement should also reflect the behavior of the orbital barycenter of the MW–LMC system.

The barycenter motion manifests to observers as a reflex motion in velocities Petersen & Peñarrubia (2021); Erkal et al. (2020). In particular, it induces a dipole pattern in the radial velocities and net vertical velocity $v_{z}$ of stars/DM in the outer halo ($>30$ kpc). The northern hemisphere appears to be redshifted while the southern is blueshifted. We quantify this dipole pattern by computing a spherical harmonic expansion; the magnitude of the contribution from the $\ell=1$ modes (i.e., the dipole modes), as a function of radius, is shown in Figure 14. The power in the dipole increases steeply after 30 kpc, consistent with our findings with the COM analysis.

While the barycenter motion and the reflex motion can be estimated to first order using simple sourcing of rigid potentials as done in earlier works (e.g: Gómez et al., 2015; Erkal et al., 2019b), these models cannot capture the higher order variations as a function of distance which arise from the complex marriage between contributions from tidal stripping and response of a halo a sinking perturber, which give rise to different behaviors in the inner and outer halo of the Galaxy.

Here we outline specific cases where the barycenter motion of the MW must be accounted for. The barycenter motion manifests as both a spatial displacement and a reflex velocity in the phase space properties of stars in the outer halo ($>$ 30 kpc), relative to an observer in the disk. The barycenter motion also varies as a function of radius and LMC infall mass (see Figures 12 and 14).

Models of the Sagittarius stream: Models of the Sagittarius stream must include the gravitational influence of the LMC. Early works by Law & Majewski (2010); Vera-Ciro & Helmi (2013) showed that the LMC can have an important contribution to the shape of the potential. Gómez et al. (2015) showed that the LMC is perturbing the orbit of Sagittarius and hence its tidal debris will not lie in the plane of the orbit.

Based on our analysis we expect that stars within the Sagittarius stream will exhibit the barycenter motion of the inner halo (both in terms of the displacement and reflex motion). The Sagittarius stream is unique in that it exists in both the inner, intermediate and outer halos, and thus may show interesting phase space properties at distances larger than 30 kpc vs. at smaller radii.

Stellar streams: We expect the above to be true of any stream that crosses the inner and outer halos. Note that for any stream outside of 30 kpc, the COM of the MW+LMC system is changing as a function of radius (see Figure 13).

Mass estimates of the Milky Way: Mass estimates of the MW that use the kinematics of outer halo tracers, e.g Watkins et al. (2010), will be impacted by the barycenter motion. The reflex motion has already been shown to impact mass estimates derived from the phase space properties of halo stars using, e.g., Jeans modeling (Erkal & Belokurov, 2020; Deason et al., 2021). We argue here that the barycenter displacement (see Figure 13) must also be accounted for to accurately assess the errors in such measurements. Note that the halo response to the LMC (DM dynamical friction wake and Collective Response) as well as the LMC debris will result in asymmetric perturbations to the kinematics on the sky (G19) that will impact mass measurements differently depending on the angular location of the tracer on the sky. As shown in G19, their Figure 16, the anisotropy parameter varies across the sky (see also Erkal & Belokurov, 2020).

Radial velocities of the outer halo: The kinematics of the outer halo are predicted to exhibit a dipole in radial velocities (Garavito-Camargo et al., 2019; Petersen & Peñarrubia, 2020; Cunningham et al., 2020; Petersen & Peñarrubia, 2021; Erkal et al., 2020). Where the northern galactic hemisphere is moving away, while the southern hemisphere is moving towards us. Here we further stress that the dipole should vary as a function of Galactocentric radius and should be maximized at large distances (110-150 kpc, depending on the mass of the LMC; see Figure 14).

Here we explore how our results scale as a as a function of LMC mass and anisotropy of the MW halo. All results until now have largely focused on our fiducial model (sim #7). We have computed the BFE for 4 different LMC models [$8-25\times 10^{11}$ M_⊙]. In two different MW halo models (isotropic and radially biased), as outlined in Table 1. Our first goal is to understand the changes in the BFE, and thus the structure of the halo response.

We start by studying the effect of the LMC mass and MW model on the displacement of the COM of the inner halo with respect to the outer regions, which manifests as a reflex motion. This is shown in Figure 15. Both MW models lead to similar results (blue vs gray shaded regions). The width of the shaded line represents results from all four LMC models. The assumed LMC mass affects the amplitude of the COM displacement only at radii $>$ 30-50 kpc. The larger impact of the assumed LMC mass on the COM displacement of the outer regions of the MW halo is expected as the bound mass of the LMC was larger at infall vs. at present day (see Figure 3). This is consistent with the results of the radial velocity dipole in Figure 14.

In our study we identify a Galactocentric radius of 30 kpc as the boundary between the inner and outer halo, defined as where both the COM of the halo is no longer coincident with the disk center and where the velocity dipole is appreciable (Figures 13 and 14). We find that this value of 30 kpc is independent of LMC infall mass. This is because although the LMC models have different infall masses, they are all constrained to have similar inner mass profiles owing to constraints from the observed rotation curve (§ 2.1 and G19).

We now study the total gravitational energy (U) in the coefficients as a function of LMC mass. We use the BFEs of the MW to compute the gravitational energy in different modes as a function of LMC mass, as illustrated in Figure 16. We find that the $l=1,m=1$ modes are the most impacted by LMC mass. These modes represent the amplitude of the barycenter motion (see middle column of Figure 8). The figure shows that the ratio of the energy of the $l=1,m=1$, relative to the fiducial LMC mass model, increases as a function of LMC mass. The increase in amplitude in these terms is not linear, in particular beyond our fiducial LMC3 model, the increase is not as dramatic. The minimum value of the dipole is 20% of that of the fiducial model. For the quadrupole modes, $m=2,l=2$, the change in amplitude with mass is less pronounced.

Using the BFE of the simulated LMC we can estimate which of the observed MW satellite galaxies are bound to the LMC at the present time, accounting for the gravitational potential of the dynamical friction wake, the Collective Response, the LMC and the LMC’s DM debris. We define whether a satellite is bound or unbound using the criteria outlined in 3.2.

The solid lines in Figure 17 illustrate our main findings for six of the proposed Magellanic satellites in the literature (e.g, Jethwa et al., 2016; Sales et al., 2017; Erkal & Belokurov, 2020; Patel et al., 2020), as a function of LMC mass.

Note that we have ignored the explicit gravitational potential of the SMC, which is instead implicitly accounted for in the range of LMC masses explored. The dashed lines show the corresponding analytic calculation. assuming a rigid spherical LMC with the same mass and scale lengths as those used in the initial conditions for the G19 simulations (i.e no mass loss over time). We find that, at the present time, the only bound ultra-faint dwarf satellite of the LMC is Phoenix 2.

The SMC is barely bound to the LMC only in our fiducial model, LMC3. Increasing the mass of the LMC (e.g. LMC4) does not improve the LMC’s ability to hold on to the SMC. This is because the LMC models were constructed to match the rotation curve, which required decreasing the concentration of the LMC halo as the infall mass increased. In the simulations, this makes it easier to unbind material (see section 3.3), which is not captured in the analytic models.

We show the present day location and velocity vectors of the six most bound satellites of the LMC, in the LMC’s reference frame, in Figure 18. The locations of the proposed Magellanic satellites are largely outside the currently bound DM distribution of the LMC, but are coincident with the distribution of the LMC’s DM Debris (Sales et al., 2011; Kallivayalil et al., 2018). This suggests that some of these satellites could been bound to the LMC in the past (Patel et al., 2020).

Interestingly, all six satellites are currently moving away from the LMC. Phoenix 2 has the lowest relative radial velocity, $v_{r}\sim 6$ km/s. Patel et al. (2020) conclude that Phoenix 2 is a recently captured satellite and Jerjen et al. (2018) identified the presence of tidal arms, possibly from the interaction with the LMC. Our results are consistent with these findings, but new orbital calculations are needed that account for the LMC’s time evolving potential.

We conclude that analytic models of the LMC that ignore mass loss will overestimate the number of satellites presently bound to the LMC. However, this does not mean that these satellites were not bound to the LMC in the past and hence a time-dependent analysis remains to be done.

Within the context of gravitational dynamics, BFEs are functional expansions that are used to compactly represent the density and potential fields of collisionless gravitational systems. BFEs could be used within gravity solvers to improve the performance of galactic N-body simulations (e.g., Clutton-Brock, 1973; Weinberg, 1989; Hernquist & Ostriker, 1992; Weinberg, 1994, 1995) or to analyze the output from N-body simulations (e.g., Lowing et al., 2011). Hence there is plenty of room for using BFEs in various astrophysical regimes, including: quantifying DM halo shapes, quantifying the morphology of galactic disks and the bar (e.g., Holley-Bockelmann et al., 2005; Petersen et al., 2019), building time-dependent potentials for galaxies and galaxy clusters, computing fast orbit integration of tracers in time-dependent potentials (Garavito-Camargo, et al. in prep), and building strong gravitational lensing maps.

There are two important aspects that one needs to take into account when using BFEs: 1) choosing the right BFE; and 2) choosing the length of the expansion (this will be discussed in detail in section B). Choosing the appropriate BFE depends on the dynamical state of the system to be studied. If the system does not drastically change in morphology within a dynamical time, one can use the BFE that is built on an analytical basis (e.g., this study, Clutton-Brock, 1973; Hernquist & Ostriker, 1992; Lilley et al., 2018b, a). However, if the system is changing rapidly, the chosen BFE will not represent the system closely and hence the convergence of the expansion is not guaranteed. To alleviate this obstacle, Weinberg (1999) applied BFE to directly solve numerically the Sturm-Liouville equation in which a pair of bi-orthogonal functions are found at every time-step. Hence the system can be represented accurately at every time-step. As such, BFEs provide a unique way to characterize and simulate potentials of galaxies accurately and efficiently, enabling broad applications.