# Constraining the Milky Way Halo Kinematics via Its Linear Response to the Large Magellanic Cloud

Simon Rozier, Benoit Famaey, Arnaud Siebert, Giacomo Monari, Christophe

Pichon, Rodrigo Ibata

► To cite this version:

Simon Rozier, Benoit Famaey, Arnaud Siebert, Giacomo Monari, Christophe Pichon, et al. Constraining the Milky Way Halo Kinematics via Its Linear Response to the Large Magellanic Cloud. The Astrophysical Journal, 2022, 933, 10.3847/1538-4357/ac7139. insu-03748249

HAL Id: insu-03748249

https://insu.hal.science/insu-03748249

Submitted on 9 Aug 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Constraining the Milky Way Halo Kinematics via Its Linear Response to the Large Magellanic Cloud

Simon Rozier<fnref n="1" />, Benoit Famaey<fnref n="1" />, Arnaud Siebert<fnref n="1" />, Giacomo Monari<fnref n="1" />, Christophe Pichon<fnref n="2" />,<fnref n="3" />, and Rodrigo Ibata<fnref n="1" />

Received 2022 January 14; revised 2022 May 9; accepted 2022 May 17; published 2022 July 7

## Abstract

We model the response of spherical, nonrotating Milky Way (MW) dark matter and stellar halos to the Large Magellanic Cloud using the matrix method of linear response theory. Our computations reproduce the main features of the dark halo response from simulations. We show that these features can be well separated by a harmonic decomposition: the large-scale over/underdensity in the halo (associated with its reflex motion) corresponds to the $l=1$ terms, and the local overdensity to the $l \ge 2$ multipoles. Moreover, the dark halo response is largely dominated by the first-order *forcing* term, with little influence from self-gravity. This makes it difficult to constrain the underlying velocity distribution of the dark halo using the observed response of the stellar halo, but it allows us to investigate the response of stellar halo models with various velocity anisotropies: a tangential (respectively radial) halo produces a shallower (respectively stronger) response. We also show that only the local wake is responsible for these variations, the reflex motion being solely dependent on the MW potential. Therefore, we identify the structure (orientation and winding) of the in-plane quadrupolar ($m=2$) response as a potentially good probe of the stellar halo anisotropy. Finally, our method allows us to tentatively relate the wake strength and shape to resonant effects: the strong radial response could be associated with the inner Lindblad resonance, and the weak tangential one with corotation.

*Unified Astronomy Thesaurus concepts:* Milky Way dark matter halo (1049); Milky Way stellar halo (1060); Milky Way dynamics (1051); Analytical mathematics (38); Perturbation methods (1215)

## 1. Introduction

The nature of dark matter (DM) is certainly one of the most pressing questions of modern physics. While DM is generally assumed to consist of a cold and collisionless component of non-baryonic particles, direct evidence of its particle nature is still lacking. While a direct detection might still take a very long time, an achievable short-term goal would be to test whether galaxies do indeed reside in DM halos that are made of a self-gravitating sea of collisionless particles, interacting with baryons and with each other through gravity. A distinctive signature of such a self-gravitating sea of particles is that it reacts to external perturbations in a predictable way: this response would in principle itself leave a signature in the dynamics of the stellar halos of galaxies, and possibly in their disk dynamics too. Since current and forthcoming large surveys of the Milky Way (MW) are mapping the kinematics of its stellar halo with unprecedented accuracy, predicting and searching for such signatures is timely. Once detected, they could also provide additional information on the DM distribution in the outskirts of the Galaxy, but also unique information on the underlying *phase-space* distribution of DM, which is difficult to get access to otherwise. Hence, studying the response of the MW DM halo to external perturbations is in principle a unique way to gain insight into both the existence and nature of DM and into its phase-space distribution.

In recent years, an array of evidence has emerged to indicate that the main first-order perturbation to the MW halo comes

from the infall of the Large Magellanic Cloud (LMC), whose total mass might represent more than ~10% of the MW mass. Early analyses of its internal dynamics concluded that the mass of the LMC (with a stellar mass of $2.7 \times 10^9 M_\odot$) had to be larger than $1.5 \times 10^{10} M_\odot$, with its enclosed mass within 8–9 kpc (Schommer et al. 1992; van der Marel & Kallivayalil 2014). It soon also became clear that the LMC was most probably on its first infall toward the MW (Besla et al. 2007), meaning that its DM halo would be essentially unstripped: combined with the requirement that the LMC and Small Magellanic Cloud (SMC) have been a long-lived binary, this first infall scenario would imply a total mass of as much as $1.8 \times 10^{11} M_\odot$ (Kallivayalil et al. 2013). This mass is also well in line with expectations from abundance matching (e.g., Behroozi et al. 2013; Moster et al. 2013). This was later confirmed by a timing constraint integrating backward the orbits of galaxies currently sitting within 3 Mpc in the Local Volume as well as the relative motion between the MW and Andromeda, with the LMC mass as a free parameter, returning a very high mass of $2.5 \times 10^{11} M_\odot$ (Peñarrubia et al. 2016).

The first study of the influence of such a massive LMC on the dynamics of the Galaxy was conducted by Gómez et al. (2015), who concluded that due to the displacement of the center of the system, the inner regions of the MW were moving with a bulk velocity w.r.t. its outskirts, a reflex motion that would create a dipole in the stellar velocity field. Tentative observational evidence for this reflex motion has been recently provided by Petersen & Peñarrubia (2021) and Erkal et al. (2021). This reflex motion is expected to be accompanied by a local wake trailing behind the LMC (see, e.g., Garavito-Camargo et al. 2019), a phenomenon also tentatively detected by Conroy et al. 2021. But, as mentioned above, the actual

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France; s.rozier@astro.unistra.fr</footnote>

<footnote n="2">²Institut d'astrophysique de Paris, CNRS and Sorbonne Université, UMR 7095, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3">³IPHT, DRF-INP, UMR 3680, CEA, Orme des Merisiers Bat 774, F-91191 Gif-sur-Yvette, France</footnote>
global response of the DM halo (and subsequently, of the stellar halo) is more subtle than that, and depends on the underlying distribution of DM in phase space, as studied by, e.g., Laporte et al. (2018), Garavito-Camargo et al. (2019), Cunningham et al. (2020), Garavito-Camargo et al. (2021a), and Tamfal et al. (2021). It has, for instance, been shown that the LMC could induce an overdensity in the orbital poles of particles moving through the halo (Garavito-Camargo et al. 2021b), although it has been subsequently shown that the effect was only present for particles with low specific angular momentum (Pawlowski et al. 2021).

Searching for all those signatures obviously requires one to depart from the assumption that the Galaxy is in equilibrium, but starting from equilibrium configurations is still useful as it can serve as a basis for analytic linear perturbation theory, which represents a powerful alternative to numerical simulations. In the present paper, we apply such methods to the specific case of the response of the MW halo to the infall of the LMC. They allow us to efficiently compute the response of several possible MW models, and importantly, to gain physical insight into the various processes and resonances at play, which are obviously more difficult to decipher in N-body simulations.

The paper is organized as follows. In Section 2, we develop the essential steps of the response matrix formalism, while some analytical details can be found in Appendices A and B, some numerical details of the computation in Appendices D and E, and a validation of the implementation in Appendix D.3. Section 3 details our models for the LMC and the MW, and shows our results for the response of the MW to the LMC in a fiducial isotropic case. In Section 4, we show that the halo's self-gravity has a minor effect, which prevents us from constraining the DM phase-space structure from the observation of the stellar halo, but allows us to focus on the forced response of the stellar halo for various anisotropies. In Section 5, we discuss the dependence of our results on the details of the LMC's orbit, and we take full advantage of our method by (i) extracting meaningful information on the MW's reflex motion and local wake from separating different multipolar components, and by (ii) identifying possible resonant effects acting to build up the wake. We conclude and summarize our results in Section 6.

## 2. Linear Response Theory and the Matrix Method

In order to analytically model the response of the MW to the infall of the LMC, we resort to the linearization of the collisionless Boltzmann equation (CBE)-Poisson system (or Vlasov-Poisson), a framework known as linear response theory (Binney & Tremaine 2008). We first describe here the guiding principles of this analysis, as well as the particular method that we use to tackle this problem, called the matrix method.

### 2.1. Linear Response Theory

The MW background potential is modeled as a spherical potential $\psi_0(r)$. It basically represents its DM halo, which dominates the potential at large radii, but note that it can also include a baryonic component, and that the general method described here is not restricted to spherical potentials. Such potentials are, however, best suited as a first step for the analysis we conduct hereafter.

In the absence of any perturber, the orbits of stars and DM particles in this potential are regular, since their dynamics

derive from the spherically symmetric Hamiltonian

$$H_0(\mathbf{x}, \mathbf{v}) = \frac{v^2}{2} + \psi_0(r). \qquad (1)$$

The orbits are then best described in action-angle phase-space coordinates. The actions $J = (J_r, L, L_z)$, where $J_r$ is the radial action, $L$ the norm of the angular momentum, and $L_z$ its projection onto the z-axis, are then fully labeling regular orbits, while the canonically conjugate angle coordinates indicate where a given particle is located along its orbit.

According to Jeans' theorem (Binney & Tremaine 2008), the phase-space distribution function (DF) of a fully phase-mixed system at equilibrium is only a function of the actions: $F(J)$. Here, we define the equilibrium phase-space DF such that $F(\mathbf{x}, \mathbf{v}) dxdv$ is the mass located in the phase-space volume $dxdv$ around $(\mathbf{x}, \mathbf{v})$. As a consequence, we have that $\int dxdv F = M_{\text{tot}}$, the total mass of the distribution of matter, the response of which we choose to study. In other words, the DF $F$ represents a distribution of orbits (which could be either DM, stars, or both) in the potential $\psi_0$, and linear response theory aims at computing the self-gravitating, collisionless response of this collection of orbits when they are perturbed, in our case, by the LMC's infall. We emphasize that the DF $F$ need not self-consistently generate the mean-field potential $\psi_0$, although we will consider that specific case in Section 3.3.

The LMC is modeled as an external perturber with potential $\psi^s$, with a small amplitude compared to the potential $\psi_0$ of the MW. The MW's response is represented as a self-induced potential perturbation $\psi^s$, and is also considered as a source of perturbations on the mean-field orbits, so that the self-gravity of the response is well taken into account. In addition to these extra forces, we consider the reference frame to be accelerated, resulting in a corresponding pseudo-force term in the Hamiltonian. Our accelerated reference frame indeed follows the motion of the MW center in the asymmetric potential generated by the perturber and the response (the reflex motion, see Section 1). This acceleration corresponds to that of a test particle, initially at rest at the center of the MW, and accelerated by the gravitational influence of the total perturbation $\psi_1 = \psi^e + \psi^s$. All these effects add up to a Hamiltonian $H_0 + \Delta H$ with

$$\Delta H(\mathbf{x}, t) = \psi^e(\mathbf{x}, t) + \psi^s(\mathbf{x}, t) + \mathbf{a}_c(t) \cdot \mathbf{x}, \qquad (2)$$

where the acceleration of the MW center is given by<fnref n="4" />

$$\mathbf{a}_c(t) = -\left.\frac{\partial\psi_1(\mathbf{x}, t)}{\partial\mathbf{x}}\right|_{\mathbf{x}=0} = \int d\mathbf{x} \rho_1(\mathbf{x}, t) \frac{G}{|\mathbf{x}|^2} \mathbf{e}_r, \qquad (3)$$

with $\mathbf{x}$ centered at the MW center, $\rho_1 = \rho^e + \rho^s$ the density of the perturber and the response, $G$ Newton's constant of gravity, and $\mathbf{e}_r = \mathbf{x}/|\mathbf{x}|$.

In this reference frame, the CBE can be linearized, with the total perturbed distribution function written as $F+f$ with $f \ll F$, giving us the linearized CBE

$$\frac{\partial f}{\partial t} + \Omega \cdot \frac{\partial f}{\partial \theta} - \frac{\partial F}{\partial \mathbf{J}} \cdot \frac{\partial \Delta H}{\partial \theta} = 0, \qquad (4)$$

<footnote n="4">⁴ Note that this acceleration is different from that used in Murali (1999). Indeed, the acceleration considered there was that of the host's barycenter. This does not correspond to the acceleration of the host's density center, which we are following here because the perturber penetrates the host, creating a shift in the motion of the host's outer parts versus its inner parts.</footnote>
where $\theta = (\theta_1, \theta_2, \theta_3)$ are the angles canonically conjugated to the mean-field actions $J$, and

$$ \Omega = (\Omega_1, \Omega_2, 0) = \frac{\partial H_0}{\partial J} \qquad (5) $$

are the corresponding mean-field orbital frequencies.

Together with the Poisson equation, $\nabla^2\psi^s = 4\pi G\rho^s$ with $\rho^s = \int dv f$, this system of partial differential equations allows for the full integration of our MW model's response to the LMC.

## 2.2. The Matrix Method

To proceed, we follow Kalnajs (1977), who undertook a projection of all perturbed quantities on a bi-orthogonal basis of the potentials and densities. This practice has gained renewed interest in recent years (see, e.g., Sanders et al. 2020; Garavito-Camargo et al. 2021a; Petersen et al. 2022) for its ability to solve the Poisson equation by construction, allowing for more natural reconstructions of the gravitational potential in simulations. In the context of linear response theory, this technique can also be used to solve the Poisson equation, while the CBE is transformed into an integral equation in a linear space.

In short, and as we detail in Appendix A, we start from a bi-orthogonal basis of potential-density pairs $(\psi_n^{(p)}, \rho_n^{(p)})$, where $(p)$ typically stands for a triplet of indices $l^p, m^p, n^p$ with $l^p \ge 0$, $|m^p| \le l^p$, and $n^p \ge 0$, with the bi-orthogonality condition

$$ \int d\mathbf{x} \, \psi^{(p)}(\mathbf{x}) \rho^{(q)}(\mathbf{x}) = -\delta_p^q, \qquad (6) $$

and we define the projections of the perturbing potentials

$$ \begin{aligned} \psi^s(\mathbf{x}, t) &= \sum_p a_p(t) \psi^{(p)}(\mathbf{x}), \\ \psi^e(\mathbf{x}, t) &= \sum_p b_p(t) \psi^{(p)}(\mathbf{x}), \end{aligned} \qquad (7) $$

where the vectors $\mathbf{a}(t)$ and $\mathbf{b}(t)$, respectively, correspond to the response and external perturbation at time $t$, projected on the bi-orthogonal basis.

Then, using the bi-orthogonality condition of Equation (6), we can write

$$ a_p(t) = - \int d\mathbf{x} \int d\mathbf{v} f(\mathbf{x}, \mathbf{v}, t) \psi^{(p)*}(\mathbf{x}), \qquad (8) $$

and we can replace $f(\mathbf{x}, \mathbf{v}, t)$ with the solution of the linearized CBE (Equation (4)), assuming that the system is unperturbed at the initial time. Making use of the fact that the integration variables can be canonically changed from $d\mathbf{x}d\mathbf{v}$ to $d\mathbf{J}d\mathbf{\theta}$, we then get (see Appendix A)

$$ \mathbf{a}(t) = \int_0^t d\tau \, \mathbf{M}(t-\tau)[\mathbf{a}(\tau) + \mathbf{b}(\tau)], \qquad (9) $$

where the response matrix $\mathbf{M}(t-\tau)$ contains the full information on the gravitational dynamics in the mean-field equilibrium as a function of time delay $t-\tau$. This matrix is given by

$$ \mathbf{M}_{pq}(t) = -i(2\pi)^3 \sum_n \int d\mathbf{J} \, e^{-i\mathbf{n}\cdot\mathbf{\Omega}t} \mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}} \times \psi_n^{(p)*}(\mathbf{J})(\psi_n^{(q)}(\mathbf{J}) + \phi_n^{(q)}(\mathbf{J})). \qquad (10) $$

Some elements of this expression derive from the Fourier transform of spatial functions w.r.t. the angles (see Equation (A1a)), with the resonance vector of integers $\mathbf{n} = (n_1, n_2, n_3)$ labeling the corresponding Fourier coefficients. In more detail, $\psi_n^{(p)}$ results from the transformation of the potential basis function with index $(p)$, and $\phi_n^{(q)}$ emerges from the transformation of the $(q)$ component of the non-inertial fictitious potential, given by

$$ \phi_n^{(q)}(\mathbf{x}) = \mathbf{x} \cdot \int \frac{d\mathbf{x}}{|\mathbf{x}|^2} \rho_n^{(q)}(\mathbf{x}) e_r. \qquad (11) $$

The details of the demonstration of Equations (9) and (10) are given in Appendix A, while the particular choice of the bi-orthogonal basis for the present study (Clutton-Brock 1973) is detailed in Appendix B.

Equation (9) highlights how this method incorporates the response's self-gravity. Indeed, not only does the response $\mathbf{a}(t)$ derive from the influence of the perturber at all time steps $\mathbf{b}$, via the dynamics encoded in the response matrix, but the response at the current time also derives from the influence of the response itself at all previous time steps, as represented by the term $\mathbf{a}(\tau)$ on the right-hand side (rhs). If self-gravity could be neglected, one would only need to replace $\mathbf{a}(\tau)$ by 0 to compute the system's response to the perturber (see Section 4.1).

While the matrix method was mostly used in its frequency version (after a Laplace transform of all time-varying quantities) to detect linear instabilities in various types of self-gravitating systems (see, e.g., Zang 1976; Polyachenko & Shukhman 1981; Weinberg 1991; Vauterin & Dejonghe 1996; De Rijcke et al. 2019; Rozier et al. 2019; Breen et al. 2021), the present time version has scarcely been used, despite its ability to solve for the self-gravitating response of perturbed stellar systems (but see Seguin & Dupraz 1994; Weinberg 1998; Murali 1999, and Pichon & Aubert 2006 for noticeable, if not unique, exceptions). We therefore hope that the present study will help in reviving the interest of the community in this approach.

## 2.3. From Time Integration to Matrix Inversion

In the form of Equation (9), the problem is not yet explicitly linear, in the sense that there remains a step to directly relate the perturber $\mathbf{b}$ and the response $\mathbf{a}$ through a linear relation. Indeed, it appears that the response at the current time, $\mathbf{a}(t)$, explicitly depends on the self-induced perturbation at all previous times, as given by the $\mathbf{a}(\tau)$ term in the integrand on the rhs. To do so, we will approximate the time integral by its Riemann sum, therefore exhibiting the intrinsic linearity of the problem.

Let us assume that we aim at computing the linear response for a full period of time $[0, T]$. Let us next divide this time interval in $K+1$ steps $0 = t_0 < \dots < t_K = T$. At each step $i \ge 1$, Equation (9) can be approximated as

$$ \mathbf{a}(t_i) = \sum_{j=0}^{i-1} \Delta t \, \mathbf{M}(t_i - t_j)[\mathbf{a}(t_j) + \mathbf{b}(t_j)], \qquad (12) $$

using the rectangular rule at the lower bound of each step, with $\Delta t = T/K$. Let us now define $\mathbf{a}$ (resp. $\mathbf{b}$) as the vector built by stacking all vectors $\mathbf{a}(t_0), \dots, \mathbf{a}(t_K)$ (resp. $\mathbf{b}(t_0), \dots, \mathbf{b}(t_K)$) on top of each other. Furthermore, the matrix $\mathbf{M}$ is defined by blocks, so
that the block in the line i and column j is given by

$$ \underline{M}_{ij} = \begin{cases} \Delta t \mathbf{M}(t_i - t_j) & \text{for } j < i, \\ \mathbf{0} & \text{for } j \ge i. \end{cases} \qquad (13) $$

Here, the vectors **a** and **b** contain the information on the external and induced perturbations over the full time interval [0, T], and the matrix **M** contains the information on the system's linear dynamics over all possible time delays. With these definitions, Equation (12) can be rewritten as a matrix product as

$$ \underline{a} = \underline{M} (\underline{a} + \underline{b}). \qquad (14) $$

According to Equation (13), the matrix $I - \underline{M}$ can always be inverted, where $I$ is the identity matrix of suitable size, and the system's response over the full time interval can be computed by a simple product of matrices through

$$ \underline{a} = [I - \underline{M}]^{-1} \underline{M} \underline{b} = ([I - \underline{M}]^{-1} - I) \underline{b}. \qquad (15) $$

In this form, the problem is evidently linear, and its resolution can be summarized in the following steps: (i) take a model for the external perturber's density at each time step, and project it on the bi-orthogonal basis to get the full perturbing vector **b**; (ii) compute the full response matrix **M** following Equations (10) and (13); (iii) perform the matrix inversion and multiplication according to Equation (15) to compute the host system's response **a** at each time step. These are the steps we take hereafter to compute the MW's response to the LMC, as detailed in the next section.

One asset of response theory is to provide an explicit linear relationship via Equation (15) between the response of the MW halo, a parametric representation of the underlying property of the unperturbed equilibrium, and the properties of the perturbation. While observing the former and the latter, one can constrain the corresponding parameter, hence, e.g., probe the internal kinematics of the MW halo.

## 2.4. Two-component System

In order to describe an MW model made of stars and DM, let us introduce the linear response of a system made of two subcomponents. For that purpose, we follow the work of Weinberg (1998) in this section.

Let us consider the MW halo to be made of a dominant DM component of mass $M_{DM}$, and a light stellar component of mass $M_*$, so that $M_{MW} = M_{DM} + M_*$ and $M_*/M_{MW} = \epsilon \ll 1$. Note that we neglect here the stellar disk, which we effectively absorb within the DM component, so that $M_*$ only represents the stellar halo, with $\epsilon \sim 10^{-3}$. Let us assume that both components follow the identical density profile, which adds up to the total potential $\psi_0$, and that they are described by the DFs $F_{DM}$ and $F_*$.

$$ \left( \begin{array}{c} \underline{a}_{DM} \\ \underline{a}_{*} \end{array} \right) = \left( \begin{array}{cc} \underline{M}_{DM} & \underline{M}_{DM} \\ \underline{M}_{*} & \underline{M}_{*} \end{array} \right) \cdot \left( \begin{array}{c} \underline{a}_{DM} \\ \underline{a}_{*} \end{array} \right) + \left( \begin{array}{cc} \underline{M}_{DM} & \mathbf{0} \\ \mathbf{0} & \underline{M}_{*} \end{array} \right) \cdot \left( \begin{array}{c} \underline{b} \\ \underline{b} \end{array} \right), \qquad (16) $$

where $\underline{a}_{DM}$ and $\underline{a}_{*}$ are the responses in each subcomponent, and the same perturber **b** is applied. The respective responses of the two subcomponents therefore verify the system of coupled equations

$$ \underline{a}_{DM} = \underline{M}_{DM} (\underline{a}_{DM} + \underline{a}_{*} + \underline{b}), \qquad (17a) $$

$$ \underline{a}_{*} = \underline{M}_{*} (\underline{a}_{DM} + \underline{a}_{*} + \underline{b}). \qquad (17b) $$

Since we assumed $\epsilon \ll 1$, then the response of the stellar halo is generically negligible in mass w.r.t. that of the DM, i.e., $\underline{a}_{*} \ll \underline{a}_{DM}$, and the system of equations becomes

$$ \underline{a}_{DM} = ([I - \underline{M}_{DM}]^{-1} - I) \underline{b}, \qquad (18a) $$

$$ \underline{a}_{*} = \underline{M}_{*} [I - \underline{M}_{DM}]^{-1} \underline{b}. \qquad (18b) $$

This last equation reveals how the kinematics of the DM halo can impact the response of the stellar halo. Indeed, if $[I - \underline{M}_{DM}]^{-1}$ is significantly different from the identity matrix, then the stellar response depends on the kinematic state of the DM via its response matrix $\underline{M}_{DM}$. As we will discuss in Section 4.1, $[I - \underline{M}_{DM}]^{-1}$ describes how much the DM's self-gravity impacts the structure of the response: the larger the impact of self-gravity, the further from identity $[I - \underline{M}_{DM}]^{-1}$ is.

# 3. MW's Response to the LMC

In this section, we describe our model for the MW and the LMC, we apply the matrix formalism to that interaction, and analyze our results.

## 3.1. Models for the MW and the LMC

In order to qualitatively compare our results to those of N-body simulations from the literature, we chose to represent the MW and the LMC with models that resemble those of the fiducial simulation of Garavito-Camargo et al. (2021a). In the present case, the MW is fully described as a self-consistent spherical halo, incorporating both DM and stellar halo components. It is modeled as a Hernquist sphere with a DF from Baes & van Hese (2007) (see Appendix C for details), which self-consistently generates its total potential. The total MW mass is taken to be $M_{MW} = 1.57 \times 10^{12} M_\odot$, and the scale radius of the Hernquist profile is $a_{MW} = 40.8$ kpc. The halo is isotropic, i.e., it has a constant $\beta$ parameter set to 0, and it does not rotate. These characteristics define the mean-field potential ($\psi_0$) and phase-space DF, $F(E)$, which are used in the computation of the response matrix (Equation (10)). Note that this setup straightforwardly enters the two-component description of Section 2.4 by considering that the DFs of the DM and the stars are proportional, with $F_* = \epsilon F$ and $F_{DM} = (1 - \epsilon) F$. In that case, the response matrices and the responses are also proportional, with the same factor $\epsilon$ between the DM and stellar quantities.

The LMC is also modeled as a Hernquist sphere, with a total mass of $M_{LMC} = 1.8 \times 10^{11} M_\odot$, and a Hernquist scale radius of $a_{LMC} = 20$ kpc. Note that only the LMC density matters here, not its internal dynamics, as it is merely considered as a gravitational perturber to the MW. This density is used to construct the vector **b**.

In order to represent the LMC's infall onto the MW, we simply integrated the orbit of a particle in the aforementioned MW potential, starting from estimates of the position and velocity of the LMC at its pericenter ($\sim 50$ Myr ago):
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Evolution of the distance from the MW center to the LMC center, as given (dashed line) by Garavito-Camargo et al. (2019), (full, black line) from a leap-frog integration in a static Hernquist MW potential, and (green) from a more realistic setup described in Section 5.1. Our fiducial orbit represented in black is relatively close to that of Garavito-Camargo et al. (2019) in terms of the evolution of the LMC’s distance to the MW, which does not exclude differences in the azimuthal path of the LMC.</figcaption></figure>

$r_{\text{p,LMC}} = 48$ kpc, $v_{\text{p,LMC}} = 340$ km. s$^{-1}$ (Salem et al. 2015). We integrated this orbit using a leap-frog algorithm with 100 time steps per gigayear, and selected the portion of the orbit that has covered the infall of the LMC since 2 Gyr ago. The distance of the LMC to the MW center is represented as a function of time in Figure 1. It may seem that this trajectory is very close to that of the LMC in the simulations from Garavito-Camargo et al. (2019); however, in our case, there is no shift between the MW center of mass (COM) and its cusp in the construction of the orbit, while they represent the distance to the MW COM. Additional differences in the shape of the MW potential also result in differences in the shape of the LMC’s orbit in its plane. We keep this simple orbit as our fiducial setup, and discuss the possible influence of the MW’s motion later in Section 5.1. In that later discussion, we build the orbit represented in green in Figure 1. From our fiducial trajectory, we extracted $K+1 = 21$ equally spaced time steps to represent the LMC’s orbit, i.e., $\Delta t = 100$ Myr. We finally have a succession of mass density profiles for the LMC in the frame centered on the MW center, $\rho_{\text{LMC}}(\mathbf{x}, t_i)$ for $0 \le i \le K$.

This time interval may seem large for N-body simulations; however, it is appropriate for linear theory. Indeed, in the case of N-body simulations, the inertial motion of a particle between two time steps is a straight line of constant velocity, as if the particle were isolated. This implies large discrepancies in the particles’ orbits when the time resolution is not high enough. For the response matrix method, three key features still occur between two time steps: (i) the system’s response still follows the orbits in the mean-field potential, (ii) the system still responds to the perturber, as if it had not moved from the previous time step, and (iii) the system still responds to the response itself (the response is self-gravitating), as if it had not moved since the previous time step. These features imply that there is a much lower number of time steps required to compute the system’s response with reasonable accuracy. Typically, the matrix method requires a rough representation of the positions that the perturber takes in its motion, while N-body simulations require a fine integration of the orbits into the host.

## 3.2. Projection of the LMC onto the Basis

In order to apply the matrix formalism to the MW–LMC interaction, we need to project the perturber (the LMC) on the

bi-orthogonal basis, to compute the perturbing vector $\mathbf{b}$. To do so, we simply apply Equation (A9) with the density of the LMC (considering the LMC’s orbital plane as the $z=0$ plane) and the basis defined by Equation (A5a) and Appendix B. This gives

$$b_p(t_i) = -\int d\mathbf{x} \rho_{\text{LMC}}(\mathbf{x}, t_i) \psi^{(p)}(\mathbf{x}). \quad (19)$$

Recall that the index $(p)$ actually stands for the three integers $m^p$, $\ell^p$, and $n^p$. In our basis, $m^p$ and $\ell^p$ identify the azimuthal and polar harmonic order of the basis element, while $n^p$ identifies its radial order. Applying the projection of Equation (19) to a series of $(p)$ yields the sub-vector $\mathbf{b}(t_i)$, which is then computed at all times $t_0 < \dots < t_K$ to give the stacked vector $\mathbf{b}$.

In theory, the potential-density basis has an infinite number of elementary functions, so that an infinite number of projection coefficients is required to reconstruct the perturber with arbitrary precision. In practice, the basis should be truncated in both its angular and radial orders. These orders of truncation are parameters that should be tested for convergence to ensure the robustness of our results. In the angular direction, our fiducial choice is $\ell_{\text{max}} = 6$, while the definition of the spherical harmonics always imposes that $-\ell \le m \le \ell$. This choice is motivated by two facts: first, Garavito-Camargo et al. (2021a) show that most of the information on the MW’s response is contained within low harmonics $\ell \le 4$. Second, as is shown in Appendix D.1, in a spherical, nonrotating system, there is no coupling in the system’s response between different angular harmonics. Therefore, reconstructing the MW’s response up to $\ell = 6$ only requires us to project the LMC up to the same harmonic order. For each of these harmonics, we restrict $m$ to values such that $0 \le m \le \ell$ and that $\ell-m$ is even because we are dealing with real fields that are symmetric w.r.t. the equatorial plane. In terms of radial truncation, we found that using a fiducial maximal order of $n_{\text{max}} = 200$ represented a good balance between the accuracy of the reconstruction and the computational time. We therefore have a total of 3216 basis functions.

Figure 2 shows how the LMC’s density in its orbital plane is reconstructed once it has been projected on this truncated basis. Strikingly, the quality of the spatial reconstruction of the LMC is not constant through time. In particular, the LMC appears as a shallow angular ripple at the beginning of its infall, when it is the farthest from the MW center. This can be mainly explained by the low value of $\ell_{\text{max}}$. Indeed, the angular resolution of the basis is approximately given by $\pi/\ell_{\text{max}}$, so that when the LMC is far from the MW center, it is too small to be well resolved. We checked, however, that the LMC’s total mass is well recovered, even in the first time steps.<fnref n="5" /> By the end of the interaction, the LMC is much closer to the MW and its angular structure can be better resolved by our basis. In terms of radial reconstruction, it seems that the position and structure of the LMC are well reconstructed after the projection. This is expected, as we have used a relatively large number of radial elements.

## 3.3. Results: Response of the MW

Once the perturbing vector $\mathbf{b}$ is computed, the response is recovered by applying Equation (15). The details of the

<footnote n="5">⁵ Indeed, the information on the LMC mass is only borne by the $\ell=0$ harmonics, so that the quality of the mass reconstruction only depends on the radial truncation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Reconstruction of the LMC density after a projection on the bi-orthogonal basis. $t = 0$ Gyr corresponds to the initial conditions, 2 Gyr ago. The full circle represents the MW scale radius $a_{\text{MW}}$, while the dotted circle is the LMC scale radius, $a_{\text{LMC}}$. Although it has a limited angular resolution, the basis reconstruction marks well the position of the LMC at all time steps. The reconstruction is more accurate at late times, when the LMC is close to the MW. The mass of the LMC is reconstructed by the $l = 0$ component; therefore, it is unaffected by the limited angular resolution.</figcaption></figure>

computation of the response matrix for a spherical nonrotating system are given in Appendix D, and we used tailor-made matrix inversion and multiplication algorithms as described in Appendix E. In particular, the response matrix is given by Equations (D13–(D15)) as an integral over the two actions $\tilde{J} = (J_r, L)$, and as a sum over the corresponding in-plane resonance vectors $\tilde{n} = (n_1, n_2)$. We checked for the convergence of the following results w.r.t. our many computational parameters; in particular, we observed little variation when we doubled or halved the basis scale radius $R_b = 11 a_{\text{MW}}$, or when we increased the maximum radial Fourier number $n_{1\text{max}} = 2$. Convergence w.r.t. the number of radial basis elements is also comfortably reached.

Once Equation (15) has been applied, the response density can be reconstructed as

$$ \rho^s(\mathbf{x}, t) = \sum_p a_p(t) \rho^{(p)}(\mathbf{x}), \qquad (20) $$

where the density basis elements $\rho^{(p)}$ are defined in Appendices A and B. Figure 3 shows the response of the MW along the LMC’s infall, in terms of the relative overdensity w.r.t. the MW background density, $\rho^s/\rho_{\text{MW}} = \rho^s/\rho_0$, in the LMC’s orbital plane. This figure, and particularly its last panel, is consistent with the corresponding figure in Garavito-Camargo et al. (2021a). This consistency indicates that linear response theory is able to realistically reproduce the self-gravitating response of the DM and stellar halos to the LMC’s infall.

In Figure 3, we can identify and follow two particular features in these density maps: on the one hand, there is an overdensity that trails behind the LMC’s trajectory. This

feature emerges as a spatially large overdensity in the first ~1.6 Gyr of the interaction, when the LMC is slow and can attract some of the MW in its vicinity. Later, the LMC starts falling faster toward the MW center, and the overdensity starts moving behind the LMC’s motion, and decays because the LMC does not stay long enough to attract more material in a single place. Since this overdensity quickly disappears as the LMC moves away, it is identified as the transient response of the MW (see, e.g., Garavito-Camargo et al. 2021a), also called the LMC’s wake. On the other hand, we can identify a large-scale dipolar over/underdensity oriented in the direction of the LMC (with the overdensity opposing the LMC), which grows in amplitude all along the LMC’s infall. This feature can be associated with the shift in the dynamics of the MW center w.r.t. that of its outer parts, i.e., the MW’s reflex motion. Since our reference frame follows the MW cusp, parts of the MW that are farther away than the LMC get shifted w.r.t. the current reference frame at each time step. This explains (i) that the dipolar feature grows in amplitude all along the interaction, as the LMC roughly stays on the same side of the MW, and (ii) that the dipolar over/underdensity populates more central regions of the MW as the LMC gets closer to the center. We will discuss these features in more detail in Section 5.

## 4. Anisotropic Stellar Halo

Now that we have shown that the matrix method is able to quantitatively compute the MW halo response to the LMC, we can take advantage of its numerical efficiency to repeat the experiment on different halo models with varying velocity anisotropy. In strongly anisotropic spheres, it has been shown (see, e.g., Henon 1973; Palmer et al. 1989; Rozier et al. 2019) that linear instabilities can occur, which makes these structures
<figure><img src="image_6.png" /><figcaption>Figure 3. Reconstruction of the response of our fiducial isotropic MW DM + stellar halo in terms of relative overdensity in the LMCs orbital plane.</figcaption></figure>

unrealistic. However, we will first show that in our fiducial model there is little impact of the MW’s self-gravity on the response, in agreement with the results of Seguin & Dupraz (1994) in a similar setup. This implies that we can treat the response of the stellar halo separately, since the absence of self-gravity means that the stellar halo is essentially insensitive to the response of the DM halo. We will therefore show the forced response of models of the stellar halo with various velocity anisotropies.

### 4.1. Influence of Self-gravity

As we previously highlighted in Section 2.2, the matrix method straightforwardly allows for the measurement of the response’s self-gravity. In practice, we can compare the self-gravitating response $a_s$, as obtained from Equation (15), to the bare response $a_b$, defined by

$$a_b = M b. \qquad (21)$$

The difference between these two responses resides in the inclusion or not of the response $a$ on the rhs of Equation (9). Formally, they satisfy the relation

$$a_s = a_b + \sum_{i=2}^{\infty} M^i b. \qquad (22)$$

The bare response only includes the response to the perturber at the first gravitational order (i.e., a single application of the response matrix). The self-gravitating one further includes the response induced in the system by this first-order response (i.e., $M(M b)$), and recursively at all other orders (i.e., all $M^i b$ with

$i \ge 3$). Alternatively (see Equation (15)), $a_s$ and $a_b$ are also related through

$$a_s = [I - M]^{-1} a_b. \qquad (23)$$

It shows that the influence of self-gravity is represented by how much $[I - M]^{-1}$ deviates from the identity matrix.

Computing the bare response follows similar steps as the self-gravitating one, as described in Appendix E. Figure 4 shows the bare response $\rho^b$ (reconstructed from $a_b$ following Equation (20)) of the MW to our LMC model at the last time step, as well as the residuals when it is subtracted to the self-gravitating one. Although it is interesting to see that the region where self-gravity seems most active is within the small-scale wake, it is striking that the bare response is the dominant component in the full self-gravitating response, while further recursive applications of the response matrix (as in Equation (22)) have smaller contributions.

This has very important consequences for the way this problem can be analyzed with linear theory, as well as in the hope of constraining the phase-space structure of the MW’s DM halo. If we want to constrain the DM halo’s response, we would need this response to influence an observable population of stars, such as the stars in the stellar halo. Coming back to our model of the MW as a two-component system (Section 2.4), with a stellar and a DM halo, what we have shown is that the total response verifies $[I - M]^{-1} \simeq I$. If we exclude the possibility that the DM and stellar components conspire into such a result, it means that the influence of the DM’s self-gravity is low. Complementary results additionally show that this conclusion is still true if the DM halo is anisotropic. Following Equation (18b), this implies that the DM kinematic
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Bare MW response (top panel) and residuals when compared with the self-gravitating one (bottom panel). The scale of the residuals is about 10 times smaller than the bare response. Self-gravity is therefore responsible for less than 10% of the total response.</figcaption></figure>

state has little influence on the stellar halo response. The latter therefore mainly corresponds to the forced impact of the perturber ($a_* \approx M_* b$). As a consequence, there is little hope for strongly constraining the kinematics of the DM halo from the response of the stellar halo to the LMC.

Let us emphasize that the minor influence of self-gravity should not be regarded as a property of the MW, but as a consequence of the merger's young age. Indeed, while the ability to amplify perturbations through self-gravity is a property of the initial equilibrium, the self-gravitating wakes still take time to build up, even in a highly responsive system. Formally, the response matrix bears the information of the halo's dynamics, but only for the finite time during which we are modeling it. When the integration time $T$ is small, the matrix $[I - M]^{-1}$ cannot get very far from the identity, whatever the MW's kinematic state, and only when $T$ grows can the matrix significantly deviate from identity, and this deviation should depend on the halo's kinematics. In this work, the LMC is on its first infall, so that the interaction is too short for such self-gravitating wakes to develop, even when the MW halo is assumed to be anisotropic. In a different setup where a satellite is on a periodic orbit around the Galaxy for a long time

(an infinite amount of time, really), Weinberg (1989) shows that self-gravity has a significant influence on the Galaxy's response, even in an isotropic case.

The fact that it is possible to neglect the effect of self-gravity on the MW's response also has positive consequences for our analysis. Indeed, it means that the response of the MW stellar halo to the LMC is essentially independent of the phase-space distribution of the MW DM halo. We can therefore modify our approach, and focus on the bare response of the stellar halo only, considering the DM halo as an inert gravitational component. We perform this study for stellar halos with different orbital anisotropies in Section 4.2 below, to evaluate how the strength of the response of the stellar halo depends on its underlying phase-space distribution.

## 4.2. Response of Stellar Halos with Different Anisotropies

Starting here, and until the end of Section 5 (except briefly in Section 5.1), we shift from studying the full self-gravitating response of the MW DM and stellar halos altogether to studying the bare response of the stellar halo only. Formally, we are assuming that the response of the stellar halo is not influenced by the DM halo's self-gravity, i.e., we are replacing $[I - M_{DM}]^{-1}$ by $I$ in Equation (18b). For that purpose, we take the DM halo as a background, inert potential, and we consider the stellar halo as a low mass component responding to the potential of the MW's DM halo, but with its own density and kinematic structure. More precisely, the background potential is still the same Hernquist sphere, but now the phase-space DF $F_*(E, L)$ only represents the stellar halo, normalized with a total mass $M_{tot} = M_* = \epsilon M_{MW}$ (with $\epsilon \ll 1$). We consider that the stellar halo is distributed according to a Hernquist density, denoted $\rho_*$, with a DF also given by Equation (C1), but rescaled by the factor $\epsilon$. Note that the total mass of the stellar halo need not be specified, since (i) we consider the bare response of the stellar halo, which has a linear dependence on its mass, and (ii) all quantities we consider are relative to the initial stellar density $\rho_*$, so the linear dependence w.r.t. $M_*$ is dropped.

Now, we let the stellar halo have a different kinematic structure by changing the value of $\beta$: due to their different formation scenarios, we can expect the stars to present different kinematics from the DM. Indeed, on the one hand, the stars in the halo either come from the tidal stripping of accreted satellites, or from star formation along gas filaments connected to the halo. In both cases, their kinematics should be imprinted by the kinematics of the gas from which they were formed, which shocks and forms stars with well-ordered motions, near the bottom of the large-scale structure's potential wells. On the other hand, the DM shell crosses but never shocks, so that particle motions are less ordered, whether the DM halo comes from the first galactic gravitational collapse, the accretion of satellites, or slower accretion from DM filaments (Pichon et al. 2011; Stewart et al. 2011; Danovich et al. 2015). We therefore test two additional values of the stellar halo's anisotropy: one tangentially anisotropic ($\beta = -0.8$), and one radially anisotropic ($\beta = 0.49$), which lies near the maximum central anisotropy allowed for a Hernquist sphere (see An & Evans 2006) and is closer to the stellar halo anisotropy measured by Bird et al. (2019). Again, the anisotropy parameter is independent of the distance to the center, and no net angular momentum is introduced. In both cases, such an anisotropic distribution may also describe a subpopulation of
<figure><img src="image_8.png" /><figcaption>Figure 5. Response of the MW stellar halo when different anisotropies of the halo are considered: tangential anisotropy ($\beta = -0.8$, top panel) and radial anisotropy ($\beta = 0.49$, bottom panel). Only the bare response of the stellar component is computed ($\rho_{*}^{b}$), and is represented relative to the initial stellar density ($\rho_{*}$). The response of the MW stellar halo is much stronger when it is radially anisotropic.</figcaption></figure>

the stellar halo, e.g., stars accreted from a single merger event (Belokurov et al. 2018), or from a collection of satellites (Riley et al. 2019).

Figure 5 presents the results of this experiment: it shows the relative overdensity $\rho_{*}^{b}$ (the bare response) of the stellar halo at the present time in the orbital plane of the LMC for both values of the stellar halo's anisotropy. These panels clearly show that the response is much stronger when the stellar halo is radially anisotropic than when it is tangentially anisotropic. In more detail, it appears that the large dipole is essentially unaffected by anisotropy, whereas the small-scale wake is strongly amplified in the radially anisotropic sphere, and strongly depleted in the tangentially anisotropic system. Finally, in the radially anisotropic system, a small-scale overdensity appears in the region of space preceding the LMC in its orbit, while this region of space presents a small-scale underdensity in the tangentially anisotropic stellar halo.

## 5. Discussion

The results of Section 3 showed that the matrix method is able to reproduce N-body experiments, and those of Section 4 clearly show that the stellar halo response depends on its initial anisotropy. Let us now discuss how these results depend on the specific orbit considered for the LMC, and how they provide physical insight into the MW's response, as well as quantitative constraints on the stellar halo's anisotropy.

### 5.1. Influence of the LMC's Orbit

We now explore the influence of the specific LMC orbit we chose for our previous computations. Indeed, the orbit was computed by integrating backward a point mass in the MW static Hernquist potential, which makes it unrealistic in several ways. First, the LMC is not a point mass, and its orbit is influenced by the fact that the LMC and the MW are two inter-penetrating extended objects. More importantly, the MW is not static, and its intricate reflex motion strongly affects the LMC's dynamics in the frame centered on the MW's cusp.

In order to build a more realistic orbit, we integrated the motion of the LMC-MW system using a leap-frog integrator with the following prescriptions:

1. The motion of the MW is represented by that of its cusp: at each time step, the MW's acceleration is that of a point mass in the LMC's potential, so that

$$a_{\text{MW}} = -\nabla\psi_{\text{LMC}}(x_{\text{MW}}). \quad (24)$$

In this expression, $\psi_{\text{LMC}}$ represents the current LMC potential, and in particular, it depends on its current position.

2. The motion of the LMC is computed within the Hernquist MW potential, but centered on its cusp as determined by the previous step. The acceleration is then

$$a_{\text{LMC}} = -\frac{1}{M_{\text{LMC}}} \int d^3x \nabla\psi_{\text{MW}}(x) \rho_{\text{LMC}}(x). \quad (25)$$

Note that these rules differ from those of Gómez et al. 2015, who also consider the LMC to be a point mass in the MW potential to compute $a_{\text{LMC}}$. The present prescriptions were chosen because they more closely match the physical processes at play, as well as the framework of the matrix method. Indeed, choosing the MW cusp as the reference for computing forces is relevant for two reasons. On the one hand, since the satellite is at its first infall, it is always sensitive to more and more central parts of the MW—the region of the cusp—while the MW's outskirts act as a spherically symmetric shell with no resultant gravitational influence. In the central region, the potential is still that of a Hernquist sphere, centered on the cusp; hence, we take that potential as that responsible for the acceleration of the LMC. On the other hand, the matrix method computes the response of the MW in the reference frame of its cusp, so this refined prescription for the LMC's orbit is more adapted to it.

Note that this prescription does not conserve momentum, as the reciprocal forces applied by one object on the other are not equal. However, it appears that it still produces quasi-periodic trajectories. Using the same values of the LMC pericentric radius and velocity as in Section 3.1, we obtain the orbit shown in Figure 1. This new orbit is consistent with the new
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Density response, in the orbital plane, at the final time step ($t_{20} = 2$ Gyr) of the fiducial MW ($eta = 0$) to the LMC following the orbit described by Equations (24) and (25). The MW's response is only weakly impacted by the details of the LMC's orbit at the early stages of its infall.</figcaption></figure>

prescriptions we applied to the problem: since the MW is now moving in the LMC's potential, the LMC should come from farther away in order to produce the same pericentric distance.

We then computed the response of the fiducial MW ($\beta = 0$) to this new orbit for the LMC. Let us emphasize that this does not require us to recompute the response matrix, but only to apply Equation (15) to the new perturbing vector. In order to compare the present case with the response in Section 3, we consider here the full self-gravitating response of the DM + stellar halo to the LMC on its new trajectory. The results are shown in Figure 6. Comparing this map with the last panel of Figure 3, it appears that the MW's response is weakly sensitive to the details of the LMC's orbit. Only at large radii, in the tail of the local wake (around $(x, y) = (300 \text{ kpc}, -100 \text{ kpc})$) and the orientation of the dipolar component can one notice small deviations in the response's shape. This result hinders the possibility of setting strong constraints on the LMC's orbit far in the past from the sole study of its present influence on the MW.

## 5.2. The Reflex Motion is Pure Potential

Let us now decompose the MW stellar halo's response to the LMC, and disentangle the physical processes responsible for its characteristic shape. In Appendix D.1, we show that the linear response of a spherical, nonrotating system involves no couplings between different $m$ and $\ell$ orders in the multipolar expansion of the perturber. This is particularly highlighted by the $\delta_{m}^{m\ell}$ and $\delta_{\ell}^{m\ell}$ term in Equation (D13). The response can therefore be decomposed into a sum of separate harmonic terms. Here, we focus on the dipolar over/underdense pattern, while the local wake is studied in Section 5.3.

Figure 7 shows the bare response of the MW halo to the $\ell=1$ components of the LMC. It appears that this harmonic fully accounts for the effect of the reflex motion of the MW: it presents the same dipolar feature, with similar amplitude and orientation. As will be shown further in Figure 9, the other harmonics do not contribute to that component in the response, but to the overdense wake trailing behind the LMC's trajectory.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Linear response of the MW stellar halo to the dipolar terms in the perturber. This response includes both the $\ell=1, m = \pm 1$ terms. These terms clearly account for the large-scale over/underdense pattern related to the reflex motion of the MW halo.</figcaption></figure>

Now, let us evaluate the bare response of stellar halos with strong anisotropy to the dipolar components of the LMC. The results for the radially and tangentially anisotropic halos are shown in Figure 8. One can clearly see the similarity between the present maps and Figure 7: both the amplitude of the dipolar response and its shape do not seem to depend on the velocity anisotropy of the MW stellar halo. This implies that the MW reflex motion is independent of its internal kinematics, and only depends on its potential and the characteristics of the LMC's orbit.

Note that our ability to study this dipolar term using the matrix method was provided by a simple yet critical improvement, i.e., by considering the motion in the reference frame of the density maximum instead of the COM. Indeed, previous studies (see related discussions in Seguin & Dupraz 1994; Murali 1999) considered the COM as the reference frame for matrix computations, so that the density maximum gets displaced from the center of the reference frame. This has two unfortunate consequences on the ability of the matrix method to accurately reconstruct dipolar effects. On the one hand, the central displacement appeared as a rather large amplitude dipolar pattern located close to the center. This made it difficult to reconstruct it using bi-orthogonal basis functions, and all the more given that there is a large contrast between the central density and its close vicinity (e.g., in a cuspy system). This central displacement could also dominate over the more subtle effects related to the host's reflex motion, which we study in more detail here. On the other hand, the matrix method makes the assumption that the potential center is at the reference frame center, any deviation being considered as a small perturbation. When the density maximum is displaced, the linear hypothesis may be violated, and all the more so for a cuspy system. Overall, working in a reference frame that follows the density maximum at each time addresses some significant problems. The problem of the density maximum being off-centered is also a concern with N-body codes using basis function expansions (e.g., the so-called self-consistent field codes). This problem is usually dealt with by recentering the origin of the expansion at the density maximum of the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Response of the MW stellar halo to the $l=1$ components of the LMC, when different anisotropies of the halo are considered: tangential anisotropy ($\beta = -0.8$, top panel) and radial anisotropy ($\beta = 0.49$, bottom panel). The dipolar response of the MW stellar halo is essentially insensitive to the anisotropy.</figcaption></figure>

particles (see, e.g., Choi et al. 2007; Meiron et al. 2014; Petersen et al. 2022), but these codes may also benefit from the method we developed here to follow more closely the density maximum and correct the force calculations accordingly.

### 5.3. Local Wake as a Probe of the Halo Anisotropy

Let us now focus on the $l \ge 2$ part of the response. The middle panel of Figure 9 shows the response of the isotropic MW to the $l \ge 2$ components of the LMC. Clearly, this part of the response encompasses the local overdensity trailing behind the LMC. One can also see traces of a slight overdensity on the other side of the MW, as well as underdense regions between these two overdensities. However, these last two features are much less significant in these regions than the dipolar pattern induced by the reflex motion. In addition to these spatial shapes, and in contrast to the behavior of the $l=1$ response, the $l \ge 2$ response is transient, in the sense that it reaches a maximum amplitude while the LMC is around its apocenter, and starts decaying while the LMC falls toward the MW. We,

therefore, suggest that the *wake* should be identified with the $l \ge 2$ part of the response. This provides a more uniform definition than the visual characterization previously used in the literature. We note that this definition relies on the fact that different $l$ harmonics are gravitationally decoupled; therefore, it may be ineffective for nonspherical or rotating models.

If we consider Figure 9 altogether, we can evaluate the influence of the stellar halo anisotropy on the aspect of the local wake. Clearly, its amplitude is increased in the radially anisotropic halo, and depleted in the tangentially anisotropic one. The wake’s structure is also different: in the tangential halo, the wake very closely follows the LMC in its orbit, while there is an underdensity ahead of the satellite. In the radial case, however, the part of the wake trailing behind the LMC does so from a distance, and is accompanied by an important overdensity directly ahead of the satellite in its orbit.

As a consequence, it appears that the anisotropy of the stellar halo can be probed by measuring the structure and amplitude of the wake in the stellar density of the halo. Let us be more specific, and focus on a possibly measurable feature: the quadrupolar ($m = \pm 2$) component of the stellar distribution in the orbital plane of the LMC. Indeed, this quantity could be recovered as a Fourier transform of the stellar halo density along the angular direction in that plane. Notice that, in our response, these terms involve all even $l \ge 2$ harmonics with $m = \pm 2$. Figure 10 shows the quadrupolar response of stellar halos with different anisotropies. Not only does the amplitude of this term vary with anisotropy, but more importantly, it is more wound in the tangential halo than in the radial one, which entails different orientations in the central parts. Within ~50 kpc, the positive part of the $m = \pm 2$ response is nearly aligned with the position of the LMC in the case of a tangential halo, while it is largely misaligned from that position if the halo is radial. Farther away, the orientation of this harmonic evolves at a greater pace in the tangential case, so that in all cases, its orientation at large distances coincides with the initial location of the LMC.

The exact angles between this quadrupolar response and the LMC are likely to also depend on the MW potential and the details of the LMC’s orbit. However, there is little doubt that this dependence of the orientation on anisotropy will remain in other models of the MW–LMC interaction. Furthermore, we also expect that other kinematic features of the stellar halo (e.g., rotation or different distributions of the anisotropy) would imprint its quadrupolar response.

### 5.4. Frequencies in the Response

A possible advantage of the matrix method, compared to N-body simulations, is to interpret the MW’s response as the combined action of orbital frequencies. Indeed, since the bare response of the stellar halo merely involves the application of the response matrix to the external perturber (see Equation (21)), it can be straightforwardly decomposed into the sum of the influence of different combinations of orbital frequencies, corresponding to each value of the Fourier numbers $n$ (also called resonance vector) in Equation (10)<fnref n="6" />. As we will see, comparing their influence yields very distinct

<footnote n="6">⁶ This analysis is coined the *restricted matrix method* in Breen et al. (2021). In the present case, the absence of self-gravity makes the interpretations easier because the resonant terms are directly additive.</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 9. Linear response of the different MW models (tangentially anisotropic, isotropic, and radially anisotropic, respectively, in the left, middle, and right panels) to the $ℓ ≥ 2$ terms in the perturber. These terms clearly account for the overdense wake trailing behind the LMC.</figcaption></figure>

pictures in the radially and tangentially anisotropic models of the MW.

The influence of each **n** component of the matrix on the perturber **b** can be interpreted as follows. Let us decompose the perturber spatially in each of its azimuthal harmonics *m*, then further decompose one of these spatial harmonic components through a time Fourier transform with a spectral continuum of frequencies $\omega_p$. The resulting pattern appears as an *m*-armed feature in the *z* = 0 plane, with pattern speed $\Omega_p = \omega_p/m$ (see, e.g., Tremaine & Weinberg 1984). A single of these components therefore follows

$$ \mathbf{b}_m(t) = \mathbf{b}_{m0} e^{-im\Omega_p t}. \qquad (26) $$

When a single **n** component of the response matrix is applied to this periodic perturber, it yields the bare response (see Equations (9) and (10))

$$ \mathbf{a}_m(t) = \int_0^t d\tau \int d\mathbf{J} e^{-i \mathbf{n} \cdot \Omega_t} e^{i(\mathbf{n} \cdot \Omega - m\Omega_p)\tau} \mathbf{N}_n(\mathbf{J}) \mathbf{b}_{m0}, \quad (27) $$

where the matrix $\mathbf{N}_n(\mathbf{J})$ is defined in terms of its components by

$$ \mathbf{N}_{n,pq}(\mathbf{J}) = -i(2\pi)^3 \mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}} \psi_n^{(p)*}(\mathbf{J}) (\psi_n^{(q)}(\mathbf{J}) + \phi_n^{(q)}(\mathbf{J})). \quad (28) $$

Now, the integral over $\tau$ in Equation (27) can be straightforwardly performed, giving

$$ \mathbf{a}_m(t) = t \int d\mathbf{J} e^{-i \frac{n \cdot \Omega + m \Omega_p}{2} t} \text{sinc}\left(\frac{n \cdot \Omega - m \Omega_p}{2} t\right) \times \mathbf{N}_n(\mathbf{J}) \mathbf{b}_{m0}. \quad (29) $$

This expression helps us identify the orbits, which will have a large impact on the response in terms of their orbital frequencies. Indeed, when the phase space is scanned by the integral over $\mathbf{J}$, the sinc function acts as a window that selects the orbits for which $\mathbf{n} \cdot \mathbf{\Omega} - m\Omega_p \simeq 0$, and all the more that time gets large. This implies that the component we selected in the perturber (with spacetime structure given by *m* and $\omega_p$) has a larger impact on regions of phase space where $\mathbf{n} \cdot \mathbf{\Omega} - m\Omega_p \simeq 0$ for one of the Fourier sets **n**. This motivates calling the Fourier numbers **n** a resonance vector, and considering the resonance

condition

$$ \mathbf{n} \cdot \mathbf{\Omega} = m\Omega_p \qquad (30) $$

in the following discussion. Additionally, this model gives us insight into how the duration of the merger impacts the response. At early times, the frequency spectrum of the perturber is rather broad, and the width of the sinc window is too, so that many different regions of the host's frequency distribution can interact with the perturber. At later times, the perturber's frequency spectrum narrows down to perhaps a single or a set of frequencies (because it follows a quasi-periodic orbit). Besides, the width of the sinc window gets narrower too. Therefore, we can expect resonances to play a more prominent role on longer timescales.

In order to quantify the importance of resonant effects in the development of the wake, we follow Garavito-Camargo et al. (2021a) and compare the gravitational potential energy of its different resonant components. We first separate the harmonic terms by values of $ℓ$ and *m*, and consider the *bare response* of the stellar halo when only a selection of resonant vectors is considered in the sum of Equation (10). Formally, we are computing

$$ \underline{\underline{a}}_{lm}^{n} = \underline{\underline{M}}_{lm}^{n} \underline{\underline{b}}_{lm}, \quad (31) $$

where the subscript *lm* means that only a single harmonic is considered in the perturber and the matrix, and the superscript **n** that only a single resonant term is considered in the matrix. For ease of presentation, we will always implicitly consider $\pm m$ terms altogether (they are overdetermined because the densities and potentials are real functions), as well as $\pm \mathbf{n}$ terms (they belong to the same family of resonances). Once the different response vectors have been computed, we consider the gravitational potential energy of the resulting mass distribution (Binney & Tremaine 2008),

$$ W_{lm}^{n}(t) = \frac{1}{2} \int d\mathbf{x} \rho_{lm}^{n}(\mathbf{x}, t) \psi_{lm}^{n}(\mathbf{x}, t). \quad (32) $$

By developing the response density and potential on a bi-orthogonal basis, we can easily show that this energy is equal
<figure><img src="image_13.png" /><figcaption>Figure 10. Linear response of the different MW models (tangentially anisotropic, isotropic, and radially anisotropic, respectively, in the left, middle, and right panels) to the m = 2 terms in the perturber. The orientation of this quadrupolar response can represent a crucial probe of the stellar halo’s anisotropy.</figcaption></figure>

to

$$W_{\ell m}^{n}(t) = -\frac{1}{2} \sum_{k} |a_{\ell mk}^{n}(t)|^2, \quad (33)$$

where k identifies the radial basis element.

In the tangential case, it appears that several combinations of frequencies add up to produce the final wake, the most important of which being the corotation frequency. The top panel of Figure 11 shows minus the gravitational potential energy at the final time step, $-W_{\ell m}^{n}(t_{20})$, of the 8 resonances with the largest energies among the 183 possible resonances with $\ell \ge 2$, assuming that the stellar halo has a mass of $M_{*} = 10^{-3} M_{\text{MW}}$. It appears that several resonances with $\ell = 2$ play an important role, among which the $(\ell, m) = (2, 2)$, $(n_1, n_2) = (0, 2)$ resonance dominates$^7$. Beyond the $\ell = 2$ terms, the largest energy is reached by the $(\ell, m) = (3, 3)$, $(n_1, n_2) = (0, 3)$ resonance. We associate those $(n_1, n_2) = (0, m)$ combinations with corotating orbits, corresponding to stars that rotate together with the perturber, at the same orbital rate $\Omega_2 = \Omega_p$ (see Equation (30)). In the bottom panel of Figure 11, we show the final wake of the MW when only the $(n_1, n_2) = (0, m)$ terms are kept in the response, for all $m \ge 2$. It appears that the contribution of this frequency accounts for the bulk of the wake, with the same shape, the same position but a lower amplitude. On top of it, a number of smaller contributions from other frequencies add up to form the full wake of Figure 9. It should be noted that the corotation frequency has a relatively high value, which could be related to the fact that (i) the wake is able to move fast and to follow the perturber closely, and it is therefore more wound around the MW, (ii) the response dissipates efficiently once the LMC enters the fastest portion of its orbit, which produces a rather shallow wake at the last time step.

In the radial halo, the wake is clearly dominated by the contribution of the inner Lindblad resonance. In the top panel of Figure 12, we show the counterpart of the top panel of Figure 11 for the radial halo. This diagram shows that the $(\ell, m) = (2, 2)$, $(n_1, n_2) = (-1, 2)$ resonance has the largest potential energy, and dominates the other terms more clearly than in the tangential halo. This resonant vector is characteristic

<figure><img src="image_14.png" /><figcaption>Figure 11. Resonant contributions to the wake in the bare response of the tangentially anisotropic stellar halo. The top panel represents minus the gravitational potential energy, $-W_{\ell m}^{n}$, of the main resonant terms, assuming that the stellar halo has a mass of $M_{*} = 10^{-3} M_{\text{MW}}$. Among the 183 resonances considered in Section 5.3, these resonant terms are the eight ones with the largest energy. The $(\ell, m) = (2, 2)$ corotation (i.e., with $(n_1, n_2) = (0, 2)$) is dominant. Among all terms with $\ell \ge 3$, the $(\ell, m) = (3, 3)$ corotation (with $(n_1, n_2) = (0, 3)$) is also dominant. In the bottom panel, we represent the response when only the corotation terms are included, i.e., $n = (0, m)$. The response’s shape is similar to the full wake, but its amplitude is much lower: other combinations of frequencies also make an important contribution to the response.</figcaption></figure>

of the inner Lindblad frequency, $\Omega_{\text{ILR}} = \Omega_2 - \Omega_1/2$, because it corresponds to $\Omega_{\text{ILR}} = \Omega_p$ in Equation (30). Interestingly enough, another member of this family is present among the

<footnote>$^7$ We remind the reader that this implicitly includes all terms with $m = \pm 2$ and $(n_1, n_2) = \pm (0, 2)$.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Resonant contributions to the wake in the bare response of the radially anisotropic stellar halo. The top panel is the counterpart of that in Figure 11. The $(\ell, m) = (2, 2)$ inner Lindblad resonance (i.e., with $(n_1, n_2) = (-1, 2)$) is dominant. The $(\ell, m) = (4, 4)$ inner Lindblad resonance (i.e., with $(n_1, n_2) = (-2, 4)$) is also present among the eight resonant terms with the highest energy. In the bottom panel, we represent the response when only those terms are included, i.e., $\mathbf{n} = (-m/2, m)$. The response's shape is similar to the full wake, and its amplitude is about 50% of the wake: the contribution of other frequencies is less important than in the tangential case.</figcaption></figure>

eight most energetic resonant terms, for $(\ell, m) = (4, 4)$, $(n_1, n_2) = (-2, 4)$. In the bottom panel of Figure 12, we show the final wake of the MW when only $(n_1, n_2) = (-m/2, m)$ terms are kept in the response, for $m=2, 4$. This fraction of the response accounts for most of the amplitude of the full wake shown in Figure 9. The fact that this specific combination of frequencies dominates tends to indicate that instead of attracting particles that can move with it, the LMC rather attracts orbits that can precess with it. This $\Omega_{ILR}$ frequency has a smaller value than the orbital frequency $\Omega_2$, which could explain the fact that (i) the wake is relatively slow and trails behind the LMC from a distance, which implies that it is less wound, and (ii) the overdensity that appears at early times, when the LMC is slow and able to attract a lot of material, never dissipates, so that the final wake has a large contribution from this epoch.

The influence of resonant frequencies should, however, be interpreted with care, for several reasons. First, the building blocks of the matrix method are orbits and frequencies, which may overemphasize the role of resonances in every linear problem. For example, dominant resonant effects can also be identified when the same analysis is applied to the reflex motion discussed in Section 5.2, although it does not appear as a resonant process. Second, the resonant interpretation is based

on the assumption that the frequency spectrum of the perturber is narrow, which is obviously not the case because the LMC is integrated for just a fraction of its orbit. Finally, the MW's response is only integrated for 2 Gyr, which is quite short compared to the orbital times in the stellar halo, whereas the effect of resonances could only truly dominate the MWs response over secular timescales. On shorter timescales, their impact is broadened by a width given the timescale of the interactions (see Equation (29)).

To summarize these arguments, we propose the following interpretation, also consistent with the results of Seguin & Dupraz (1994) and Weinberg (1989). In the very early stages of the interaction ($t \lesssim 1$ Gyr), the stars in the MW halo behave like pure particles and are attracted by the gravitational field of the slowly moving LMC. In an idealized future, the LMC will have made several rotations around the MW in its orbit, and the commensurabilities between the orbits of stars in the MW and the LMC's orbit could be the main driver of the MW stellar halo's response. In the meantime, the situation is mixed: the initial impulsive effects are still important, but the stars also started to be influenced by the frequency structure of the LMC's motion.

## 6. Conclusion

Recent photometric and spectroscopic surveys focusing on the MW's stellar halo have shown that it is dynamically perturbed by the influence of the LMC. These observations require detailed modeling of this infall, in order to extract constraints on its characteristics: density and kinematics of the MW's DM and stellar halos, mass, and orbit of the LMC. While this problem is mostly studied by means of N-body simulations, here we took a complementary approach relying on the matrix method from linear response theory, which yields an explicit estimator of the impact of the LMC on the structure of the stellar and DM halos. We first showed that this analytical tool is able to quantitatively reproduce the results of N-body simulations, opening the possibility of framing the match to the observed response of the stellar halo as an inverse problem for its internal kinematics. In addition, the matrix method provides interesting insights into the physical processes at play, which informs us on the effect of minor mergers on galaxies in general. Our main results are the following:

(i) At that stage of the interaction, the response's self-gravity has little influence. As a consequence, the stellar halo can be considered as a massless component, and its response is weakly sensitive to the DM halo kinematics.

(ii) The response of the MW stellar halo is strongly dependent on its internal kinematics. Here, we focused on its anisotropy, but it is expected that global angular momentum would also impact its response.

(iii) The MW's response is weakly dependent on the details of the LMC's orbit. We compared the response when the LMC's orbit is integrated in a static MW potential, with that when the MW center is also allowed to move in the LMC's potential, yielding very similar results in the density response of the MW halo.

(iv) The reflex motion of the MW corresponds to the dipolar $\ell=1$ part of the response. Its structure is very weakly sensitive to the internal kinematics of the MW halo, and mainly depends on its potential.
(v) The local wake, as visually identified in simulations, corresponds well with the $l \ge 2$ part of the response. We suggest that in a spherical, nonrotating model, the local wake should be formally defined based on these $l \ge 2$ terms. We found that its structure is strongly dependent on the stellar halo's anisotropy, both in amplitude and in spatial distribution. Typically, a tangentially anisotropic halo produces a shallow wake and an underdensity ahead of the LMC in its orbit, while a radially anisotropic halo produces a strong wake and an overdensity ahead of the LMC.

(vi) The winding and orientation of the quadrupolar in-plane response ($m=2$) of the halo represent a novel probe of its anisotropy. Typically, a tangentially anisotropic halo produces a tightly wound $m=2$ response where the overdensity is nearly aligned with the LMC within 50 kpc of the MW, while it is more loosely wound and rather misaligned with the LMC in the radial halo.

(vii) Resonances seem to matter, in particular in the radially anisotropic case with radial orbits at the inner Lindblad resonance with the LMC.

These first results call for further work relying on the linear response formalism. On the one hand, the structure and kinematics of the MW could be modified, in order to account for the details of the DM halo density profile (e.g., different central and outer slopes, triaxiality, or clumpiness of the halo), and of the complexity of the stellar halo kinematics (e.g., global angular momentum, spatial distribution of the anisotropy). In particular, it is possible that even a mild level of rotation could significantly impact the MW halo's response. On the other hand, the tidal evolution of the LMC could be included in the model, in order to evaluate how the evolution of its mass influences the MW's response. Surely, such studies will help refine our models of the LMC infall.

The method developed here should also be adapted to disk-halo configurations, in order to allow for the detailed analytical study of the response of the MW disk-halo system to the accretion of the Sagittarius dwarf galaxy. In this merger, the effect of the self-gravity of the DM halo could potentially play a much more important role because of the multiple wraps that the Sagittarius dwarf has already made around the MW halo. More generally, it would be important to investigate how and when self-gravity becomes important, depending on the kinematic structure of the host galaxy.

We thank the anonymous referee for their helpful suggestions. This work has been supported by the Agence Nationale de la Recherche (ANR project GaDaMa ANR-18-CE31-0006). It has also received funding from the project ANR-19-CE31-0017 and from the European Research Council (ERC grant agreement No. 834148). We used the High Performance Computing facility of the Observatoire astronomique de Strasbourg for our computations, and we thank Mathieu Misslin for running it smoothly.

The data underlying this article is available through a reasonable request to the corresponding author. The code is distributed on GitHub at the URL: https://github.com/simrozier/LiRGHaM.

# Appendix A
## Derivation of the Matrix Method

Starting from the linearized CBE (Equation (4)), let us follow a path similar to that of Murali (1999) and Pichon & Aubert (2006) to derive the response matrix formalism in its temporal version.

### A.1. Angular Fourier Transform

First, let us expand all perturbed quantities of Equation (4) in a Fourier series of the angles. Since each quantity should be $2\pi$ periodic in the angles, these series can be written as

$$f(\theta, \mathbf{J}, t) = \sum_{n} f_n(\mathbf{J}, t) e^{in\cdot\theta}, \quad (A1a)$$

$$f_n(\mathbf{J}, t) = \int \frac{d\theta}{(2\pi)^3} f(\theta, \mathbf{J}, t) e^{-in\cdot\theta}, \quad (A1b)$$

and similarly for the perturbing Hamiltonian

$$\Delta H(\theta, \mathbf{J}, t) = \sum_{n} \Delta H_n(\mathbf{J}, t) e^{in\cdot\theta}. \quad (A2)$$

Here, $n \in \mathbb{Z}^3$ is the triplet labeling each Fourier coefficient. Multiplying Equation (4) by $e^{-in\cdot\theta}$ and integrating over the angles, each Fourier component separately satisfies

$$\frac{\partial f_n}{\partial t} + i \mathbf{n} \cdot \Omega f_n = i \mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}} \Delta H_n. \quad (A3)$$

Now, Equation (A3) takes the form of an integro-differential equation on $f_n$. Assuming that the system is unperturbed at the initial time, i.e., $f_n(t=0)=0$, the solution satisfies the integral equation

$$f_n(\mathbf{J}, t) = i \mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}} \int_0^t d\tau \Delta H_n(\mathbf{J}, \tau) e^{-i \mathbf{n} \cdot \Omega (t-\tau)}, \quad (A4)$$

where $\Delta H_n$ itself depends on $f_n$ through the Poisson equation.

### A.2. Basis Function Expansion

To make this dependence explicit, we project the perturbing quantities on a bi-orthogonal basis of potentials and densities. We can assume that this basis takes the form

$$\psi^{(p)}(\mathbf{x}) = \psi_{lmn}(r, \theta, \phi) = Y_{\ell}^{m}(\theta, \phi) U_{n}^{\ell}(r), \quad (A5a)$$

$$\rho^{(p)}(\mathbf{x}) = \rho_{lmn}(r, \theta, \phi) = Y_{\ell}^{m}(\theta, \phi) D_{n}^{\ell}(r), \quad (A5b)$$

where a given basis element is characterized by three indices, $\ell \ge 0$, $|m| \le \ell$, and $n \ge 0$, and $Y_{\ell}^{m}$ is a spherical harmonic following the normalization convention

$$\int d\theta \, d\phi \, \sin(\theta) |Y_{\ell}^{m}(\theta, \phi)|^2 = 1. \quad (A6)$$

In Equation (A5a), we also introduced the radial parts of the respective bases, $U_n^\ell$ and $D_n^\ell$, which are normalized so that (see also Equation (6))

$$\int d\mathbf{x} \, \psi^{(p)}(\mathbf{x}) \rho^{(q)}(\mathbf{x}) = -\delta_{p,q}. \quad (A7)$$

Several choices for these radial functions are given in the literature, e.g., by Clutton-Brock (1973), Fridman et al. (1984), Hernquist & Ostriker (1992), Bertin et al. (1994), Zhao (1996),
Rahmati & Jalali (2009), and Lilley et al. (2018). Unlike that required to fully reconstruct the perturber and the host (see, e.g., Weinberg 1999), the matrix method only demands an accurate reconstruction of the perturber, so that there is no need in the present study for a basis with cuspy elements. In Appendix B, we briefly describe our choice of basis, which is that of Clutton-Brock (1973) and has a Plummer profile (Plummer 1911) as the first element.

We define the projections of the perturbing potentials $a_p$ and $b_p$ so that (see also Equation (7))

$$ \psi^s(\mathbf{x}, t) = \sum_p a_p(t) \psi^{(p)}(\mathbf{x}), \quad (A8a) $$

$$ \psi^e(\mathbf{x}, t) = \sum_p b_p(t) \psi^{(p)}(\mathbf{x}). \quad (A8b) $$

Using such an expansion, it becomes clear that the purpose of the matrix method will be to compute the coefficients $a_p(t)$ (the system's response), given a certain set of coefficients $b_p(t)$ (the external perturber). We can now derive the response equation that relates these quantities.

### A.3. The Response Matrix

Using the bi-orthogonality condition, we can invert Equation (A8a) to get

$$ a_p(t) = - \int d\mathbf{x} \, \rho^s(\mathbf{x}, t) \psi^{(p)*}(\mathbf{x}). \quad (A9) $$

Since the DF perturbation is related to the density response in the system through

$$ \rho^s(\mathbf{x}, t) = \int d\mathbf{v} f(\mathbf{x}, \mathbf{v}, t), \quad (A10) $$

Equation (A9) can be developed as Equation (8).

Now, the DF perturbation, $f$, can itself be developed in angular Fourier elements as in Equation (A1a), and the integration variables can be canonically changed from $d\mathbf{x}\mathbf{v}$ to $d\mathbf{J}\mathbf{d}\theta$, with a Jacobian equal to 1 owing to phase-space volume conservation. We then have

$$ \begin{aligned} a_p(t) &= -\sum_n \int d\mathbf{J} f_n(\mathbf{J}, t) \left( \int d\mathbf{\theta} e^{-i\mathbf{n}\cdot\mathbf{\theta}} \psi^{(p)}(\mathbf{\theta}, \mathbf{J}) \right)^* \\ &= -(2\pi)^3 \sum_n \int d\mathbf{J} f_n(\mathbf{J}, t) \psi_n^{(p)*}(\mathbf{J}), \end{aligned} \quad (A11) $$

where the last expression was obtained thanks to Equation (A1a).

While for now, we only used the definition of the projection coefficient and some field equations, let us include the dynamics through the CBE, and in particular Equation (A4), to get

$$ a_p(t) = -i(2\pi)^3 \int_0^t d\tau \sum_n \int d\mathbf{J} e^{-i\mathbf{n}\cdot\Omega(t-\tau)} \mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}} \times \psi_n^{(p)*}(\mathbf{J}) \Delta H_n(\mathbf{J}, \tau). \quad (A12) $$

To proceed further, we develop the perturbation to the Hamiltonian in more detail. Starting from Equation (2), let us first expand all potentials ($\psi^e$ and $\psi^s$) and densities ($\rho_1$ in

Equation (3)) in the potential-density basis elements. This gives

$$ \Delta H(\mathbf{x}, t) = \sum_q (a_q + b_q) \left[ \psi_n^{(q)}(\mathbf{x}) + \mathbf{x} \cdot \int \frac{d\mathbf{x}}{|\mathbf{x}|^2} \rho^{(q)}(\mathbf{x}) \mathbf{e}_r \right], $$

where the time dependence is fully borne by the projection coefficients $a_q$ and $b_q$. Defining the new set of functions of Equation (11), the angular Fourier transform of the perturbing Hamiltonian is given by

$$ \Delta H_n(\mathbf{J}, t) = \sum_q (a_q + b_q) [\psi_n^{(q)}(\mathbf{J}) + \phi_n^{(q)}(\mathbf{J})]. \quad (A13) $$

Using Equations (A12) and (A13), we get

$$ a_p(t) = \int_0^t d\tau \sum_q M_{pq}(\tau - \tau_q) [a_q(\tau) + b_q(\tau)], \quad (A14) $$

where the response matrix is defined as in Equation (10). Equation (A14) is also another version of Equation (9).

## Appendix B
### The Clutton-Brock Bi-orthogonal Basis

In this appendix, we detail our choice of basis functions, which was first constructed by Clutton-Brock (1973). The potential and density elements of this basis are given by Equation (A5a) with

$$ U_n^\ell(r) = A_n^\ell \frac{(r/R_b)^\ell}{(1 + (r/R_b)^2)^{\ell+1/2}} C_n^{(\ell+1)}(\chi) \quad (B1a) $$

$$ D_n^\ell(r) = B_n^\ell \frac{(r/R_b)^\ell}{(1 + (r/R_b)^2)^{\ell+5/2}} C_n^{(\ell+1)}(\chi), \quad (B1b) $$

where $R_b$ is the basis scale radius, $C_n^{(\alpha)}$ are the Gegenbauer polynomials, and the renormalized radius is given by

$$ \chi = \frac{(r/R_b)^2 - 1}{(r/R_b)^2 + 1}. \quad (B2) $$

In Equation (B1a), we defined the normalization constants

$$ A_n^\ell = - \sqrt{\frac{G}{K_n^\ell R_b}} 2^{2\ell+3} \ell! \sqrt{\frac{(n+\ell+1)n!}{(n+2\ell+1)!}}, \quad (B3a) $$

$$ B_n^\ell = - \frac{K_n^\ell}{4\pi G R_b^2} A_n^\ell, \quad (B3b) $$

where $K_n^\ell$ is defined as

$$ K_n^\ell = 4n(n + 2\ell + 2) + (2\ell + 1)(2\ell + 3). \quad (B4) $$

The radial basis is therefore defined for $0 \le n \le n_{max}$.

## Appendix C
### The Baes-van Hese Equilibrium Distribution Function

For the choice of anisotropic phase-space DFs with Hernquist density, we relied on the work of Baes & van Hese 2007 (their Equations (92) and (93)). More specifically, we focused on the particular case of spheres with a constant
anisotropy parameter $\beta$, so that the DF is given by

$$F(E, L) = \frac{1}{(2\pi)^{5/2}} \frac{1}{(GMa)^{3/2}} \left(-\frac{aE}{GM}\right)^{5/2-2\beta} \times \frac{\Gamma(5-2\beta)}{\Gamma(1-\beta)\Gamma(\frac{7}{2}-\beta)} \left(-\frac{L^2}{2a^2E}\right)^{-\beta} \times_2 F_1\left(5-2\beta, 1-2\beta, \frac{7}{2}-\beta, -\frac{aE}{GM}\right), \quad (C1)$$

where $a$ stands for the scale radius of the MW, $a_{\text{MW}}$, and $M$ for the total mass of the MW, $M_{\text{MW}}$. In order to produce a non-negative DF, the anisotropy parameter is restricted to $\beta \le 0.5$ (see An & Evans 2006). Note that, in the case where the DF only represents the stellar halo (Sections 4 and 5), this DF should merely be rescaled by the factor $M_{\text{tot}}/M_{\text{MW}}$, without rescaling the energy or the angular momentum (but see Section 4.2).

An extra step is required to consider this DF as a function of the actions, so that it can be input in Equation (10). Indeed, once the spherically symmetric potential is specified, a bijective relation exists between the sets of conserved quantities, which are $(E, L)$ and $(J_r, L)$. In practice, all quantities that are entered into Equation (10) are actually computed on a third set of conserved quantities, the peri- and apocenters $(r_p, r_a)$, which make the coordinate transforms more straightforward. Some of the technical details of these coordinate transforms are given in Appendix D.

## Appendix D
### Computation of the Response Matrix

For a spherical, nonrotating mean-field DF, the formula giving the matrix method can be simplified. We perform these simplifications in the following section. Later, we describe in some detail the numerical techniques we developed for the computation of the matrix, and validate the implementation by recovering unstable modes from the literature.

#### D.1. Matrix of a Spherical, Nonrotating Equilibrium

Here, we consider the special case where $\Omega_3 = 0$ and $\partial F/\partial L_z = 0$, i.e., that of a spherical, nonrotating system. Let us first use the derivation from Tremaine & Weinberg (1984) for the Fourier-transformed basis functions,

$$\psi_n^{(p)}(\mathbf{J}) = \delta_{mn}^3 i^{m_p-n_2} Y_{\ell_p \ell}^{n_2}(\frac{\pi}{2}, 0) R_{n_2 m}^{\ell p}(\beta) W_{\ell' n' \ell}^{\tilde{n}}(\tilde{\mathbf{J}}), \quad (D1)$$

where $n = (n_1, n_2, n_3)$ is the resonance vector associated with each Fourier coefficient, $\beta$ is the inclination angle of the orbit associated with $\mathbf{J}$, defined so that $\cos(\beta) = L_z/L$, and $\tilde{\mathbf{J}} = (J_r, L)$, $\tilde{n} = (n_1, n_2)$. Additionally, the rotation matrix $R_{nm}^{\ell}(\beta)$ is defined as

$$R_{nm}^{\ell}(\beta) = \sum_{t} (-1)^{t} \frac{\sqrt{(\ell + n)! (\ell - n)! (\ell + m)! (\ell - m)!}}{(\ell - m - t)! (\ell + n - t)! t! (t + m - n)!} \times [\cos(\beta/2)]^{2\ell + n - m - 2t} [\sin(\beta/2)]^{2t + m - n}, \quad (D2)$$

where the sum over $t$ is restricted to the values such that the arguments of the factorials are positive, i.e., $t_{\min} \le t \le t_{\max}$, with $t_{\min} = \text{Max}[0, n - m]$ and $t_{\max} = \text{Min}[\ell - m, \ell + n]$. In Equation (D1), the Fourier-transformed in-plane radial

coefficients $W_{\ell n}^{\tilde{n}}(\tilde{\mathbf{J}})$ are defined as

$$W_{\ell n}^{\tilde{n}}(\tilde{\mathbf{J}}) = \frac{1}{\pi} \int d\theta_1 U_n^{\ell}(r(\theta_1)) \cos[n_1\theta_1 + n_2(\theta_2 - \xi)], \quad (D3)$$

which are real for real radial basis functions. In this integral, the radial dependence of the angles $\theta_1$ and $(\theta_2 - \xi)$ is given by

$$\theta_1 = \Omega_1 \int_C dr \frac{1}{\sqrt{2(E - \psi_0(r)) - L^2/r^2}}, \quad (D4a)$$

$$\theta_2 - \xi = \int_C dr \frac{\Omega_2 - L/r^2}{\sqrt{2(E - \psi_0(r)) - L^2/r^2}}, \quad (D4b)$$

where $\xi$ is the angle between the ascending node and the current position, measured in the orbit plane along the orbital motion, $E, L$ are the energy and angular momentum of the orbit, and $C$ is the integration contour going from the pericenter $r_p$ up to the current position $r = r(\theta_1)$ along the radial oscillation.

In order to simplify Equation (10), we also need to decompose the inertial term $\phi_n^{(p)}(\mathbf{J})$. Using Equation (11), we have

$$\phi_n^{(p)}(\mathbf{J}) = G \left( \int \frac{d\theta}{(2\pi)^3} \mathbf{x} e^{-i n \cdot \theta} \right) \cdot \left( \int \frac{d\mathbf{x}}{|\mathbf{x}|^2} \rho^{(p)}(\mathbf{x}) e_r \right). \quad (D5)$$

Let us first focus on the integral on the left-hand side of the dot product, which we will rewrite as $\mathbf{x}_n(\mathbf{J})$, as it is the angular Fourier transform of the position vector. We can rewrite $\mathbf{x}$ as a Cartesian vector in terms of spherical coordinates $(r, \theta, \phi)$ as

$$\mathbf{x} = r \mathbf{e}_r = r \sqrt{\frac{2\pi}{3}} \begin{pmatrix} Y_1^{-1}(\theta, \phi) - Y_1^1(\theta, \phi) \\ i[Y_1^{-1}(\theta, \phi) + Y_1^1(\theta, \phi)] \\ \sqrt{2} Y_1^0(\theta, \phi) \end{pmatrix}. \quad (D6)$$

In each of the Cartesian directions, we therefore have to perform an angular Fourier transform of a function that is separable in terms of a linear combination of spherical harmonics $Y$, times a function that depends on the radius only (in the present case, $r$ itself). This is precisely the context in which Equation (D1) can be applied. As a result, we have

$$\mathbf{x}_n(\mathbf{J}) = \sqrt{\frac{2\pi}{3}} Y_1^{n_2}(\frac{\pi}{2}, 0) X^{\tilde{n}}(\tilde{\mathbf{J}}) i^{-n_2} \times \begin{pmatrix} -i[\delta_{-1}^{n_3} R_{n_2-1}^{1}(\beta) + \delta_1^{n_3} R_{n_2-1}^{1}(\beta)] \\ \delta_{-1}^{n_3} R_{n_2-1}^{1}(\beta) - \delta_1^{n_3} R_{n_2-1}^{1}(\beta) \\ \sqrt{2} \delta_0^{n_3} R_{n_2-0}^{1}(\beta) \end{pmatrix}, \quad (D7)$$

where we defined new in-plane radial coefficients as

$$X^{\tilde{n}}(\tilde{\mathbf{J}}) = \frac{1}{\pi} \int d\theta_1 r(\theta_1) \cos[n_1\theta_1 + n_2(\theta_2 - \xi)]. \quad (D8)$$

Similarly, we can express the integral on the rhs of the dot in Equation (D5) as a Cartesian vector in terms of spherical harmonics. Given the harmonic dependence of the basis functions (see Equation (A5a)), as well as the orthogonality of the spherical harmonics, most of the basis functions will
yield a vanishing integral. The only nonzero terms give

$$ \int \frac{d\mathbf{x}}{|\mathbf{x}|^2} \rho^{(p)}(\mathbf{x}) e_r = \sqrt{\frac{2\pi}{3}} \delta^{1}_{\ell p} \begin{pmatrix} \delta_{m^p}^{-1} - \delta_{m^p}^{1} \\ -i(\delta_{m^p}^{-1} + \delta_{m^p}^{1}) \\ \sqrt{2} \delta_{m^p}^{0} \end{pmatrix} d_{n^p}, \quad (\text{D9}) $$

where we defined the radial integral of the basis functions over the whole radial range as

$$ d_n = \int_0^\infty dr \, D_n^1(r). \quad (\text{D10}) $$

If we now perform the scalar product of Equations (D7) and (D9), which are both written in Cartesian coordinates, we have

$$ \phi_n^{(p)}(\mathbf{J}) = \frac{4\pi G d_{n^p}}{3} \delta_{\ell p}^{1} \delta_{m^p}^{n_3} i^{m^p - n_2} Y_1^{n_2}\left(\frac{\pi}{2}, 0\right) \times R_{n_2 m^p}^{1}(\beta) X^n(\mathbf{J}). \quad (\text{D11}) $$

The fact that this term accounts for the translation of the reference frame is recovered, since it is only present in dipolar harmonics $\ell^p = 1$.

We can now use Equations (D1) and (D11) to simplify Equation (10). In Equation (D11), the dependence on the third action, $L_z$, is only borne by the rotation matrices $R(\beta)$. We can therefore make use of their orthogonality relation,

$$ \int_0^\pi d\beta \sin(\beta) R_{n_2 m^p}^{\ell p}(\beta) R_{n_2 m^p}^{q}(\beta) = \delta_{\ell p}^{q} \frac{2}{2\ell^p + 1}. \quad (\text{D12}) $$

Once this simplification is performed, we end up with the final form of the response matrix for spherical, nonrotating systems,

$$ M_{pq}(t) = \delta_{m^p}^{m_q} \delta_{\ell p}^{\ell q} \sum_{\tilde{n}} C_{\ell \tilde{n}}^{n_2} P_{\ell \tilde{n} p n_q}^{\tilde{n}}, \quad (\text{D13}) $$

where we define the coefficients $C_\ell^n$ as

$$ C_\ell^n = -2i(2\pi)^3 \frac{\left|Y_\ell^n\left(\frac{\pi}{2}, 0\right)\right|^2}{2\ell + 1}, \quad (\text{D14}) $$

and the functions $P_{\ell \tilde{n} p n_q}^{\tilde{n}}(t)$ as

$$ P_{\ell \tilde{n} p n_q}^{\tilde{n}}(t) = \int d\tilde{\mathbf{J}} L \tilde{\mathbf{n}} \cdot \frac{\partial F}{\partial \tilde{\mathbf{J}}} e^{-i \tilde{n} \cdot \tilde{\Omega} t} W_{\ell \tilde{n} p n_q}^{\tilde{n}}(\tilde{\mathbf{J}}) \times \left[ W_{\ell \tilde{n} p n_q}^{\tilde{n}}(\tilde{\mathbf{J}}) + \delta_{\ell p}^{1} \frac{4\pi G}{3} d_{n^q} X^{\tilde{n}}(\tilde{\mathbf{J}}) \right]. \quad (\text{D15}) $$

One can notice the similarity of this equation to Equation (23) in Murali (1999), the main difference being our definition of $d_n^q$ (their $p_j^{lm}$), which stems from considering the reference frame of the cusp instead of that of the COM. Interestingly enough, the response matrix element $M_{pq}$ is proportional to $\delta_{m^p}^{m_q} \delta_{\ell p}^{\ell q}$. This means that there is no coupling between different angular harmonics in the system's response: each angular harmonic effect in the response is only induced by the corresponding harmonic cause in the perturber, mediated by the corresponding harmonic term in the matrix. As shown by Rozier et al. (2019), this characteristic is specific to nonrotating spheres. Let us now detail the numerical methods that we used to compute the response matrix.

## D.2. Numerical Methods

To compute the action space integral of Equation (D15), we carefully analyzed the different terms of the integrand. A critical feature appears when this integrand is rewritten as $g(\tilde{\mathbf{J}}) e^{i h(\tilde{\mathbf{J}})}$. In this form, the integrand is written as a slowly varying function of the actions, $g$, times a fast trigonometric oscillation. The argument of this oscillating term itself, $h$, is also a slowly varying function of the actions. We therefore choose to divide the action space into small surfaces, in which both functions $g$ and $h$ are well approximated by their first-order expansion.

In order to reach a better sampling of action space, we relied on the same change of variables as in Rozier et al. (2019): the integration variables are changed to $(u, v)$, which are written as functions of the orbits' peri- and apocenters. This change of variables allows for a logarithmic sampling of the orbits that are either close to the host's center or close to circular, while the other orbits are sampled linearly in terms of peri- and apocenter. This usually leads to better sampling of the regions where the integrand of Equation (D15) reaches a significant amplitude. The function $g$ can be redefined to include the transformation's Jacobian, and the resulting functions $g(u, v)$ and $h(u, v)$ are still slowly varying functions of their arguments, as compared to the fast trigonometric oscillation.

In the end, we compute Equation (D15) as the sum over a grid on the $(u, v)$ surface of the integral $\Re(g, \frac{\partial g}{\partial u}, \frac{\partial g}{\partial v}, h, \frac{\partial h}{\partial u}, \frac{\partial h}{\partial v})$ defined by

$$ \Re(a, b, c, d, e, f) = \int \int_{-\frac{\Delta u}{2}}^{\frac{\Delta u}{2}} du dv \, (a + bu + cv) e^{i(d+eu+fv)}, \quad (\text{D16}) $$

where $g$, $h$, and their derivatives are evaluated at the center $(u_0, v_0)$ of each square of side $\Delta u$. This integral can be renormalized as

$$ \Re(a, b, c, d, e, f) = \Delta u^2 a \, e^{id} \Re_D \left( \frac{b \Delta u}{a}, \frac{c \Delta u}{a}, \frac{e \Delta u}{d}, \frac{f \Delta u}{d} \right), \quad (\text{D17}) $$

where the normalized integral is defined as

$$ \Re_D(b, c, e, f) = \int \int_{-\frac{1}{2}}^{\frac{1}{2}} dx dy \, (1 + bx + cy) e^{i(ex+fy)}. \quad (\text{D18}) $$

Finally, we found an analytical expression for this last integral as

$$ \Re_D(b, c, e, f) = \begin{aligned} & \mathrm{sinc}\left(\frac{e}{2}\right) \mathrm{sinc}\left(\frac{f}{2}\right) \\ & -i \frac{b}{e} \mathrm{sinc}\left(\frac{f}{2}\right) \left(\cos\left(\frac{e}{2}\right) - \mathrm{sinc}\left(\frac{e}{2}\right)\right) \\ & -i \frac{c}{f} \mathrm{sinc}\left(\frac{e}{2}\right) \left(\cos\left(\frac{f}{2}\right) - \mathrm{sinc}\left(\frac{f}{2}\right)\right), \end{aligned} \quad (\text{D19}) $$

where $\mathrm{sinc}(x) = \sin(x)/x$.

Let us now explain in more detail how we evaluate the functions $g$ and $h$, and their partial derivatives. In general, most of the functions involved in $g$ and $h$ can be expressed as functions of $(r_p, r_a)$, the orbit's peri- and apocenter. Such functions can later be considered as functions of $(u, v)$, owing
to the analytical relations $r_p(u)$ and $r_a(u, v)$ (see Rozier et al. 2019). In particular, the energy (required in the phase-space DF $F(E, L)$) and angular momentum are given by

$$E = \frac{r_a^2 \psi_0(r_a) - r_p^2 \psi_0(r_p)}{r_a^2 - r_p^2}, \qquad (D20a)$$

$$L = \sqrt{\frac{2(\psi_0(r_a) - \psi_0(r_p))}{r_p^2 - r_a^2}}, \qquad (D20b)$$

and the orbital frequencies are given by

$$\Omega_1 = \left[ \frac{1}{\pi} \int_{r_p}^{r_a} dr \frac{1}{\sqrt{2(E - \psi_0(r)) - L^2/r^2}} \right]^{-1}, \qquad (D21a)$$

$$\Omega_2 = \frac{\Omega_1}{\pi} \int_{r_p}^{r_a} dr \frac{L/r^2}{\sqrt{2(E - \psi_0(r)) - L^2/r^2}}. \qquad (D21b)$$

Owing to these relations, the functions $g$ and $h$ can be computed, as well as their partial derivatives through explicit analytical expressions. In the case of both $g$ and its partial derivatives, special treatment should be mentioned in the computation of $W$, $X$, and their partial derivatives (see Equations D3 and D8). Indeed, these functions a priori involve nested integrals of the form

$$\int dr S[r, \theta_1(r), (\theta_2 - \xi)(r)], \qquad (D22)$$

where $\theta_1(r)$ and $(\theta_2 - \xi)(r)$ themselves are integrals (see Equation (D4a)). Besides, the integrands involved are unbound at the edges of the integration region, which could be the source of issues when performing derivatives. To cure these two problems, we first regularize the integrals at their edges by applying the same effective anomalies as in Henon 1971 (see also Rozier et al. 2019). Then, following Rozier et al. (2019), instead of directly computing the nested integrals, we transform the problem into the single integration of a multicomponent vector. These tricks allow us to compute $W$, $X$, as well as their partial derivatives, as simple well-posed integrals using an RK4 integration scheme.

### D.3. Validation of the Implementation

In order to validate our implementation of the response matrix, here, we present the response of a radially anisotropic isochrone sphere. In Saha (1991), it is shown that such a sphere undergoes a radial orbit instability. We therefore compute the response matrix for a spherical system with an isochrone potential (Binney & Tremaine 2008), and an Osipkov–Merritt DF (Osipkov 1979; Merritt 1985) with a transition radius from the isotropic center to the radially anisotropic outskirts taken at $R_a = 1.0b$, where $b$ is the scale radius of the isochrone potential. Since the instability is expected to emerge as a quadrupolar ($l=2$) mode, we focus on this component of the response matrix and we perturb the system with a potential presenting the corresponding level of symmetry. The computation uses 100 basis functions as defined in Appendix B, with a scale radius $R_b = 20b$, and a maximum radial resonance number of $n_{l \max} = 5$. In this section, all quantities are rescaled so that $G = M_{\text{tot}} = b = 1$, where $M_{\text{tot}}$ is the total mass of the sphere.

<figure><img src="image_16.png" /><figcaption>Figure D1. Time evolution of the norm of the coefficients of the perturber $|b_n|$ (dashed lines) and the response $|a_n|$ (solid lines). Different colors indicate different values of $n$. Note that the different $|b_n|$ values were artificially shifted in order to avoid overlapping. Despite the early disappearance of the external perturber, the system spontaneously develops an exponentially growing mode at the rate $\eta = 0.0245$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure D2. Radial profile of the instability's potential, as computed by Saha (1991) (dashed line) and by our method (solid line). The vertical scale of each curve is arbitrary because the instability grows exponentially in time, starting from a seed that is specific to the particular realization of the initial conditions. Therefore, we decided to match the maximum values of both curves. The shapes are very similar, which further validates our implementation.</figcaption></figure>

By definition, the unstable mode is a property of the background sphere; as a consequence, it is independent of the perturber we apply. We therefore arbitrarily choose an axisymmetric perturber ($m=0$), and instead of choosing a potential for the perturber that would later be projected onto the functional basis to get the vector $b$, we directly fix the value of the vector as well as its time evolution. For simplicity, we choose to give equal values to all vector coefficients, with a rapid cutoff in time of the form $e^{-t^2/40}$. The time evolution of these coefficients is shown in Figure D1.

Figure D1 also shows the time evolution of the coefficients of the response vector $a$, when the response matrix formalism is applied according to Equation (15) and evolved up to $t=300$. Although the perturber is steeply cutoff after $t \sim 10$, it has clearly excited an instability that later grows at an exponential rate: indeed, above $t \sim 25$, all coefficients of $a$ grow at the same exponential pace. When we compute the common slope of these curves, we can estimate the growth rate of the identified instability to be $\eta = 0.0245$, which favorably compares to the computation from Saha (1991) at $\eta = 0.025$. In addition to the
norm of these coefficients, we measured their phase angle, and did not detect any variation of it. This indicates that the mode is not oscillating, which is also the conclusion of Saha (1991).

Focusing on the spatial shape of the instability, Figure D2 compares the radial profile of the mode's potential from our computation at $t=300$ to the same measurement from Saha (1991) (both using an arbitrary normalization of the potential's amplitude). The two profiles look very much alike, despite the use of different computation parameters as well as different detection methods (time space in our case versus frequency space in theirs).

Overall, this comparison suggests that our algorithm is valid for our purposes.

## Appendix E

### Optimized Operations with the Response Matrix

Equation (15) involves the inversion of a large matrix, as well as its product with a large vector. In order to perform these operations, we developed fast algorithms that take into account their specificities.

As defined by Equation (13), the matrix $M$ is a block-triangular-Toeplitz matrix. Such a matrix is defined by the following particular shape:

$$ \begin{pmatrix} M_0 & & & & \\ M_1 & \mathbf{0} & & & \\ & \ddots & \ddots & & \\ & & M_2 & M_1 & M_0 \\ & & & M_K & \end{pmatrix}, \qquad (E1) $$

where the diagonal blocks are denoted as $M_i$. One can easily show that the inverse of such a matrix is also a block-triangular-Toeplitz matrix. This inverse can be computed by the following recurrence. Let us first coin $A_i$ the blocks of the inverse matrix. The first block is straightforwardly given by $A_0 = M_0^{-1}$, where this inversion is computed using a standard matrix inversion scheme. Then, for $1 \le i \le K$, the matrix $A_i$ is computed thanks to the relation

$$ A_i = -M_0^{-1} \sum_{k=0}^{i-1} M_{k+1} A_{i-1-k}. \qquad (E2) $$

One can easily show that such an algorithm indeed yields the inverse of the original matrix. Note that, when inverting $I - M$, the inversion of the first term is straightforward, since it is equal to the identity itself.

A similar algorithm can be developed for the product of such a block-triangular-Toeplitz matrix with a vector. When the matrix (E1) is multiplied with a vector defined by the stacked sub-vectors $b_0, \dots, b_K$, then the stacked sub-vectors of the product, which we denote as $c_0, \dots, c_K$ can be computed via

$$ c_i = \sum_{k=0}^{i} M_k b_{i-k}. \qquad (E3) $$

### ORCID iDs

Simon Rozier https://orcid.org/0000-0001-6994-6708

Benoit Famaey https://orcid.org/0000-0003-3180-9825

Rodrigo Ibata https://orcid.org/0000-0002-3292-9709

### References

An, J. H., & Evans, N. W. 2006, ApJ, 642, 752

Baes, M., & van Hese, E. 2007, A&A, 471, 419

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

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

Bertin, G., Pegoraro, F., Rubini, F., & Vesperini, E. 1994, ApJ, 434, 94

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

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: 2nd ed. (Princeton, NJ: Princeton Univ. Press)

Bird, S. A., Xue, X.-X., Liu, C., et al. 2019, AJ, 157, 104

Breen, P. G., Rozier, S., Heggie, D. C., & Varri, A. L. 2021, MNRAS, 502, 4762

Choi, J.-H., Weinberg, M. D., & Katz, N. 2007, MNRAS, 381, 987

Clutton-Brock, M. 1973, Ap&SS, 23, 55

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

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

Danovich, M., Dekel, A., Hahn, O., Ceverino, D., & Primack, J. 2015, MNRAS, 449, 2087

De Rijcke, S., Fouvry, J.-B., & Pichon, C. 2019, MNRAS, 484, 3198

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

Fridman, A. M., Polyachenko, V. L., Aries, A. B., & Poliakoff, I. N. 1984, Physics of Gravitating Systems, Vol. II (Berlin: Springer-Verlag), 289

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

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

Garavito-Camargo, N., Patel, E., Besla, G., et al. 2021b, ApJ, 923, 140

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

Henon, M. 1971, Ap&SS, 14, 151

Henon, M. 1973, A&A, 24, 229

Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375

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

Kalnajs, A. J. 1977, ApJ, 212, 637

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

Lilley, E. J., Sanders, J. L., & Evans, N. W. 2018, MNRAS, 478, 1281

Meiron, Y., Li, B., Holley-Bockelmann, K., & Spurzem, R. 2014, ApJ, 792, 98

Merritt, D. 1985, AJ, 90, 1027

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Murali, C. 1999, ApJ, 519, 580

Osipkov, L. P. 1979, SvAL, 5, 42

Palmer, P. L., Papaloizou, J., & Allen, A. J. 1989, MNRAS, 238, 1281

Pawlowski, M. S., Oria, P.-A., Taibi, S., Famaey, B., & Ibata, R. 2022, ApJ, 932, 70

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

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

Petersen, M. S., Weinberg, M. D., & Katz, N. 2022, MNRAS, 510, 6201

Pichon, C., & Aubert, D. 2006, MNRAS, 368, 1657

Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493

Plummer, H. C. 1911, MNRAS, 71, 460

Polyachenko, V. L., & Shukhman, I. G. 1981, SvA, 25, 533

Rahmati, A., & Jalali, M. A. 2009, MNRAS, 393, 1459

Riley, A. H., Fattahi, A., Pace, A. B., et al. 2019, MNRAS, 486, 2679

Rozier, S., Fouvry, J. B., Breen, P. G., et al. 2019, MNRAS, 487, 711

Saha, P. 1991, MNRAS, 248, 494

Salem, M., Besla, G., Bryan, G., et al. 2015, ApJ, 815, 77

Sanders, J. L., Lilley, E. J., Vasiliev, E., Evans, N. W., & Erkal, D. 2020, MNRAS, 499, 4793

Schommer, R. A., Suntzeff, N. B., Olszewski, E. W., & Harris, H. C. 1992, AJ, 103, 447

Seguin, P., & Dupraz, C. 1994, A&A, 290, 709

Stewart, K. R., Kaufmann, T., Bullock, J. S., et al. 2011, ApJ, 738, 39

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

Tremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729

van der Marel, R. P., & Kallivayalil, N. 2014, ApJ, 781, 121

Vauterin, P., & Dejonghe, H. 1996, A&A, 313, 465

Weinberg, M. D. 1989, MNRAS, 239, 549

Weinberg, M. D. 1991, ApJ, 368, 66

Weinberg, M. D. 1998, MNRAS, 299, 499

Weinberg, M. D. 1999, AJ, 117, 629

Zang, T. A. 1976, PhD thesis, MIT

Zhao, H. 1996, MNRAS, 278, 488