# Collisionless relaxation from near equilibrium configurations: Linear theory and application to tidal stripping

Simon Rozier, Raphaël Errani

► To cite this version:

Simon Rozier, Raphaël Errani. Collisionless relaxation from near equilibrium configurations: Linear theory and application to tidal stripping. The Astrophysical Journal, 2024, 971 (1), pp.91. 10.3847/1538-4357/ad4c6e. hal-04480658

HAL Id: hal-04480658

https://hal.science/hal-04480658v1

Submitted on 17 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

# Collisionless Relaxation from Near-equilibrium Configurations: Linear Theory and Application to Tidal Stripping

Simon Rozier¹ and Raphaël Errani²,³

¹ School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Kings Buildings, Edinburgh, EH9 3FD, UK;
simon.rozier@ed.ac.uk

² McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

³ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France

*Received 2024 February 6; revised 2024 May 14; accepted 2024 May 14; published 2024 August 7*

## Abstract

Placed slightly out of dynamical equilibrium, an isolated stellar system quickly returns toward a steady virialized state. We study this process of collisionless relaxation using the matrix method of linear response theory. We show that the full phase-space distribution of the final virialized state can be recovered directly from the disequilibrium initial conditions, without the need to compute the time evolution of the system. This shortcut allows us to determine the final virialized configuration with minimal computational effort. Complementing this result, we develop tools to model the system's full time evolution in the linear approximation. In particular, we show that moments of the velocity distribution can be efficiently computed using a generalized moment matrix. We apply our linear methods to study the relaxation of energy-truncated Hernquist spheres, mimicking the tidal stripping of a cuspy dark matter subhalo. Comparison of our linear predictions against controlled, isolated N-body simulations shows agreement at percent level for the parts of the system where a linear response to the perturbation is expected. We find that relaxation generates a tangential velocity anisotropy in the intermediate regions, despite the initial disequilibrium state having isotropic kinematics. Our results also strengthen the case for relaxation depleting the amplitude of the density cusp, without affecting its asymptotic slope. Finally, we compare the linear theory against an N-body simulation of tidal stripping on a radial orbit, confirming that the theory still accurately predicts density and velocity dispersion profiles for most of the system.

*Unified Astronomy Thesaurus concepts: Cold dark matter (265); Tidal interaction (1699); Galaxy dynamics (591); Perturbation theory (1215); Analytical mathematics (38); N-body simulations (1083)*

## 1. Introduction

On the long term, stellar systems evolve through a succession of dynamical (quasi-) equilibria. This evolution can either be driven by internal processes (such as secular relaxation or instabilities), or be externally induced (like, for example, through tidal forces or the accretion of mass). Analytical methods that model the return of a stellar system toward equilibrium in response to a perturbation can be broadly categorized into two classes, depending on whether the perturbation is assumed to act on a timescale that is long with respect to the system's internal dynamics, or short.

For slow perturbations, Young (1980) developed a method based on the *adiabatic invariance* of the orbital actions to study the response of a star cluster to the (slow) mass growth of a central black hole. Sellwood & McGaugh (2005) applied the same formalism to the contraction of a dark matter (DM) halo due to a growing baryonic disk, and more recently Stücker et al. (2023a) used it to model the tidal evolution of DM subhalos. The Young (1980) method allows for a computationally less costly calculation of the new equilibrium configuration than an N-body simulation would. However, the underlying assumptions do limit the range of problems that the method can be applied to. First, the assumption of a slow perturbation suppresses any reference to time evolution in the computation of the final state, so that there is no way to

understand how the system reaches its final equilibrium. Second, the *adiabatic invariance* of the actions is not based on a solid theoretical background (Weinberg 1994a, 1994b), and may be challenged as soon as one departs from the purely spherically symmetric system, i.e., as soon as the *quasi-periodic* nature of the orbits comes into play.<fnref n="4" />

Methods developed to study the effect of a very fast perturbation with respect to the system's internal dynamics are generally based on the *impulse approximation*. These methods consider that the perturbation transfers an instantaneous velocity kick to each particle in the system (see, e.g., Kundic & Ostriker 1995; Gnedin et al. 1999), which then settles into a new equilibrium, following some analytical prescriptions (such as, for example, the conservation of the energy of individual shells of material in Dutton et al. 2016). Methods that make use of the impulse approximation have been used to estimate the effect of feedback-driven outflows on the central structure of DM halos (e.g., Pontzen & Governato 2012; Freundlich et al. 2020; Li et al. 2023), as well to model the tidal evolution of galaxies and DM subhalos (e.g., Gnedin & Ostriker 1999; Taylor & Babul 2001; Drakos et al. 2020; Benson & Du 2022). Not all parameters for these models follow directly from first principles, but some require calibration to simulation results (e.g., the adiabatic parameter of Gnedin & Ostriker 1999, or the diffusion parameter of Benson & Du 2022).

<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. 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="4">⁴ We note that in Stücker et al. (2023a), the method based on adiabatic conservation of the actions is much more accurate when the tide is spherically symmetric than when it is not.</footnote>
The application of either the impulse or the adiabatic approximation becomes challenging as soon as the perturbed system spans a wide range of dynamical timescales, causing some parts of the system to perceive a given perturbation as adiabatic, and others as impulsive (e.g., Weinberg 1994a, 1994b; Gnedin & Ostriker 1999),

In the present study, we develop an analytical method to model the process of collisionless relaxation for the specific case of weak perturbations to an initial equilibrium configuration. In this context, our method overcomes some of the limitations of existing tools, while still benefiting from the advantages of an analytical framework. Our model is rooted in the Kalnajs (1977) matrix method of linear response theory (LRT), which allows us to explicitly solve for the evolution of the phase-space distribution function (DF) at linear order in the perturbation, via the linearization of the collisionless Boltzmann-Poisson system of equations.

As an example, we apply our method to study the response of Hernquist (1990) spheres with initially isotropic kinematics to energy-based mass removal (similar to the initial conditions studied in Drakos et al. 2020, Amoros 2021, Errani et al. 2022, and Stücker et al. 2023a), approximating the type of perturbation typical for weak tidal encounters of cold DM subhalos. A theoretical understanding of the detailed phase-space structure of tidally stripped subhalos is of relevance, for example, for the modeling of Milky Way dwarf galaxies (see, e.g., Peñarrubia et al. 2008, 2010), the search for dark substructures through strong gravitational lensing (e.g., Vegetti & Koopmans 2009; Vegetti et al. 2010; Despali et al. 2022), as well as the modeling of potential DM annihilation or decay signals (e.g., Goerdt et al. 2007; Stref et al. 2019; Facchinetti et al. 2022; Stücker et al. 2023b; Delos & White 2023; Lovell 2024).

The detailed evolution of the inner regions of cold DM halos under the effect of tides remains a matter of active debate. Multiple studies that use N-body simulations to model the tidal stripping of cuspy subhalos strongly suggest a depletion in the amplitude of the central density cusp (e.g., Hayashi et al. 2003; Peñarrubia et al. 2008, 2010; Green & van den Bosch 2019; Errani & Navarro 2021), while recent studies based on analytical arguments suggest that the very center should remain virtually unaffected by tides (Drakos et al. 2020, 2022; Stücker et al. 2023a). Motivated by these conflicting predictions, we will address the tidal evolution of the central density cusp using LRT.

The paper is structured as follows. First, we show that we can cast the collisionless Boltzmann equation (CBE) so as to directly infer the characteristics of the final state from the initial conditions, in a single computational step. This allows us to compute the DF in the final equilibrium with minimal computational effort. The method is presented in Section 2, which also includes a comparison against the results of controlled N-body simulations. Second, in Section 3, we turn our attention to the time evolution of the system, and develop a method to efficiently compute the velocity moments of the evolving DF. Again, we compare the linear model predictions against controlled N-body simulations. Finally, in Section 4, we discuss the application of this model to the relaxation of cuspy DM subhalos subject to tidal stripping. We show that the model is particularly suitable to study the tidal evolution of the very inner regions of a DM cusp, and that collisionless relaxation drives the depletion in amplitude of the DM cusp during tidal

stripping. Finally, we summarize our main results and conclusions in Section 5.

## 2. The New Equilibrium in a Single Computation

### 2.1. Model for a Weakly Perturbed Equilibrium

We aim at analyzing a relatively broad class of physical problems, in which a self-gravitating system initially close to an equilibrium relaxes to reach a steady state. We specify here the assumptions that we make on the initial state out of equilibrium.

Our first assumption, that the system is initially close to an equilibrium, implies that there exists a decomposition of the initial phase-space DF, $F_i$, and the corresponding gravitational potential, $\psi_i$, into a sum:

$$ 
\begin{aligned}
F_i(\mathbf{x}, \mathbf{v}) &= F_0(\mathbf{x}, \mathbf{v}) + u(\mathbf{x}, \mathbf{v}), \\
\psi_i(\mathbf{x}) &= \psi_0(\mathbf{x}) + \psi_u(\mathbf{x}),
\end{aligned} 
\quad (1) $$

where $F_0$ is the DF of an equilibrium configuration with $\psi_0$ its gravitational potential, and $u$ represents a perturbation of the DF on top of this equilibrium, $|u| \ll F_0$, generating a potential $\psi_u$.<fnref n="5" /> We also introduce a small external perturbing potential $\psi_{ext}$, with $|\psi_{ext}| \ll |\psi_0|$, which is allowed to vary in time and is added on top of $\psi_i$. Note that the decomposition of $F_i$ is not unique, given that an equilibrium is usually surrounded by other equilibria in a continuous way. Here, we assume the potential $\psi_0$ to be integrable, hence we can define a set of angle-action variables $(\theta, J)$ and apply Jeans' theorem (Binney & Tremaine 2008), so that $F_0$ is only a function of $J$. We further assume that one of these decompositions lets us express $u$ too as a function of $J$ only. This assumption means that in one of these equilibria, the perturbation that relates the equilibrium to our disequilibrium DF is independent of the orbital phase. We can then define a new function $g = 1 + u/F_0$, so that

$$ F_i(J) = g(J)F_0(J). \quad (3) $$

Assuming that the external perturber $\psi_{ext}$ is asymptotically time stationary, the initial disequilibrium configuration $F_i$ is expected to reach a new steady state.<fnref n="6" /> In the following sections, we show that this final state can be efficiently computed using LRT.

These assumptions on the initial state still describe a relatively broad class of physical problems. We can consider the following scenario: a system in dynamical equilibrium undergoes a change of its potential, from an external source, which then remains stationary in time. In that case, $u = 0$, $\psi_u = 0$, and the external perturber is represented by $\psi_{ext}$. This scenario can describe a rapid infall of gas in a galaxy, or on the contrary a fast loss of mass due to stellar feedback in a young star cluster. Another scenario, which we will follow more closely in this paper, is that of a change in the equilibrium DF. In that case, $\psi_{ext} = 0$ and the change in the DF is represented by $u$. This model can represent a system that is still going through

$$ \int dx dv F_0(\mathbf{x}, \mathbf{v}) = M_{\text{tot}}. \quad (2) $$

<footnote n="5">⁵ By convention, we normalize the DF to the total mass of the equilibrium,</footnote>

<footnote n="6">⁶ Implicit here is the fact that the system does not support any neutral mode, which could theoretically drive it into infinite oscillations.</footnote>
the consequences of a past perturbation, once the perturber has disappeared. We will show in Section 4.1 that this model can describe the way a DM subhalo evolves after tidal stripping.

## 2.2. Linear Response Theory

### 2.2.1. A Boltzmann Equation for the Final State

In order to model the revirialization process from the initial state ($F_i, \psi_i$), we rely on the linearization of the Vlasov-Poisson system of equations. Here, we show that in our case of interest, computing the final equilibrium can be done very efficiently.

When the initial conditions of Equation (1) are evolved in time, the initial DF and potential turn into a new state given by

$$ 
\begin{aligned}
F_i(\mathbf{x}, \mathbf{v}, t) &= F_i(\mathbf{x}, \mathbf{v}) + f(\mathbf{x}, \mathbf{v}, t), \\
\psi_i(\mathbf{x}, t) &= \psi_i(\mathbf{x}) + \psi_{\text{ext}}(\mathbf{x}, t) + \psi^s(\mathbf{x}, t),
\end{aligned}
\quad (4) $$

where $f$ represents the evolution of the DF and $\psi^s$ the corresponding evolution of the potential. These responses in the DF and in the potential are related through

$$ \psi^s(\mathbf{x}, t) = -G \int d\mathbf{x}' d\mathbf{v}' \frac{f(\mathbf{x}', \mathbf{v}', t)}{|\mathbf{x}' - \mathbf{x}|}. \quad (5) $$

The dynamical evolution of the system is then described by the CBE, which in an inertial reference frame can be written as

$$ \frac{\partial(F_i + f)}{\partial t} + \frac{\partial(F_i + f)}{\partial \theta} \cdot \frac{\partial(H_0 + \psi^e + \psi^s)}{\partial \mathbf{J}} - \frac{\partial(F_i + f)}{\partial \mathbf{J}} \cdot \frac{\partial(H_0 + \psi^e + \psi^s)}{\partial \theta} = 0, \quad (6) $$

where $H_0 = \frac{v^2}{2} + \psi_0(r)$ is the Hamiltonian in the equilibrium potential $\psi_0$, and $\psi^e = \psi_u + \psi_{\text{ext}}$. We can simplify this equation, as we know that neither $H_0$ nor $F_i$ depend on the angles $\theta$ (since $H_0$ is conserved along the orbits, and because of Equation (3)), and that $F_i$ is static in time. We also apply here the linear approximation, considering that $|f| \ll F_i$ and $|\psi^e|, |\psi^s| \ll |\psi_0|$, and neglecting second-order terms. Noting $\Omega(\mathbf{J}) = \partial H_0 / \partial \mathbf{J}$ the vector of orbital frequencies, the revirialization is consequently described by the linearized CBE:

$$ \frac{\partial f}{\partial t} + \Omega(\mathbf{J}) \cdot \frac{\partial f}{\partial \theta} - \frac{\partial F_i}{\partial \mathbf{J}} \cdot \frac{\partial (\psi^e + \psi^s)}{\partial \theta} = 0. \quad (7) $$

This equation, together with the linearized Poisson equation $\nabla^2 \psi^s = 4\pi G \rho^s$, allow for the computation of the full revirialization process in the linear regime (see Section 3).

In this section, we are merely interested in the final state after revirialization, and not in the details of the time evolution. Knowing that this final state is time stationary, it can be found by considering the time asymptotic limit of Equation (7), when $\partial f_i / \partial t = 0$, the subscript “f” denoting the characteristics of the final equilibrium. Hence, we have

$$ \Omega(\mathbf{J}) \cdot \frac{\partial f_i}{\partial \theta} - \frac{\partial F_i}{\partial \mathbf{J}} \cdot \frac{\partial (\psi^e + \psi^s)}{\partial \theta} = 0. \quad (8) $$

It should be noted that $(\theta, \mathbf{J})$ are still the angle-action coordinates in the equilibrium $\psi_0$.

It is interesting to see that the final state can be computed, independently of the specific path that the system took to reach it. In particular, it does not depend on the way the perturbation $\psi_{\text{ext}}$ is turned on, but only on its time asymptotic value.

Consequently, in the linear approximation, the structure of the final equilibrium does not depend on whether the perturber is turned on impulsively or adiabatically. This may explain why collisionless relaxation in the linear limit can be studied within both the fast (Amorisco 2021) and the slow (Stücker et al. 2023a) approximations. Two remarks can be made about this finding. First, remaining in the linear limit implies that the perturbation $\psi_{\text{ext}}$ is always small, including at intermediate times, and that the system's response does not involve nonlinear resonant interaction, in particular interaction with neutral or weakly damped modes (Weinberg 1994c; Heggie et al. 2020). Second, the fact that the final state is independent from the actual time evolution is true only at the level of the global properties, i.e., the total DF and derived quantities. The fate of individual particles, however, may depend on how the perturber forces the system. Li et al. (2023) illustrate that behavior: They show that when the perturber is instantaneous, the action distribution after relaxation can be conserved even if the action of each individual particle is not.

### 2.2.2. Final Potential

To solve the CBE-Poisson system of equations, we use the matrix method, as devised by Kalnajs (1977), which we adapt here to the specificities of Equation (8). First, for any integer vector $\mathbf{n}$, we multiply Equation (8) by $e^{-i\mathbf{n}\cdot\theta}$ and integrate over the angles to get

$$ f_{\mathbf{n},\mathbf{n}}(\mathbf{J}) = \frac{\mathbf{n} \cdot \partial F_i / \partial \mathbf{J}}{\mathbf{n} \cdot \Omega} (\psi_{\mathbf{n},\mathbf{n}}^e + \psi_{\mathbf{n},\mathbf{n}}^s), \quad (9) $$

where we defined the angle Fourier transform of any function $h$ of phase space to be

$$ h_{\mathbf{n}}(\mathbf{J}) = \int \frac{d\theta}{(2\pi)^3} h(\theta, \mathbf{J}) e^{-i\mathbf{n}\cdot\theta}. \quad (10) $$

Then, the perturbing potentials and densities are projected on a biorthogonal basis with potential elements $\psi^{(p)}$ and density elements $\rho^{(p)}$ related through the Poisson equation, so that

$$ 
\begin{aligned}
\psi^s(\mathbf{x}) &= \sum_{p} a_p \psi^{(p)}(\mathbf{x}), \\
\psi^e(\mathbf{x}) &= \sum_{p} b_p \psi^{(p)}(\mathbf{x}),
\end{aligned}
\quad (11) $$

with the biorthogonality condition

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

Following these definitions, the projections can be computed as, e.g.,

$$ b_{\mathbf{f},\mathbf{p}} = -\int d\mathbf{x} \, \psi_{\mathbf{f}}^{e}(\mathbf{x}) \rho^{(p)}(\mathbf{x}) = -\int d\mathbf{x} \, \rho_{\mathbf{f}}^{e}(\mathbf{x}) \psi^{(p)}(\mathbf{x}). \quad (13) $$

When applied to the final response with coefficients $a_{\mathbf{f},\mathbf{p}}$, the projection gives

$$ 
\begin{aligned}
a_{\mathbf{f},\mathbf{p}} &= -\int d\mathbf{x} \, \rho_{\mathbf{f}}^{s}(\mathbf{x}) \psi^{(p)}(\mathbf{x}) \\
&= -\int d\mathbf{x} d\mathbf{v} \, f_{\mathbf{f}}(\mathbf{x}, \mathbf{v}) \psi^{(p)}(\mathbf{x}),
\end{aligned}
\quad (14) $$

where we used the fact that

$$ \rho_{\mathbf{f}}^{s}(\mathbf{x}) = \int d\mathbf{v} \, f_{\mathbf{f}}(\mathbf{x}, \mathbf{v}). \quad (15) $$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Top panel: phase-space distribution function (DF) of the ergodic Hernquist sphere $F_0$ (black, Equation (26)), and of three near-equilibrium models $F_i$ (models 1, 2, and 3 in blue, orange, and red, respectively; see Equations (3) and (30) using parameters as listed in Table 1), as a function of the orbital energy $E$ (defined in the original Hernquist potential). The DFs are expressed in units of $M_{\text{tot}}r_s^{-3}(-E_0)^{-3/2}$. Arrows indicate the respective energy truncation scale $E_{\text{tr}}$ (Equation (30)) for the three different models. Bottom panel: initial differential energy distributions $dN_i/dE$ (Equation (28)) corresponding to the cases shown above. In each panel, blue dotted curves show the energy-truncated (disequilibrium) DF and energy distribution measured from an N-body simulation of a Hernquist subhalo in a tidal field (see Section 4.1.2 for details).</figcaption></figure>

From Equation (10), and because any function of phase space is $2\pi$ periodic in the angles, the Fourier series of $f_f$ follows the convention

$$f_f(\theta, \mathbf{J}) = \sum_n f_{f,n}(\mathbf{J}) e^{in\theta}, \quad (16)$$

which we can inject into Equation (14). We can now canonically exchange integration variables from $dxdv$ to $d\theta d\mathbf{J}$, with a Jacobian equal to 1. Therefore, we have

$$\begin{aligned} a_{f,p} &= -\sum_n \int d\mathbf{J} f_{f,n}(\mathbf{J}) \left(\int d\theta \, e^{-in\theta} \psi^{(p)}(\theta, \mathbf{J})\right)^* \\ &= -(2\pi)^3 \sum_n \int d\mathbf{J} f_{f,n}(\mathbf{J}) \psi_n^{(p)*}(\mathbf{J}), \end{aligned} \quad (17)$$

where $\psi^{(p)}(\theta, \mathbf{J}) = \psi^{(p)}(\mathbf{x}(\theta, \mathbf{J}))$. We can now use Equation (9) to get

$$a_{f,p} = -(2\pi)^3 \sum_n \int d\mathbf{J} \frac{n \cdot \partial F_i / \partial \mathbf{J}}{n \cdot \Omega} \psi_n^{(p)*}(\mathbf{J}) (\psi_{f,n}^e(\mathbf{J}) + \psi_{f,n}^s(\mathbf{J})). \quad (18)$$

Finally, we can expand the perturbing potentials in the last parenthesis onto the biorthogonal basis following Equation (11),

yielding

$$\begin{aligned} a_{f,p} = &-(2\pi)^3 \sum_q \sum_n \int d\mathbf{J} \frac{n \cdot \partial F_i / \partial \mathbf{J}}{n \cdot \Omega} \psi_n^{(p)*}(\mathbf{J}) \psi_n^{(q)}(\mathbf{J}) \\ &\times (b_{f,q} + a_{f,q}). \end{aligned} \quad (19)$$

This equation is a matrix equation:

$$a_f = M_f(a_f + b_f), \quad (20)$$

where the matrix $M_f$ is defined by

$$M_{f,pq} = -(2\pi)^3 \sum_n \int d\mathbf{J} \frac{n \cdot \partial F_i / \partial \mathbf{J}}{n \cdot \Omega} \psi_n^{(p)*}(\mathbf{J}) \psi_n^{(q)}(\mathbf{J}). \quad (21)$$

We can also rearrange Equation (20) to show the linear relation between the perturber and the response:

$$a_f = ([\mathbf{I} - M_f]^{-1} - \mathbf{I})b_f. \quad (22)$$

This relation between the time asymptotic response and the perturber highlights again the fact that the final equilibrium that the system reaches is already encoded in the initial conditions, and is therefore independent from the details of the system's temporal evolution, in the linear approximation. In other words, a transient linear perturbation cannot have an ever-lasting effect on the system.

### 2.2.3. Final Phase Space Distribution Function

From a computation of the final potential, we can recover the full phase-space DF of the final equilibrium, $F_f = F_i + f_f$. A priori, it can be obtained by combining Equations (16) and (9) to get

$$F_f(\theta', \mathbf{J}') = F_i(\mathbf{J}') + \sum_n \frac{n \cdot \partial F_i / \partial \mathbf{J}}{n \cdot \Omega(\mathbf{J})} [\psi_{f,n}^e(\mathbf{J}') + \psi_{f,n}^s(\mathbf{J}')] e^{in\theta'}, \quad (23)$$

where $(\theta, \mathbf{J})$ are the angle-action coordinates computed in the potential $\psi_0$, of the same phase-space point of angle-action coordinates $\theta'$, $\mathbf{J}'$ when computed in the final potential $\psi_f = \psi_0 + \psi_f^e + \psi_f^s$. We additionally know that the final state is a new equilibrium, so that it should be independent of the phase angles $\theta'$. Because of the linear approximation, full phase mixing of the final state is only approximately reached by LRT. In order to improve the quality of our reconstruction of the final equilibrium, we force phase mixing in the final state by considering the final DF to be

$$\bar{F}_f(\mathbf{J}') = \int \frac{d\theta'}{(2\pi)^3} F_f(\theta', \mathbf{J}'). \quad (24)$$

This DF can be used, together with the final potential $\psi_f$, to compute velocity moments in the final state.

### 2.3. Comparison of the Final State with N-body Simulations

In order to validate our theory for collisionless relaxation close to an equilibrium, we compare the outcome of linear theory to N-body simulations, in a set of near-equilibrium systems which we define in the following section.

#### 2.3.1. Model for a Near-equilibrium System

We construct our initial DF $F_i$ in the following way. First, we choose the equilibrium DF–potential pair, $F_0 - \psi_0$; then, we
choose a function g as in Equation (3) to perturb the equilibrium in agreement with our model assumptions.

We consider an equilibrium following the Hernquist potential (Hernquist 1990), i.e., with

$$ \psi_0(r) = -\frac{GM_{\text{tot}}}{r + r_s}, \quad (25) $$

where $M_{\text{tot}}$ is the total mass of the model and $r_s$ its scale radius. Furthermore, we consider the equilibrium DF to be ergodic, i.e.,

$$ F_0(E) = \frac{M_{\text{tot}}}{8\sqrt{2}\pi^3 r_s^3 (-E_0)^{3/2}} \left[ \frac{3\arcsin\left(\sqrt{\frac{E}{E_0}}\right)}{\left(1 - \frac{E}{E_0}\right)^{5/2}} + \frac{\sqrt{\frac{E}{E_0}}\left(1 - 2\frac{E}{E_0}\right)\left(8\left(\frac{E}{E_0}\right)^2 - 8\frac{E}{E_0} - 3\right)}{\left(1 - \frac{E}{E_0}\right)^2} \right], \quad (26) $$

where the scale energy is given by

$$ E_0 = \psi_0(0) = -\frac{GM_{\text{tot}}}{r_s}. \quad (27) $$

This DF is represented as a black line in the top panel of Figure 1. The corresponding differential energy distribution $dN/dE$ can be obtained by multiplying Equation (26) by the density of states $G_0(E)$:

$$ dN/dE = F_0(E) G_0(E), \quad (28) $$

where

$$ G_0(E) = \frac{2\sqrt{2}\pi^2 r_s^{5/2} \sqrt{GM_{\text{tot}}}}{3\left(\frac{E}{E_0}\right)^{5/2}} \left[ 3\left(8\left(\frac{E}{E_0}\right)^2 - 4\frac{E}{E_0} + 1\right) \times \arccos\left(\sqrt{\frac{E}{E_0}} - \sqrt{\frac{E}{E_0} - \left(\frac{E}{E_0}\right)^2}\right) \left(4\frac{E}{E_0} - 1\right) \left(2\frac{E}{E_0} + 3\right) \right]. \quad (29) $$

This differential energy distribution is represented as a black line in the bottom panel of Figure 1.

To perturb this equilibrium, we consider a function g of the following form:<fnref n="7" />

$$ g(E) = \frac{1}{2} \operatorname{erfc} \left( \frac{E - E_{\text{tr}}}{A} \right), \quad (30) $$

where “erfc” is the complementary error function. When $E = E_0$, we typically have $g \simeq 1$, and typically $g \simeq 0$ when $E = 0$. This family has two free parameters, $E_{\text{tr}}$ and $A$, which respectively control the position in energies where the transition from 1 to 0 occurs, and the width of this transition region. This behavior is illustrated in Figure 1, where we show the DF and differential energy distribution of three cases with $E_{\text{tr}}/|E_0| = -0.176, -0.244$, and $-0.312$ (models 1, 2, and 3; see Table 1). For all

<table><caption>**Table 1**
Initial Conditions of the Energy-truncated Hernquist Models Considered in This Study</caption><thead><tr><th></th><th><i>E</i><sub>tr</sub>/|<i>E</i><sub>0</sub>|</th><th><i>M</i><sub>tr</sub>/<i>M</i><sub>tot</sub></th><th><i>M</i><sub>tr</sub>(&lt;<i>r</i><sub>s</sub>)/(0.25 <i>M</i><sub>tot</sub>)</th></tr></thead><tbody><tr><td>model 1</td><td>-0.176</td><td>0.56</td><td>0.98</td></tr><tr><td>model 2</td><td>-0.244</td><td>0.43</td><td>0.95</td></tr><tr><td>model 3</td><td>-0.312</td><td>0.33</td><td>0.90</td></tr></tbody></table>

**Notes.** We denote by $E_{\text{tr}}$ the truncation energy as defined in Equation (30) (with $A/|E_0| = 0.077$ for all three models), and by $M_{\text{tr}}$ and $M_{\text{tot}}$ the total masses of the truncated model and the original Hernquist model, respectively. The third column lists the mass enclosed within the Hernquist scale radius $r_s$ normalized by the enclosed mass in the original Hernquist model $M_{\text{H}}(<r_s) = 0.25 M_{\text{tot}}$.

models considered in this study, we use the same value of $A/|E_0| = 0.077$.

Physically, the equilibrium is perturbed in the following way: Nearly all orbits with energies below the transition region are conserved in the perturbed initial state, while hardly any orbit with energy above the transition region is kept. This is aimed at modeling a system that has lost its high-energy orbits instantaneously due to an external process, which we show relates to the tidal stripping of DM substructure in Section 4.1. We already see from Figure 1 that the linear approximation is not verified in parts of the phase space with high energy, a point which we discuss in Section 2.3.4.

### 2.3.2. Application of the Matrix Method

Once the initial disequilibrium DF $F_1$ is defined, we apply the matrix method in the following steps:

1. Construct the perturbing density as

$$ \rho^e(\mathbf{x}) = \rho_u(\mathbf{x}) = \int d\mathbf{v} (g(\mathbf{x}, \mathbf{v}) - 1) F_0(\mathbf{x}, \mathbf{v}). \quad (31) $$

2. Project it onto the basis following Equation (13), to get the vector $b_f$. For these projections, we used a separable spherical basis with spherical harmonics as the angular functions, and radial functions from Clutton-Brock (1973), with $n_{\text{max}} = 201$ (i.e., 202 basis functions) and characteristic radius $R_b = 3 r_s$.<fnref n="8" />

3. Independently compute the response matrix of Equation (21) from $F_i$. The action-space integral considers all orbits with $0.002 r_s < r_p < r_a < 20 r_s$. The sum over angular Fourier numbers is performed with $n_2 = n_3 = 0$ (spherical symmetry of the background and the perturber) and $|n_1| \le 10$. We confirmed the convergence of our results with respect to (w.r.t.) changes in these parameters.

4. Apply Equation (22) to get the vector $a_f$. Then, the response density and potential can be recovered by deprojection using Equation (11).

5. Recover the full phase-space distribution using Equation (24). Owing to spherical symmetry, the integral over $\theta'$ in that equation reduces to an integral over $\theta'_1$. After a change of variables $\theta'_1 \to r$, we compute the $(v_r, v_t)$ coordinates of the phase-space point of coordinates $(r, E', L')$. Then, the values of $(E, L)$ can be easily found, while that of $\theta'_1$ is performed via a 1D integral (see, e.g.,

<footnote n="7">⁷ The function $g(E)$ takes the role of the “filter function” in Errani et al. (2022); see their Figure 4 and Equation (9). We opt for a different parameterization which is better suited to approximate the shape of the tidal energy truncation typical for mildly perturbed Hernquist spheres; see Section 4.1.2.</footnote>

<footnote n="8">⁸ Because of spherical symmetry, there is no actual angular dependence, i.e., the multipolar expansion is restricted to the first term, $\ell = m = 0$.</footnote>
Tremaine & Weinberg 1984). To compute moments of the velocity distribution, we interpolate the DF $\bar{F}_f$ on a grid of ($E'$, $L'$), then we perform standard integrals of the interpolated function.

The bottleneck of this procedure is the computation of the response matrix. We use numerical methods straightforwardly adapted from Rozier et al. (2019) and summarized in Appendix A.2, which we complemented with a method to include edge terms in the action-space integration, as highlighted in Appendix B.<fnref n="9" /> In Appendix A.3, we discuss the numerical complexity of the calculation. We show the results of the whole procedure in Section 2.3.4.

### 2.3.3. Isolated Simulations Setup

We now briefly outline the numerical setup of the isolated N-body runs, which we use to compare the matrix results against. Our simulation setup is strongly motivated by the type of controlled simulations studied in Amorisco (2021).

N-body models. We generate an equilibrium N-body realization with $N = 10^7$ particles of a spherical Hernquist model with isotropic velocity dispersion, using the implementation of Errani & Peñarrubia (2020).<fnref n="10" /> In a subsequent step, we build three different disequilibrium N-body models by truncating the Hernquist realization in energy, following a tapered truncation as parameterized by Equation (30) with $A/|E_0| = 0.077$ and $E_{tr}/|E_0| = -0.176, -0.244$, and $-0.312$ (models 1, 2, and 3; see Table 1). These disequilibrium models are obtained by rejection sampling. They consist of $5.6 \times 10^6$, $4.3 \times 10^6$, and $3.3 \times 10^6$ N-body particles, respectively.

Particle-mesh code. We compute the dynamical evolution of our N-body models using the particle-mesh code SUPERBOX (Fellhauer et al. 2000). The code employs a high- and a medium-resolution cubic grid with 256³ cells each, centered on and comoving with the N-body model, with resolutions of $\Delta x_1 \approx 0.008$ $r_s$ and $\Delta x_2 \approx 0.08$ $r_s$, respectively. A third, static grid with a lower resolution of $\Delta x_3 \approx 0.8$ $r_s$ contains the entire simulation volume. The time integration makes use of a leapfrog scheme with fixed time step $\Delta t = T_s/1600$, where

$$T_s = 4\pi r_s^{3/2} (GM_{\text{tot}})^{-1/2} \qquad (32)$$

denotes the period of a circular orbit with radius $r_s$ in the original Hernquist potential. We follow the virialization of the disequilibrium N-body model for a total time of $15 T_s$.

### 2.3.4. Comparison of the Final State between Matrix and Isolated Simulations

Figure 2 compares the final relaxed states computed with the matrix method against our N-body models, after evolving the latter in isolation for a duration of $t = 15 T_s$. The N-body models are shown as filled circles, with error bars indicating the level of Poisson noise. Analytical initial values, as well as the results of calculations performed using the matrix method, are shown as dashed and solid curves, respectively.

Dashed curves in the top panel of Figure 2 show the initial density profiles of the three disequilibrium models in blue, orange, and red (for models 1, 2, and 3, respectively; see Table 1). Densities are normalized by the Hernquist density,

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top panel: density $\rho(r)$ in the initial disequilibrium and final relaxed states, normalized by the Hernquist density $\rho_H(r)$. Blue, orange, and red dashed curves correspond to the energy-truncated initial conditions of models 1, 2, and 3 with parameters as listed in Table 1. Results of the matrix method are shown as solid curves, whereas the corresponding N-body simulations (at $t = 15 T_s$) are shown as filled circles. Error bars indicate the expected N-body Poisson noise in each bin. Middle panel: as before, but showing the radial velocity dispersion $\sigma_r(r)$. Bottom panel: as before, but showing the anisotropy parameter $\beta = 1 - \sigma_r^2/(2\sigma_r^2)$.</figcaption></figure>

$\rho_H$. As expected, the truncation in energy introduces a depletion in density in the outer regions. In the inner regions ($r \lesssim r_s$), the initial disequilibrium profiles approach the original Hernquist density. The corresponding final equilibrium density profiles are shown as solid curves of lighter shade. In the final equilibrium state, the inner density is slightly depleted compared to the initial disequilibrium state (the behavior in the very center is discussed in Section 4.2), whereas the density in the outer regions slightly increases. At radii where the results of the matrix method are likely affected by resolution limits of the underlying basis ($r \lesssim r_{bas}$), or where nonlinear evolution is expected and the matrix method is not applicable ($r \gtrsim r_{lin}$), the matrix results are shown in gray. We choose $r_{bas} = 3 R_b/n_{max}$,

<footnote n="9">⁹ Note the similarity between Equation (21) and their Equation (9).</footnote>

<footnote n="10">¹⁰ Available at https://github.com/rerrani/nbopy.</footnote>
which is approximately the position of the second node of the last basis function (i.e., that with the highest radial resolution), and $r_{\text{lin}}$ is the position at which $\rho^e(r) = 0.5 \rho_{\text{H}}$ in each case, i.e., already far from a linear perturbation. Remarkably, the density in the final state as computed with the matrix method is in excellent agreement with the results of our N-body simulations at all radii, including regions where the perturbation is strongly nonlinear.

We now turn our attention to the kinematics of the systems in question. The central panel of Figure 2 shows the radial velocity dispersion $\sigma_r$ of our models, in units of the isotropic Hernquist radial velocity dispersion, $\sigma_{H,r}$. As before, the initial disequilibrium systems are shown as dashed curves in red, orange, and blue, while the final states are shown as solid curves in lighter shades of the same colors. The results of the matrix computation are in excellent agreement with our N-body models at those radii where the system is expected to respond linearly to the initial perturbation ($r \le r_{\text{lin}}$). At larger radii, the dispersions measured in the N-body models are slightly larger than suggested by the (linear) matrix calculation. From these curves, we clearly see that our case-dependent criterion for linearity of the perturbation is valid: The velocity dispersion computed from LRT starts to depart from the nonlinear N-body profile as soon as, in each case, the perturber's density represents a dominant fraction of the initial density.

Finally, we compare the anisotropy measured in the N-body models against the matrix calculations. We quantify the anisotropy through the usual anisotropy parameter $\beta = 1 - \sigma_t^2/(2\sigma_r^2)$, where subscripts t and r denote tangential and radial components, respectively. Both the matrix calculation and the N-body model show that the relaxed equilibrium system remains isotropic in the center, but develops a slight tangential anisotropy ($0 \ge \beta \ge -0.1$) at radii of $1 \lesssim r_s \lesssim 3$. At larger radii, the matrix method is not applicable; the N-body models suggest that $\beta$ rises again at radii beyond a few $r_s$.

To summarize, the results of the matrix calculation are in excellent agreement with those of controlled N-body simulations at radii where a linear response to the perturbation is expected ($r \le r_{\text{lin}}$). Remarkably, the density profiles as computed with the matrix methods remain consistent with the N-body models also at radii beyond $r_{\text{lin}}$.

### 3. Evolution Toward the Relaxed State

As derived in the previous section, the final equilibrium configuration may be computed directly from the initial disequilibrium distribution, without explicitly modeling the system's time evolution. Complementing this result, we now turn our attention to model in detail how the equilibrium configuration is reached, and use the matrix method of LRT to predict the full time evolution of the system in response to a perturbation.

#### 3.1. Evolution of the Density

The evolution of a perturbed equilibrium's density with LRT has already been developed in the literature (Seguin & Dupraz 1994; Weinberg 1998; Murali 1999; Pichon & Aubert 2006; Dootson & Magorrian 2022). We briefly summarize the main procedure here, following Rozier et al. (2022). We start with the linearized CBE-Poisson system of partial differential equations (see Equation (7) and the following paragraph). We

then multiply the CBE by $e^{-in \cdot \theta}$ and integrate over the angles, which gives us a differential equation for each set of Fourier numbers $n$. Each equation can be solved for $f_n(J, t)$, yielding

$$f_n(J, t) = i \mathbf{n} \cdot \frac{\partial F_i}{\partial \mathbf{J}} \int_0^t d\tau [\psi_n^e(J, \tau) + \psi_n^s(J, \tau)] e^{-i n \cdot \Omega(t-\tau)} \quad (33)$$

Independently, we can project the response density $\rho^s(\mathbf{x}, t)$ onto the same basis as defined in Equation (11), with projection coefficients $a_p(t)$. The exact same mathematical steps as Equations (14)–(17) can be applied to $a_p(t)$, replacing the final state with the evolving state. Here, the Fourier-transformed DF should be replaced with Equation (33), yielding a matrix equation with a remaining time integral:

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

where the response matrix is defined by

$$M_{pq}(t) = -i(2\pi)^3 \sum_n \int d\mathbf{J} \, e^{-i\mathbf{n} \cdot \Omega_t} \mathbf{n} \cdot \frac{\partial F_i}{\partial \mathbf{J}} \psi_n^{(p)*}(\mathbf{J}) \psi_n^{(q)}(\mathbf{J}). \quad (35)$$

Finally, we can fully linearize the problem by discretizing the time integration. To that end, we fix a time interval $[0, T]$ for computing the evolution, we evenly divide it into $K+1$ steps $0 = t_0 < \dots < t_K = T$, and define new vectors $\mathbf{a}$ and $\mathbf{b}$ obtained by stacking all vectors $\mathbf{a}(t_0), \dots, \mathbf{a}(t_K)$ and $\mathbf{b}(t_0), \dots, \mathbf{b}(t_K)$, respectively. With these definitions, Rozier et al. (2022) show that the response $\mathbf{a}$ is linearly related to the perturber $\mathbf{b}$ through

$$\mathbf{a} = ([\mathbf{I} - \mathbf{M}]^{-1} - \mathbf{I}) \mathbf{b}, \quad (36)$$

where $\mathbf{I}$ is the identity matrix of suitable size, and the matrix $\mathbf{M}$ is defined by blocks, so that the block in the line $i$ and column $j$ (with $1 \le i \le K+1$ and $1 \le j \le K+1$) is given by

$$\underline{\underline{\mathbf{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} \quad (37)$$

with $\Delta t = T/K$.

#### 3.2. Evolution of the Velocity Distribution

We are also interested in the evolution of some moments of the velocity distribution during revirialization, such as the mean radial velocity or the velocity dispersions. We note that we can recover the full phase-space DF of the evolved system, using a method from Dury et al. (2008). Here, we show that moments of the velocity distribution can be computed efficiently in the framework of the matrix method.

Let us consider the computation of a moment $\langle \mu(\mathbf{v}) \rangle(\mathbf{x}, t)$ of the velocity distribution at a given position $\mathbf{x}$ and a given time $t$ of the system's evolution. By definition, we have

$$\langle \mu(\mathbf{v}) \rangle(\mathbf{x}, t) = \frac{1}{\rho(\mathbf{x}, t)} \int d\mathbf{v} \, \mu(\mathbf{v}) F_1(\mathbf{x}, \mathbf{v}, t), \quad (38)$$

where $\rho = \rho_0 + \rho_u + \rho^s$ is the system's total density, excluding the external part of the perturbation, and $F_1$ is defined by Equation (4). The integral in the right-hand side can be decomposed into the initial contribution in $F_i$ and the evolutionary contribution in $f$. While the initial contribution can be straightforwardly computed from the initial DF, the evolutionary part requires knowledge of the full DF perturbation $f$, which is
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top panel: difference between the initial disequilibrium density profile $\rho(r, t = 0)$, and the evolving density profile at subsequent times $\rho(r, t > 0)$, normalized by the Hernquist density $\rho_H(r)$. The initial disequilibrium profile corresponds to a truncation in energy parameterized through $\{E_{tr}/E_0, A/[E_0]\} = \{-0.176, 0.077\}$ in Equation (30) (model 1; see Table 1). Profiles computed using the matrix method are drawn as solid curves, while results of isolated N-body simulation are shown as filled circles. Error bars indicate the level of Poisson noise expected in the simulation. Second panel: as before, but showing the difference between the radial velocity dispersion profile in the initial disequilibrium state $\sigma_r(r, t = 0)$, and the evolved profile $\sigma_r(r, t > 0)$. Third panel: anisotropy parameter $\beta = 1 - \sigma_t^2/(2\sigma_r^2)$. Bottom panel: mean radial velocity $\langle v_r \rangle$.</figcaption></figure>

not as direct for now. So, let us consider the term defined by

$$ \rho(\mathbf{x}, t) \langle \mu \rangle_f(\mathbf{x}, t) = \int d\mathbf{v} \, \mu(\mathbf{v}) f(\mathbf{x}, \mathbf{v}, t). \quad (39) $$

Since it is a function of $\mathbf{x}$ and $t$, it can be projected onto the biorthogonal basis of densities of Equation (11), with

projection coefficients $[\mu]_p(t)$ given by

$$ [\mu]_p(t) = -\int d\mathbf{x} \, \rho(\mathbf{x}, t) \langle \mu \rangle_f(\mathbf{x}, t) \psi^{(p)}(\mathbf{x}). \quad (40) $$

Replacing with Equation (39), we have

$$ [\mu]_p(t) = -\int d\mathbf{x} \, d\mathbf{v} \, \mu(\mathbf{v}) f(\mathbf{x}, \mathbf{v}, t) \psi^{(p)}(\mathbf{x}). \quad (41) $$

Again, we can follow the same steps as Equations (14)–(17), yielding

$$ [\mu]_p(t) = -(2\pi)^3 \sum_n \int d\mathbf{J} \, f_n(\mathbf{J}, t) [\mu\psi^{(p)}]_n^*(\mathbf{J}). \quad (42) $$

When we inject Equation (33) into Equation (42), we see that it can be transformed into a similar matrix equation as Equation (34), i.e.,

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

where $\mathbf{a}$ and $\mathbf{b}$ are the same response and perturber vectors as in Equation (34), $[\mu]$ is the vector made of $[\mu]_p$, and $\mathbf{M}[\mu]$ is the moment response matrix, defined by

$$ \mathbf{M}[\mu]_{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_i}{\partial \mathbf{J}} \times [\mu \psi^{(p)}]_n^* \mathbf{J} \, \psi_n^{(q)} \mathbf{J}. \quad (44) $$

Equation (43) shows that the moment $\mu$ derives from integrating the dynamics in the fully perturbed potential $\mathbf{a} + \mathbf{b}$, keeping track of the value of the moment along the evolution through the moment response matrix $\mathbf{M}[\mu]$.

Finally, we can apply the same time discretization as in Section 3.1, define the stacked vector $[\mu]$ in the same way as $\mathbf{a}$ and $\mathbf{b}$, and the matrix $\mathbf{M}[\mu]$ in the same way as Equation (37), to get

$$ [\mu] = \mathbf{M}[\mu](\mathbf{a} + \mathbf{b}). \quad (45) $$

This equation shows that the full time evolution of the function $\rho \langle \mu \rangle_f$, i.e., the density-weighted perturbation to the velocity moment, is linearly related to the total potential perturbation to the system, $\psi^e + \psi^s$, through the moment response matrix $\mathbf{M}[\mu]$.

### 3.3. Application of the Matrix Method

Again, let us briefly detail how the matrix method is applied:

1. Use the same perturbing density as Equation (31).

2. Project it onto the basis following Equation (13), and stack the resulting projection in order to build the vector $\mathbf{b}$. This stacking of the same vector accounts for the fact that, in our case of interest, the perturber is time stationary. For the projections, we used the same basis as in Section 2.3.2, but now only $n_{max} = 100$ with characteristic radius $R_b = 3 r_s$. Indeed, integrating the time evolution increases the computational cost of the algorithm. Hence, we had to decrease its spatial resolution.

3. Independently compute the response matrix of Equation (35) from $F_i$. The action-space integral considers all orbits with $0.002 r_s < r_p < r_a < 20 r_s$. The sum over angular Fourier numbers is performed with $n_2 = n_3 = 0$ (spherical symmetry of the background and
the perturber) and $|n_1| \le 2$. The dynamics is integrated for a total time of $10 T_s$, sampled with 600 time steps. We confirmed the convergence of our results w.r.t. changes in these parameters. Rozier et al. (2022) discuss why such a small number of time steps is sufficient in the framework of LRT.

4. Apply Equation (36) and deproject the response using Equation (11).

The bottleneck of this procedure is, again, the computation of the response matrix. We use the LiRGHaM code, which we complemented with a refined treatment of integral edges, as detailed in Appendix B.<fnref n="11" /> On the fly, we also compute the moment matrices for the mean radial velocity $M[v_r]$ and the mean squared velocities $M[v_r^2]$ and $M[v_i^2]$, using Equation (44) and the simplifications detailed in Appendix A. We show the results of the procedure in Section 3.4.

### 3.4. Comparison of the Time Evolution with Simulations

Having developed the framework necessary to model the time evolution of the perturbed system with the matrix method, we now apply this framework to the example disequilibrium system introduced in Section 2.3.1. As previously, we compare the results of the matrix calculation against controlled N-body models, as described in Section 2.3.3.

Figure 3 shows the initial ($t=0$) disequilibrium system of model 1 ($|E_{tr}|/|E_0|, A/|E_0| = \{-0.176, 0.077\}$ in Equation (30)) in dark blue, as well as in subsequently lighter shades of blue its state after evolving in isolation for a duration of $t = 0.3 T_s, 0.5 T_s, 3 T_s$, and $10 T_s$, where $T_s$ is defined in Equation (32). N-body snapshots are shown as filled circles, with error bars indicating the level of Poisson noise. The analytical properties in the initial state, as well as the results of the matrix calculation are shown as solid curves. As before, regions where the matrix method is limited by the resolution of the underlying basis ($r \le r_{bas}$), and regions where nonlinear effects may become important ($r \ge r_{lin}$) are shown in gray.

The top panel of Figure 3 shows the difference in density between the initial disequilibrium state $\rho(r, t=0)$ and the evolved density at subsequent times $\rho(r, t > 0)$, normalized by the Hernquist density $\rho_H(r)$. The time evolution as computed with the matrix method is in agreement with the N-body results with percent-level accuracy for most radii. The agreement between the N-body model and the matrix calculation further improves in the later stages of the evolution, i.e., as the final equilibrium state is approached. Note the depletion in density in the inner region as a consequence of the relaxation process, occurring very early on in the integration. This feature is further discussed in Section 4.2.

The second panel shows the difference between the radial velocity dispersion in the initial disequilibrium state $\sigma_r(r, t=0)$ and the evolved states $\sigma_r(r, t > 0)$, in units of the Hernquist dispersion $\sigma_H(r)$. The time evolution as predicted by the matrix method is in good agreement with the N-body models in the inner regions. In the outer regions ($r > 1.5 r_s$), while qualitatively similar, the dispersions computed with the matrix method and measured in the N-body models differ by several percent. Also for the case of the velocity dispersion, the agreement between N-body model and matrix calculation improves as the final equilibrium state is approached.

The third panel focuses on the time evolution of the anisotropy parameter $\beta = 1 - \sigma_r^2/(2\sigma_r^2)$. In the inner regions, the system remains isotropic as it returns to equilibrium. At radii of $1 \le r/r_s \le 3$, the system develops a mild tangential anisotropy, as seen already in the calculation of the final relaxed state. The matrix calculation and the N-body models give qualitatively similar results, but the differences between the results is most pronounced here in the calculation of the anisotropy parameter compared to all other properties studied.

Finally, the bottom panel of Figure 3 shows the time evolution of the average radial velocity $\langle v_r \rangle$, in units of the peak circular velocity of the original Hernquist model:

$$V_s = [GM_{tot}/(4 r_s)]^{1/2}. \qquad (46)$$

At early times, a large fraction of the system shows an outwards motion, as a reaction to the abrupt disappearance of a fraction of the particles. Progressively, the region where $\langle v_r \rangle \approx 0$ grows outwards as the system returns to equilibrium. Indeed, each region of the system relaxes at a rate that follows the local dynamical time. Overall, we see that for radii where the evolution is expected to be linear ($r \le r_{lin}$), the results of the matrix method and the N-body models are in excellent agreement.

## 4. Discussion

### 4.1. Comparison with the Tidal-stripping Scenario

In recent years, models of the evolution of DM subhalos in a tidal field have challenged the historical picture based on tidal heating. Amorisco (2021) and Stücker et al. (2023a), in particular, showed that the remnant subhalo can be accurately modeled by considering that the only effect of the tidal field is to strip a fraction of the system, while the surviving bound remnant is solely evolving due to the loss of these particles, without any direct impact of the tidal field. Here, we reproduce tests of this hypothesis using N-body simulations of a Hernquist DM halo in a galactic tidal field, which we compare with results from LRT.

#### 4.1.1. Numerical Setup for the Tidal N-body Simulation

**N-body model.** We generate an N-body realization with $N = 10^7$ particles of an isotropic Hernquist sphere, as described in Section 2.3.3.

**Host.** We simulate the tidal evolution in an analytical, spherical, and static host system with isothermal potential:

$$\Phi_{\text{host}}(r) = V_0^2 \ln(r/r_0), \qquad (47)$$

where $V_0 \approx 33 V_s$ (see Equation (46)) denotes the circular velocity, and $r_0$ is an arbitrary scale radius. The Hernquist model is placed at apocentre on an orbit with a pericentric distance of $r_{peri} \approx 100 r_s$ and a pericentre-to-apocentre ratio of 1:20. The radial orbital period equals $T_{orb} \approx 24 T_s$ (see Equation (32)). We chose this radial orbit to ensure that the tidally perturbed model has sufficient time available to return to an (approximate) equilibrium configuration at apocentre. As shown in Figure 3, within $\approx 10 T_s$, the perturbed model evolved in isolation has returned to a near-equilibrium state ($\langle v_r \rangle \approx 0$) at those radii where the matrix method is applicable—comparable to the time it takes the tidally perturbed model to return to apocentre after the pericentric passage.

<footnote n="11">¹¹ https://github.com/simrozier/LiRGHaM</footnote>
**Units.** For the sole purpose of illustration, we provide an example rescaling of our model to physical units. For a host circular velocity of $V_c = 220$ km s$^{-1}$ and a pericentre distance of $r_{\text{peri}} = 10$ kpc, we get $r_{\text{apo}} = 200$ kpc and $T_{\text{orb}} \approx 2.3$ Gyr. For the Hernquist model, $M_{\text{tot}} = 4 \times 10^6 M_{\odot}$, $r_s \approx 0.10$ kpc, $V_s \approx 6.6$ km s$^{-1}$, and $T_s \approx 0.093$ Gyr.

**Particle-mesh code.** We use the same simulation setup as in Section 2.3.3, albeit with a simulation box that encloses the entire orbit. We run the simulation for the duration of a full orbital period within the host potential.

### 4.1.2. Initial Conditions for the Matrix Calculation

We use the results of the N-body simulation to inform us about a suitable form of perturbation that may serve as initial conditions to the matrix calculation. For this purpose, we identify the subset of particles located in the final snapshot within a sphere of $5 r_s$ around the tidally stripped N-body model. This value is chosen purely empirically to approximately define what we consider the virialized remnant system after the pericentric passage, excluding tidal tails and other disequilibrium features.

To define initial conditions for the matrix calculation, we trace back these selected particles to the simulation's initial conditions. We note that their distribution does not satisfy the hypotheses of Equation (3), because their distribution has a residual phase dependence and velocity anisotropy. However, we can artificially phase mix and isotropize it with the following procedure.<fnref n="12" /> First, we phase mix the distribution by computing the differential energy distribution $dN/dE$ of this subset of particles in the initial Hernquist model. Then, we isotropize it by defining the truncated DF using Equation (28). The truncated DF and differential energy distribution are shown using blue filled circles in Figure 1. Finally, we fit the truncated DF using the empirical formula Equation (30), finding the parameters $\{E_{\text{tr}}/E_0\}$, $A/E_0\} = \{-0.176, 0.077\}$ to best match the N-body results. These initial conditions are identical to the ones used previously as an example in Section 2.3 (model 1; see Table 1); the matrix calculation we compare in the following against the N-body run in the presence of a tidal field is therefore completely identical to the model shown in Figure 2.

### 4.1.3. Comparison of the Matrix Calculation with a Tidal N-body Simulation

We are now able to compare the final equilibrium state as computed using the matrix method with the tidally stripped N-body model at apocentre. Figure 4 shows density $\rho(r)$, radial velocity dispersion $\sigma_r(r)$, and anisotropy parameter $\beta(r)$ for the evolved N-body model (filled circles) and the matrix calculation (solid curve), with units and notations as in Figure 2. For reference, we also show the initial N-body Hernquist realization when injected into the host potential. Note that the initial disequilibrium profile used for the matrix calculation is identical to the one shown already in Figure 2 and is omitted here for clarity.

In the inner regions ($r \lesssim r_s$), the simple empirical energy truncation to the initial conditions, evolved to equilibrium

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison of model 1 (see Figure 2 and Table 1) against an N-body realization of a Hernquist profile evolved in the presence of a tidal field. Format and units are identical to Figure 2. The “initial” N-body snapshot shows the model immediately after being injected into the tidal field. The model is injected at the apocentre of an orbit with a pericentre-to-apocentre ratio of 1:20. The “final” snapshot shows the N-body again at apocentre after a full orbital period ($T_{\text{orb}} \approx 24.2 T_s$). In the inner regions, the matrix calculation matches that of the tidally stripped N-body model. Deviations are visible in the outer regions, as discussed in the text.</figcaption></figure>

using the matrix method, matches the density and radial velocity dispersion as measured in the tidal N-body simulation. The results of the N-body simulation are consistent with a centrally isotropic velocity distribution as predicted by the matrix method. The outer anisotropy profile is mildly tangentially biased for both the matrix and the N-body model.

Beyond $r_s$, density and radial velocity dispersion as computed using the matrix method start to deviate from the profiles measured in the N-body simulation. Multiple effects may play a role in shaping the outer density profile of the tidally stripped system which are not captured by our simplistic model consisting of an energy truncation and subsequent return to equilibrium. Here, we assume that stripped particles are removed instantaneously from the system, and that the sole response of the system can be captured by its return to equilibrium in the absence of the stripped particles. Our model does not include (i) particles which were stripped and later recaptured, (ii) energy injection/tidal heating, which may be

<footnote n="12">¹² Enforcing initial conditions with isotropic kinematics is a modeling choice we made for the sake of simplicity in this first application of the matrix method to tidal stripping. The matrix method can generally be applied also to anisotropic initial conditions.</footnote>
important in particular for those particles energetically close to
the truncation shown in Figure 1, (iii) any angular momentum
dependence when selecting stripped and remaining particles in
the initial conditions, (iv) nonlinear evolution of particles
beyond $r > r_{\text{lin}}$.

## 4.2. Depletion of the Central Density Cusp

In this section, we apply the matrix method to model the
tidal evolution of the central density cusp of cold DM subhalos.
The detailed properties of the central density profile of DM
subhalos are of particular relevance for the study of potential
signals of DM self-annihilation and decay, which scale with
$\rho^2(r)$ and $\rho(r)$, respectively (see, e.g., Stref et al. 2019;
Facchinetti et al. 2022; Stücker et al. 2023b; Delos & White
2023; Lovell 2024). Models of the expected signals need to
take into account the evolution of subhalo DM density profiles
in the tidal field of the Milky Way. Studies using controlled N-
body simulations suggest a depletion of the central density
(e.g., Hayashi et al. 2003; Peñarrubia et al. 2010; Green & van
den Bosch 2019; Errani & Navarro 2021), thereby attenuating
any potential emissions from DM self-annihilation.

The interpretation of N-body simulation results, however, is
complicated by the demanding resolution requirements to
accurately model the tidal survival of dark substructures.
Insufficient (temporal, spatial, and particle) resolution of N-
body simulations is known to heavily impact the tidal evolution
of N-body subhalos, affecting both their structure and their
overall abundance (see, e.g., van den Bosch & Ogiya 2018,
who coined the term artificial disruption for the resolution-
driven depletion of dark substructures in simulations). Limited
spatial resolution is shown to lead to the formation of a
constant-density core even in initially cuspy N-body realiza-
tions of cold DM halos (Errani & Peñarrubia 2020). Tidal
evolution further amplifies this issue: Insufficient spatial and/or
particle resolution leads to a systematic underestimation of a
subhalo's characteristic density, which in turn renders the
subhalo even more vulnerable to the effects of tides (Errani &
Navarro 2021). Given these resolution challenges, it is no
surprise that the depletion of the central density in N-body
realizations of cuspy DM halos has been suspected to be a mere
numerical artifact, in particular in light of recent analytical
approaches predicting virtually unaltered central densities in
such halos (Drakos et al. 2020, 2022; Stücker et al. 2023a).

Clearly, the matrix method is not subject to the complica-
tions arising from discreteness noise, and hence particularly
suitable to study the centers of DM subhalos. The matrix
method resolves orbits within the subhalo essentially all the
way to the smallest timescales involved. There is, however, a
spatial resolution limit intrinsic to the matrix method, arising
from the choice of basis. Therefore, we perform convergence
tests where we gradually increase the number of basis elements
to understand the systematic impact on the central density
profile. The results of these tests are shown in Figure 5. We plot
the central density profile in the final equilibrium for model 2
(orange) and model 3 (red); see Table 1 for parameters. The
profiles are normalized by the initial Hernquist density at all
radii. We also plot the initial disequilibrium density profiles as
orange and red dashed curves (for model 2 and model 3,
respectively), which both converge to the initial Hernquist
density for $r \to 1$.

We gradually increase the number of basis functions used in
the matrix computation, ranging from $n_{\text{max}} = 50$ to

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Convergence study. Final density profiles computed using the matrix method for gradually increasing spatial resolution, i.e., number of basis elements $n_{\text{max}} = 50, 75, 100$, and 300. The results obtained for model 2 and model 3 (see Table 1) are shown in orange and red, respectively. The profiles are normalized by the original Hernquist density at each radius. For reference, we also plot the corresponding initial disequilibrium density profiles as dashed curves. Radii below the resolution limit $r_{\text{bas}}$ are shown in gray. Within the resolved range of radii, the inner regions of the normalized profiles are well approximated by linear functions. To guide the eye, these are shown as dotted black lines. This convergence test suggests that collisionless relaxation decreases the amplitude of the perturbed Hernquist cusp, even though the initial energy-truncated profile hardly affected the inner regions.</figcaption></figure>

$n_{\text{max}} = 300$. The results of this experiment are shown as solid curves in Figure 5. In each curve the resolution limit is taken to be $r_{\text{bas}} = 3 R_{\text{b}}/n_{\text{max}}$, this time with $R_{\text{b}} = 6 r_{\text{s}}$ because it appears to better reconstruct the perturber. Because of the *cored* nature of the basis (every element of the density basis has a finite value at $r=0$), each curve converges to 1 at $r=0$ by construction, i.e., our choice of basis forces the density of the final state at $r=0$ (below the resolution limit $r_{\text{bas}}$) to be identical to the central density of the initial Hernquist model. However, at those radii which we consider well resolved, the density in the final state is always less than in the initial energy-truncated disequilibrium profile. For both *model 2* and *model 3*, the final density in the resolved region is well approximated by linear fits ($\rho(r)/\rho_{\text{H}}(r) \approx 0.98 - 0.12 r/r_{\text{s}}$ for *model 2* and $\rho(r)/\rho_{\text{H}}(r) \approx 0.96 - 0.22 r/r_{\text{s}}$ for *model 3*), shown in Figure 5 using back dotted lines.

This convergence test has important implications on how we
can interpret the effect of the imposed mass loss. First, the
dominant term in the central density is still a cusp with slope
$d \ln \rho/d \ln r = -1$, so the mass loss did not modify the nature
of the central cusp. Second, the amplitude of the cusp has
decreased, meaning that the revirialization phase was able to
displace an amount of material with divergent central density
($\rho^s \sim 0.02 \rho_{\text{H}}(r) \propto r^{-1}$ and $\rho^s \sim 0.04 \rho_{\text{H}}(r)$ for *model 2* and *model 3*, respectively), even though the perturber $\rho^e$ has a
finite density at the center. Coupling this observation with the
discussion of Section 4.1, our results appear consistent with
those obtained with high-resolution N-body simulations of DM
subhalos in smooth tidal fields (see, e.g., Green & van den
Bosch 2019; Errani & Navarro 2021): Revirialization decreases
the cusp's amplitude, however it does not modify its slope.<fnref n="13" />

Finally, this test case also helps us learn about dynamical
studies based on basis function expansions in general. In the

<footnote n="13">¹³ Note that here we refer to a progressive tidal-stripping process along many pericentric passages in a smooth tidal field. Extreme tidal shocks cannot be modeled with the matrix method as discussed here.</footnote>
matrix method, the quality of the choice of basis is evaluated in light of how well the perturber and the response are reconstructed with as few basis functions as possible. If a basis set is taken to be biorthogonal, then usually the central behavior will be the same for all radial density basis elements. Our problem here is that the central behavior of the perturber and of the response are very different: The perturber is finite at the center, while the response diverges as $r^{-1}$. Consequently, there is no single biorthogonal basis that is perfectly suited to reconstruct both the perturber and the response. Here, we chose a basis that is finite at the center, which makes it difficult to detect the response's central divergence; had we chosen a cuspy basis set (e.g., Hernquist & Ostriker 1992; Lilley et al. 2018), our problem would transform into accurately projecting the perturber. It is likely that the problem persists for tailor-made bases, such as the ones from Weinberg (1999), Petersen et al. (2022), Lilley & van de Ven (2023), and Petersen et al. (2024).

## 5. Conclusions

In the present work, we develop an analytical, linear model for the collisionless relaxation of self-gravitating systems from near equilibrium configurations. We apply it to a spherical Hernquist model, which we perturb gently by removing a fraction of the system in the initial conditions. We choose a perturbation in energy, modeled to approximate the type of perturbations typical of weak tidal encounters. To model the return to equilibrium, we solve the linearized system of Collisionless Boltzmann-Poisson equations using a formalism rooted in the Kalnajs (1977) matrix method, which we test against controlled N-body simulations. Below, we summarize our main results.

1. We show that the final equilibrium state can be computed directly (i.e., in a single computational step) using the matrix method, at very modest computational cost. The equilibrium density profiles agree at sub-percent level with the results of controlled N-body simulations.

2. The formalism presented here allows us to compute the full phase-space DF, and thereby the velocity anisotropy profile in the final equilibrium state. Comparison against velocity dispersion- and anisotropy profiles measured in controlled N-body simulations shows agreement at percent level.

3. The matrix formalism can also be used to predict the full time evolution from the initial disequilibrium to the final equilibrium, including all intermediate evolutionary stages. The method allows to model the evolution of the density profile, as well as moments of the velocity distribution and derived quantities (such as mean radial velocity, radial velocity dispersion and anisotropy profile). Also here, the evolution predicted by the linear theory is consistent with the results of controlled N-body models.

4. We compare the results of the matrix calculation also against N-body simulations that include the time-varying tidal field as experienced on a radial orbit in an isothermal host potential. We show that a simple perturbation in energy space, used as initial condition to the matrix method, allows to predict final equilibrium density and velocity dispersion profiles as measured in the N-body simulations at radii below the initial Hernquist scale

radius. At larger radii, the matrix method and the N-body simulations remain qualitatively in agreement and predict a slight tangential bias to the velocity dispersion, though the degree of anisotropy measured in the N-body models is larger than the value predicted by our linear theory. As the initial disequilibrium configuration studied here has isotropic kinematics, the tangential anisotropy in the final state is a consequence of the collisionless relaxation alone.

5. The matrix model predicts a decrease in the amplitude of the central density cusp of the tidally perturbed Hernquist sphere, in agreement with the results of previous simulation studies. Crucially, the depletion of the central density is driven by the collisionless relaxation, and not by the initial tidal energy truncation.

The methods developed in this work constitute a computationally highly efficient framework to study the response of gravitational systems to (small) perturbations. We believe the methods may find application in numerous astrophysical contexts.

In the context of tidal stripping, the methods may be used to efficiently explore the systematics in the response to tidal perturbations for different families of initial DFs, such as those of cuspy and cored models with different anisotropy profiles, or rotation. Iterative application of the matrix method may allow to approximate the response beyond the linear regime. We will explore these topics in a future contribution.

Furthermore, the method may be applied to study perturbations to the central regions of stellar systems, such as central mass accumulation due to black hole growth in star clusters (e.g., Young 1980), or gas accretion onto galaxies (e.g., Sellwood & McGaugh 2005). Similarly, the method can be used to model the response to central mass depletion due to stellar feedback, which applies both to the study of cusp-core transformations in dwarf galaxies (Freundlich et al. 2020), or the (collisionless) response to stellar escape in star clusters.

Finally, the analytical developments presented here challenge the adiabatic versus impulsive distinction that is often made in the analysis of the response of stellar systems to (tidal) perturbations. Our results show that in the linear approximation, the time asymptotic equilibrium does not depend on how the perturbation is applied, but solely on its time asymptotic value. These systematics may be fundamental to the emergence of virtually identical tidal evolutionary tracks for cuspy stellar systems subject to impulsive- or adiabatic tidal mass removal.

## Acknowledgments

The authors would like to thank Anna Lisa Varri and Benoit Famaey for insightful discussions. S.R. acknowledges support from the Agence Nationale de la Recherche (ANR project GaDaMa ANR-18-CE31-0006) and from the Royal Society (Newton International Fellowship grant No. NIF\R1\221850). R.E. acknowledges support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement number 834148), and from the National Science Foundation (NSF) grant AST-2206046. Furthermore, the authors gratefully acknowledge the access to computational resources provided by the University of Victoria.
# Appendix A
## Velocity Moments of a Spherical, Nonrotating System

### A.1. Simplifying the Matrix

In the case of a spherical potential, the angle Fourier-transformed functions $[\mu\psi^{(p)}]_n(\mathbf{J})$ and $\psi_n^{(q)}(\mathbf{J})$ appearing in Equation (44) can be simplified, following developments from Polyachenko & Shukhman (1981), Fridman et al. (1984), and Tremaine & Weinberg (1984). First, we can explicitly develop a basis element as a spherical harmonic times a radial basis function, so that

$$[\mu\psi^{(p)}]_n(\mathbf{J}) = \int \frac{d\theta}{(2\pi)^3} \mu(\mathbf{v}) U_{n\ell}^{\ell p}(r) Y_{\ell\ell'}^{m\ell'}(\theta, \phi) e^{-i\mathbf{n}\cdot\mathbf{\theta}}. \quad (A1)$$

Here, let us highlight that $\theta$ is the polar angle in spherical coordinates, while $\theta = (\theta_1, \theta_2, \theta_3)$ are the angle phase-space variables. We can now use the fact that in a spherical potential, a phase-space point with coordinates $(r, \theta, \phi, \mathbf{v})$ verifies

$$\begin{aligned} Y_{\ell}^{m}(\theta, \phi) &= e^{im\theta_3} Y_{\ell}^{m}(\theta, \phi - \theta_3) \\ &= e^{im\theta_3} \sum_{k=-\ell}^{\ell} \mathcal{R}_{km}^{\ell}(\beta) y_{\ell}^{k} e^{ik\psi} i^{m-k}, \end{aligned} \quad (A2)$$

where $\cos(\beta) = L_z/L$, $\mathcal{R}$ is Wigner's rotation matrix, $y_{\ell}^k = Y_{\ell}^k(\frac{\pi}{2}, 0)$ and $\psi$ is the angle measured in the orbital plane from the ascending node to the current position. At that point, useful simplifications arise when we assume that $\mu(\mathbf{v})$ is only a function of the angle-action variables through the actions $\mathbf{J}$ and the radial angle $\theta_1$. In a spherical potential, this is the case if $\mu$ is a function of $v_r$ and $v_t$, so it applies to all moments of the radial and tangential velocities. Indeed, when we have $\mu(\mathbf{J}, \theta_1)$, then

$$\begin{aligned} [\mu\psi^{(p)}]_n(\mathbf{J}) = & \sum_{k=-\ell^p}^{\ell^p} \mathcal{R}_{km}^{\ell^p}(\beta) y_{\ell^p}^k i^{m^p-k} \\ & \times \int \frac{d\theta_1}{2\pi} \mu(\mathbf{J}, \theta_1) U_{n\ell}^{\ell p}(r) e^{-in_1\theta_1} \\ & \times \int \frac{d\theta_2}{2\pi} \frac{d\theta_3}{2\pi} e^{-i(n_2\theta_2+n_3\theta_3)} e^{im\theta_3} e^{ik\psi}. \end{aligned} \quad (A3)$$

Now, because $\psi$ does not depend on $\theta_3$, we can carry out the integration w.r.t. $\theta_3$, yielding $\delta_{n_3}^{m\ell}$. In addition, we can use the fact that $(\psi - \theta_2)$ is a function of $\theta_1$ only to rewrite

$$\begin{aligned} [\mu\psi^{(p)}]_n(\mathbf{J}) = & \delta_{n_3}^{m\ell} \sum_{k=-\ell^p}^{\ell^p} \mathcal{R}_{km}^{\ell^p}(\beta) y_{\ell^p}^k i^{m^p-k} \\ & \times \int \frac{d\theta_1}{2\pi} \mu(\mathbf{J}, \theta_1) U_{n\ell}^{\ell p}(r) e^{-i[n_1\theta_1 + k(\theta_2 - \psi)]} \\ & \times \int \frac{d\theta_2}{2\pi} e^{-i(n_2-k)\theta_2}. \end{aligned} \quad (A4)$$

Again, the $\theta_2$ integral gives $\delta_{n_2}^k$, which restricts the sum over $k$ to a single term. Let $\tilde{n} = (n_1, n_2)$, $\tilde{\mathbf{J}} = (J_r, L)$, and $g$ be a function of $(\tilde{\mathbf{J}}, \theta_1)$, we can define two operators,<fnref n="14" /> $W\tilde{n}[g]$ and

$Z^{\tilde{n}}[g]$, equal to

$$W^{\tilde{n}}[g](\tilde{\mathbf{J}}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} d\theta_1 g(\tilde{\mathbf{J}}, \theta_1) \cos[n_1\theta_1 + n_2(\theta_2 - \psi)], \quad (A5)$$

$$Z^{\tilde{n}}[g](\tilde{\mathbf{J}}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} d\theta_1 g(\tilde{\mathbf{J}}, \theta_1) \sin[n_1\theta_1 + n_2(\theta_2 - \psi)]. \quad (A6)$$

With these definitions, we have

$$\begin{aligned} [\mu\psi^{(p)}]_n(\mathbf{J}) = & \delta_{n_3}^{m\ell} \mathcal{R}_{n_2m}^{\ell\ell'}(\beta) y_{\ell\ell'}^{n_2} i^{m\ell'-n_2} \\ & \times (W^{\tilde{n}}[\mu U_{n\ell}^{\ell\ell'}(\tilde{\mathbf{J}})] - i Z^{\tilde{n}}[\mu U_{n\ell'}^{\ell\ell'}(\tilde{\mathbf{J}})]). \end{aligned} \quad (A7)$$

From Equations (A5) and (A6), we see that because $\theta_2 - \psi$ is odd in $\theta_1$, then $W^{\tilde{n}}[g] = 0$ (resp. $Z^{\tilde{n}}[g] = 0$) for any $g$ that is odd (resp. even) in $\theta_1$. Decomposing $g = g_+ + g_-$ in even and odd parts, we therefore have

$$W^{\tilde{n}}[g](\tilde{\mathbf{J}}) = \frac{1}{\pi} \int_{0}^{\pi} d\theta_1 g_{+}(\tilde{\mathbf{J}}, \theta_1) \cos[n_1\theta_1 + n_2(\theta_2 - \psi)], \quad (A8)$$

$$Z^{\tilde{n}}[g](\tilde{\mathbf{J}}) = \frac{1}{\pi} \int_{0}^{\pi} d\theta_1 g_{-}(\tilde{\mathbf{J}}, \theta_1) \sin[n_1\theta_1 + n_2(\theta_2 - \psi)]. \quad (A9)$$

We also note that they verify the symmetry relations

$$W^{-\tilde{n}} = W^{\tilde{n}}; \quad Z^{-\tilde{n}} = -Z^{\tilde{n}}. \quad (A10)$$

Finally, since $r$ is an even function of $\theta_1$, we have

$$\begin{aligned} [\mu\psi^{(p)}]_n(\mathbf{J}) = & \delta_{n_3}^{m\ell} \mathcal{R}_{n_2m}^{\ell\ell'}(\beta) y_{\ell\ell'}^{n_2} i^{m\ell'-n_2} \\ & \times (W^{\tilde{n}}[\mu_+ U_{n\ell}^{\ell\ell'}(\tilde{\mathbf{J}})] - i Z^{\tilde{n}}[\mu_- U_{n\ell'}^{\ell\ell'}(\tilde{\mathbf{J}})]). \end{aligned} \quad (A11)$$

We conveniently recover Equation (A4) of Rozier et al. (2019) by taking $\mu(\mathbf{v}) = 1$.

Further simplifications are obtained when we assume that $F_i$ is a function of the actions only through $\tilde{\mathbf{J}}$, implying that the system is nonrotating. Then, the only $L_z$ dependence in the integrand of Equation (44) resides in the $\mathcal{R}_{n_2m}^{\ell\ell}(\beta) \mathcal{R}_{n_2m}^{\ell\ell'}(\beta)$ terms from the angular Fourier transforms. Owing to the orthogonality relations of the Wigner rotation matrices, the $L_z$ integral can then be simplified. Thus, the response matrix can be rewritten

$$\mathcal{M}[\mu]_{pq}(t) = \delta_{m\ell}^{n\ell} \delta_{\ell\ell'}^{q\ell} \sum_{\tilde{n}} C_{\ell\ell'}^{n\ell} P[\mu]_{\ell\ell'\nu\nu'}^{\tilde{n}}(t), \quad (A12)$$

where the coefficients $C$ are given by

$$C_{\ell}^{n} = 2(2\pi)^3 \frac{|y_{\ell}^{n}|^2}{2\ell + 1}, \quad (A13)$$

and the functions $P$ by

$$\begin{aligned} P[\mu]_{\ell\ell'\nu\nu'}^{\tilde{n}}(t) = & \int d\tilde{\mathbf{J}} e^{-i\tilde{n}\cdot\tilde{\Omega}t} L \tilde{n} \cdot \frac{\partial F_i}{\partial \tilde{\mathbf{J}}} W^{\tilde{n}}[U_{n\ell}^{\ell\ell'}] \\ & \times (Z^{\tilde{n}}[\mu_- U_{n\ell}^{\ell\ell'}] - i W^{\tilde{n}}[\mu_+ U_{n\ell}^{\ell\ell'}]). \end{aligned} \quad (A14)$$

A final computational speedup is achieved by noticing that

$$C_{\ell}^{-n} = C_{\ell}^{n}; \quad P^{-\tilde{n}} = (P^{\tilde{n}})^*, \quad (A15)$$

so that the sum in Equation (A12) can be easily reduced to half of the possible $\tilde{n}$ pairs.

<footnote n="14">¹⁴ Tremaine & Weinberg (1984) only define the $W$ operator, because they focus on a choice of $g$ that is merely a function of $r$.</footnote>
## A.2. Numerical Implementation

Here, we summarize the numerical implementation of the response matrix computation, referring the reader to Rozier et al. (2019) for more details.

The numerical effort required to compute the response matrix of Equation (A12) is essentially contained in the computation of the functions P of Equation (A14). Three specific difficulties must be highlighted: (i) performing the action-space integration; (ii) computing the orbital characteristics (actions, frequencies) of each orbit evaluated; (iii) computing the angle-space integrals of Equations (A8) and (A9).

**Action-space integration.** In order to simplify the computation of the orbital elements, we change variables of the action-space integrals to the orbital peri- and apocentre ($r_p, r_a$). Because of the structure of the integrand, we perform an extra change of variables into a new set $(u, v)$ that is denser in two regions: when $r_a \ll r_s$ (central region) and when $r_a - r_p \ll r_p$ (nearly circular orbits). For numerical efficiency, $r_p(u, v)$ and $r_a(u, v)$ are explicit functions that are easy to compute and derive. Hence, if we formally write Equation (A14) as

$$P = \int d\tilde{J} p(\tilde{J}), \quad (A16)$$

then we re-express the integral so that

$$P = \int du dv \left| \frac{\partial \tilde{J}}{\partial r_p, r_a} \right| \left| \frac{\partial r_p, r_a}{\partial u, v} \right| p(\tilde{J}(r_p(u, v), r_a(u, v))), \quad (A17)$$

where we included the respective Jacobian determinants of the changes of variables.

Then, the integration on the $(u, v)$ plane is performed via a tailor-made method that is adapted to the specific structure of the integrand. In more detail, the integral can be rewritten as

$$P = \int du dv g(u, v) e^{i(h(u,v))}, \quad (A18)$$

where both $g$ and $h$ are slowly varying functions of their arguments. Therefore, we divide the $(u, v)$ plane in small squares of side $\Delta u$ centered around grid points $(u_i, v_i)$, and within each square we perform a bilinear expansion of $g$ and $h$ around $(u_i, v_i)$. The resulting integrand is

$$\int_{-\frac{\Delta u}{2}}^{\frac{\Delta u}{2}} dx \int_{-\frac{\Delta u}{2}}^{\frac{\Delta u}{2}} dy (g + g_u x + g_v y) e^{i(h+h_u x+h_v y)}, \quad (A19)$$

where the subscripts $u$ and $v$ denote respective derivatives w.r.t. $u$ and $v$, and all functions are evaluated at the point $(u_i, v_i)$. This integral has an analytical expression as a function of $g(u_i, v_i)$, $g_u(u_i, v_i)$, $g_v(u_i, v_i)$, $h(u_i, v_i)$, $h_u(u_i, v_i)$, $h_v(u_i, v_i)$. Once these six evaluations are made for each grid point, we can simply add the contribution from each square to compute the total action-space integral. All our results are checked for convergence w.r.t. the grid step size, $\Delta u$.

**Orbital characteristics.** In order to compute the values of $g$, $h$ and their derivatives on the $(u, v)$ grid, one needs to specify how the orbital characteristics (actions, frequencies and their derivatives) are computed. At each grid point, we straightforwardly obtain the values of $(r_p(u_i, v_i), r_a(u_i, v_i))$, from which all other functions can be derived. Specifically, we obtain the

orbital energy and angular momentum through

$$E = \frac{r_a^2 \psi_0(r_a) - r_p^2 \psi_0(r_p)}{r_a^2 - r_p^2}; \quad L = \sqrt{\frac{2(\psi_0(r_a) - \psi_0(r_p))}{r_p^{-2} - r_a^{-2}}}, \quad (A20)$$

the radial action through

$$J_r = \frac{1}{\pi} \int_{r_p}^{r_a} dr \sqrt{2(E - \psi_0(r)) - L^2/r^2}, \quad (A21)$$

and the frequencies through

$$\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}, \quad (A22)$$

$$\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}}. \quad (A23)$$

To perform the radial integrals, we first change variables to an orbital anomaly $x$, i.e., any integral becomes

$$Y(r_p, r_a) = \int_{r_p}^{r_a} dr \, y(r, r_p, r_a) = \int_{-1}^{1} dx \frac{dr}{dx} y(r(x), r_p, r_a). \quad (A24)$$

We tune the change of variables so that all integrals run from $-1$ to $1$ and any integrable divergence is cured. Then, any derivative of $Y$ w.r.t. $r_p$ or $r_a$ can be straightforwardly transported under the integral sign, so that they can also be written in that integral form. Then, we perform the integration through a Gauss-Kronrod method.

**Angle-space integrals.** The angle-space integral requires particular care, because it is the bottleneck of the whole matrix evaluation. Let us focus on Equation (A8). Again, the most direct way of computing the integral is to change orbital variables from radial angle to radius, because we do not have an explicit expression for the azimuthal anomaly $\theta_2 - \psi$ as a function of $\theta_1$. Hence, we rewrite the function $W$ as

$$W = \int_{0}^{\pi} d\theta_1 w(\theta_1, \theta_2) = \int_{r_p}^{r_a} dr \frac{d\theta_1}{dr} w(\theta_1(r), \theta_2(r)). \quad (A25)$$

Here, the angles have an explicit integral expression as a function of $r$, as

$$\theta_1(r) = \int_{r_p}^{r} dr' \frac{\Omega_1}{\sqrt{2(E - \psi_0(r')) - L^2/r'^2}}, \quad (A26)$$

$$(\theta_2 - \psi)(r) = \int_{r_p}^{r} dr' \frac{\Omega_2 - L/r'^2}{\sqrt{2(E - \psi_0(r')) - L^2/r'^2}}. \quad (A27)$$

Because the integration boundaries of $\theta_1$, $\theta_2$ match the integration variable of $W$, the integral for $W$ can be computed as a single 1D integral instead of a 2D nested integral, via a method detailed in Rozier et al. (2019). Similarly to the orbital characteristics, the integrals over $r$ are performed via a change of variables with an orbital anomaly $x$, which also helps computing the derivatives of $W$ w.r.t. $r_p, r_a$.

## A.3. Numerical Complexity

We quantify here the numerical cost of the matrix method. We first focus on one of the runs described in Section 2.3.4, and we give scaling relations that help estimating the complexity of other runs in the paper.
Computing the response matrix for model 1 in Figure 2 involved $2 \times 10^{15}$ floating-point operations, and was one of the most involved runs of the paper, as it required both a sufficient dynamic range in radii and resolution in the center. The radial resolution can be modified by tuning $n_{max}$, the numerical complexity changing as $n_{max}^4$. This dependence is due to the fact that $n_{max}$ impacts both the size of the response matrix (as $n_{max}^2$) and the resolution of the action-space integration grid (as $n_{max}^2$).

In the runs of Section 3.4, the number of operations is impacted in an affine way by the number of time steps. By lowering the radial resolution (through $n_{max}$) and the maximum order of the angular Fourier number $n_1$, we were able to keep the number of floating-point operations at around $2 \times 10^{15}$ as well.

## Appendix B
### Edge Integrals

#### B.1. Definition

In our case of interest, a significant contribution to the response matrix can come from integrating Equation (35) at the edges of the integration domain.<fnref n="15" /> As shown in Jalali & Hunter (2005; see also Polyachenko & Shukhman 2015), the edge integrals can be considered by replacing $F_i$ in the integrand with $F_i \times \mathcal{I}(J)$, where $\mathcal{I}$ is the indicator function of the integration domain (yielding 1 inside and 0 outside). Consequently, a term has to be added to $M_{pq}(t)$, equal to

$$ -i(2\pi)^3 \sum_n \int d\mathbf{J} e^{-i\mathbf{n} \cdot \Omega t} F_i(\mathbf{J}) \mathbf{n} \cdot \frac{\partial \mathcal{I}}{\partial \mathbf{J}} \psi_n^{(p)*}(\mathbf{J}) \psi_n^{(q)}(\mathbf{J}). \quad (B1) $$

In practice, as explained in Section A.2, the action-space integral reduces to a 2D integral over $\tilde{J} = (J_r, L)$, and it is performed via successive changes of variables, from $J_r, L$ through $E, L$, then $r_p, r_a$ and eventually $u, v$. Hence, the function $\mathcal{I}$ is defined in terms of these $u, v$ variables. Simple but tedious algebra helps us compute the derivatives of $\mathcal{I}$ w.r.t. the actions from its derivatives w.r.t. $u, v$, yielding

$$ \tilde{\mathbf{n}} \cdot \frac{\partial \mathcal{I}}{\partial \tilde{\mathbf{J}}} = \left| \frac{\partial E, L}{\partial r_p, r_a} \right|^{-1} \left| \frac{\partial r_p, r_a}{\partial u, v} \right|^{-1} \times \left( \tilde{\mathbf{n}} \cdot \tilde{\Omega} \left[ \left( \frac{\partial \mathcal{I}}{\partial u} \frac{d(r_a - r_p)}{dv} - \frac{\partial \mathcal{I}}{\partial v} \frac{dr_p}{du} \right) \frac{\partial L}{\partial r_a} - \frac{\partial \mathcal{I}}{\partial v} \frac{dr_p}{du} \frac{\partial L}{\partial r_p} \right] - n_2 \left[ \left( \frac{\partial \mathcal{I}}{\partial u} \frac{d(r_a - r_p)}{dv} - \frac{\partial \mathcal{I}}{\partial v} \frac{dr_p}{du} \right) \frac{\partial E}{\partial r_a} - \frac{\partial \mathcal{I}}{\partial v} \frac{dr_p}{du} \frac{\partial E}{\partial r_p} \right] \right). \quad (B2) $$

Here, the two first terms in the right-hand side are the Jacobian determinants of the transformations $(r_p, r_a) \rightarrow (E, L)$ and $(u, v) \rightarrow (r_p, r_a)$, and we made use of the fact that in our specific change of variables, $r_a(u, v) = r_p(u) + (r_a - r_p)(v)$, so that $\partial r_p / \partial v = 0$ and $\partial r_a / \partial u = \partial r_p / \partial u$. It therefore remains to specify the derivatives of $\mathcal{I}$ w.r.t. $u, v$ at the edges of the integration domain.

Four cases arise, depending on the specific edge considered: (i) $u = u_{min}(v)$, where locally $\mathcal{I} = \Theta(u - u_{min}(v))$, (ii) $u = u_{max}(v)$, where $\mathcal{I} = \Theta(u_{max}(v) - u)$, (iii) $v = v_{min}(u)$,

where $\mathcal{I} = \Theta(v - v_{min}(u))$, and (iv) $v = v_{max}(u)$, where $\mathcal{I} = \Theta(v_{max}(u) - v)$. More details on the way these edges relate to the input parameters can be found in Rozier et al. (2019), we summarize here their main characteristics.

(i) $u = u_{min}(v)$ is reached when $r_p(u)$ is maximal at fixed $v$, corresponding to when $r_a$ is maximal: $r_a(u_{min}(v), v) = r_{max}$. So $u_{min}(v) = r_p^{-1}(r_{max} - (r_a - r_p)(v))$. As a consequence,

$$ \frac{du_{\min}}{dv} = - \frac{\frac{d(r_a - r_p)}{dv}}{\frac{dr_p}{du}(u_{\min})}. \quad (B3) $$

Hence, on that edge, we have

$$ \begin{aligned} \frac{\partial \mathcal{I}}{\partial u} &= \delta(u - u_{\min}(v)); \\ \frac{\partial \mathcal{I}}{\partial v} &= \frac{\frac{d(r_a - r_p)}{dv}}{\frac{dr_p}{du}(u_{\min})} \delta(u - u_{\min}(v)). \end{aligned} \quad (B4) $$

(ii) $u = u_{max}(v)$ is reached when $r_p(u) = r_{min}$. Hence,

$$ \frac{du_{\max}}{dv} = 0, \quad (B5) $$

and we have

$$ \begin{aligned} \frac{\partial \mathcal{I}}{\partial u} &= -\delta(u_{\max}(v) - u); \\ \frac{\partial \mathcal{I}}{\partial v} &= 0. \end{aligned} \quad (B6) $$

(iii) $v = v_{min}(u)$ is reached when $r_a(u, v_{min}(u)) = r_{max}$. So $v_{min}(u) = (r_a - r_p)^{-1}(r_{max} - r_p(u))$. As a consequence,

$$ \frac{dv_{\min}}{du} = - \frac{\frac{dr_p}{du}}{\frac{d(r_a - r_p)}{dv}(v_{\min})}. \quad (B7) $$

Hence, on that edge, we have

$$ \begin{aligned} \frac{\partial \mathcal{I}}{\partial u} &= \frac{\frac{dr_p}{du}}{\frac{d(r_a - r_p)}{dv}(v_{\min})} \delta(v - v_{\min}(u)); \\ \frac{\partial \mathcal{I}}{\partial v} &= \delta(v - v_{\min}(u)). \end{aligned} \quad (B8) $$

(iv) Finally, $v_{max}$ is reached for $(r_a - r_p)(v) = r_{min}$, so that

$$ \begin{aligned} \frac{\partial \mathcal{I}}{\partial u} &= 0; \\ \frac{\partial \mathcal{I}}{\partial v} &= -\delta(v_{\max}(u) - v). \end{aligned} \quad (B9) $$

#### B.2. Numerical Implementation

On each of the edges of the $(u, v)$ integration domain, we have to perform a 1D integral of the form

$$ \mathcal{J} = \int dv g(u_{\min}(v), v) e^{i h(u_{\min}(v), v)}. \quad (B10) $$

In the same way as our method for computing 2D integrals with a fast trigonometric oscillation, we cut the integration domain

<footnote n="15">¹⁵ The exact same methods are applied to compute the time asymptotic matrix of Equation (21), or the moment response matrix of Equation (44).</footnote>
in small segments of width $\Delta v$ centered on successive values $v_i$. Then, we expand the slowly varying functions $g, h$ to first order in $v$ around $v_i$, so that we approximate

$$ 
\begin{aligned}
\mathcal{J} &= \sum_{v_i} \int_{-\frac{\Delta v}{2}}^{\frac{\Delta v}{2}} dy \left[ g + \left( g_u \frac{du_{\min}}{dv} + g_v \right) y \right] e^{i \left[ h + \left( h_u \frac{du_{\min}}{dv} + h_v \right) y \right]} \\
&= \sum_{v_i} g e^{ih} \Re \left( \frac{g_u \frac{du_{\min}}{dv} + g_v}{g}, h_u \frac{du_{\min}}{dv} + h_v, \Delta v \right), \quad (B11)
\end{aligned}
$$

where all functions are evaluated at $(u_{\min}(v_i), v_i)$ and $\Re$ is given by

$$ \Re(a, b, \Delta v) = \int_{-\frac{\Delta v}{2}}^{\frac{\Delta v}{2}} dv (1 + av) e^{ibv}. \quad (B12) $$

This can be further simplified by considering the integral with normalized bounds, i.e.,

$$ \Re(a, b, \Delta v) = \Delta v \Re_D(a\Delta v, b\Delta v), \quad (B13) $$

where

$$ \Re_D(\alpha, \beta) = \int_{-1/2}^{1/2} dx (1 + \alpha x) e^{i\beta x}. \quad (B14) $$

This integral has a straightforward analytical expression in terms of trigonometric functions of $\alpha$ and $\beta$. The same steps can be followed to adapt this development to the other edges of the integration domain.

## References

* Amorisco, N. C. 2021, arXiv:2111.01148

* Benson, A. J., & Du, X. 2022, MNRAS, 517, 1398

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

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

* Delos, M. S., & White, S. D. M. 2023, JCAP, 2023, 008

* Despali, G., Vegetti, S., White, S. D. M., et al. 2022, MNRAS, 510, 2480

* Dootson, D., & Magorrian, J. 2022, arXiv:2205.15725

* Drakos, N. E., Taylor, J. E., & Benson, A. J. 2020, MNRAS, 494, 378

* Drakos, N. E., Taylor, J. E., & Benson, A. J. 2022, MNRAS, 516, 106

* Dury, V., de Rijcke, S., Debattista, V. P., & Dejonghe, H. 2008, MNRAS, 387, 2

* Dutton, A. A., Macciò, A. V., Dekel, A., et al. 2016, MNRAS, 461, 2658

* Errani, R., & Navarro, J. F. 2021, MNRAS, 505, 18

* Errani, R., Navarro, J. F., Ibata, R., & Peñarrubia, J. 2022, MNRAS, 511, 6001

* Errani, R., & Peñarrubia, J. 2020, MNRAS, 491, 4591

* Facchinetti, G., Stref, M., & Lavalle, J. 2022, arXiv:2201.09788

* Fellhauer, M., Kroupa, P., Baumgardt, H., et al. 2000, NewA, 5, 305

* Freundlich, J., Dekel, A., Jiang, F., et al. 2020, MNRAS, 491, 4523

* Fridman, A. M., Polyachenko, V. L., Aries, A. B., & Poliakoff, I. N. 1984, Physics of Gravitating Systems. II. Nonlinear Collective Processes: Nonlinear Waves, Solitons, Collisionless Shocks, Turbulence (Berlin: Springer)

* Gnedin, O. Y., Hernquist, L., & Ostriker, J. P. 1999, ApJ, 514, 109

* Gnedin, O. Y., & Ostriker, J. P. 1999, ApJ, 513, 626

* Goerdt, T., Gnedin, O. Y., Moore, B., Diemand, J., & Stadel, J. 2007, MNRAS, 375, 191

* Green, S. B., & van den Bosch, F. C. 2019, MNRAS, 490, 2091

* Hayashi, E., Navarro, J. F., Taylor, J. E., Stadel, J., & Quinn, T. 2003, ApJ, 584, 541

* Heggie, D. C., Breen, P. G., & Varri, A. L. 2020, MNRAS, 492, 6019

* Hernquist, L. 1990, ApJ, 356, 359

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

* Jalali, M. A., & Hunter, C. 2005, ApJ, 630, 804

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

* Kundic, T., & Ostriker, J. P. 1995, ApJ, 438, 702

* Li, Z., Dekel, A., Mandelker, N., Freundlich, J., & François, T. L. 2023, MNRAS, 518, 5356

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

* Lilley, E. J., & van de Ven, G. 2023, A&A, 672, A91

* Lovell, M. R. 2024, MNRAS, 529, 4050

* Murali, C. 1999, ApJ, 519, 580

* Peñarrubia, J., Benson, A. J., Walker, M. G., et al. 2010, MNRAS, 406, 1290

* Peñarrubia, J., Navarro, J. F., & McConnachie, A. W. 2008, ApJ, 673, 226

* Petersen, M. S., Roule, M., Fouvry, J.-B., Pichon, C., & Tep, K. 2024, MNRAS, 530, 4378

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

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

* Polyachenko, E. V., & Shukhman, I. G. 2015, MNRAS, 451, 601

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

* Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464

* Rozier, S., Famaey, B., Siebert, A., et al. 2022, ApJ, 933, 113

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

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

* Sellwood, J. A., & McGaugh, S. S. 2005, ApJ, 634, 70

* Stref, M., Lacroix, T., & Lavalle, J. 2019, Galax, 7, 65

* Stücker, J., Ogiya, G., Angulo, R. E., Aguirre-Santaella, A., & Sánchez-Conde, M. A. 2023a, MNRAS, 521, 4432

* Stücker, J., Ogiya, G., White, S. D. M., & Angulo, R. E. 2023b, MNRAS, 523, 1067

* Taylor, J. E., & Babul, A. 2001, ApJ, 559, 716

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

* van den Bosch, F. C., & Ogiya, G. 2018, MNRAS, 475, 4066

* Vegetti, S., & Koopmans, L. V. E. 2009, MNRAS, 400, 1583

* Vegetti, S., Koopmans, L. V. E., Bolton, A., Treu, T., & Gavazzi, R. 2010, MNRAS, 408, 1969

* Weinberg, M. D. 1994a, AJ, 108, 1398

* Weinberg, M. D. 1994b, AJ, 108, 1403

* Weinberg, M. D. 1994c, ApJ, 421, 481

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

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

* Young, P. 1980, ApJ, 242, 1232