# Predicting the linear response of self-gravitating stellar spheres and discs with LinearResponse.jl

Michael Petersen, Mathieu Roule, Jean-Baptiste Fouvry, Christophe Pichon, Kerwann Tep

► To cite this version:

Michael Petersen, Mathieu Roule, Jean-Baptiste Fouvry, Christophe Pichon, Kerwann Tep. Predicting the linear response of self-gravitating stellar spheres and discs with LinearResponse.jl. Monthly Notices of the Royal Astronomical Society, 2024, 530 (4), pp.4378-4394. 10.1093/mnras/stae732. hal-04733014

HAL Id: hal-04733014

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

Submitted on 11 Oct 2024

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>
# Predicting the linear response of self-gravitating stellar spheres and discs with LinearResponse.jl

Michael S. Petersen<sup>1,2,★</sup> Mathieu Roule,<sup>2</sup> Jean-Baptiste Fouvry,<sup>2</sup> Christophe Pichon<sup>2,3,4</sup> and Kerwann Tep<sup>2</sup>

Accepted 2024 February 21. Received 2024 February 21; in original form 2023 November 17

## ABSTRACT

We present LinearResponse.jl, an efficient, versatile public library written in JULIA to compute the linear response of self-gravitating (three-dimensional spherically symmetric) stellar spheres and (two-dimensional axisymmetric razor-thin) discs. LinearResponse.jl can scan the whole complex frequency plane, probing unstable, neutral and (weakly) damped modes. Given a potential model and a distribution function, this numerical toolbox estimates the modal frequencies as well as the shapes of individual modes. The libraries are validated against a combination of previous results for the spherical isochrone model and Mestel discs, and new simulations for the spherical Plummer model. Beyond linear response theory, the realm of applications of LinearResponse.jl also extends to the kinetic theory of self-gravitating systems through a modular interface.

**Key words:** gravitation – instabilities – software: public release – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

The response of self-gravitating systems to perturbations has been studied for more than half a century, mostly using the *matrix method* (Kalnajs 1971). As discussed in detail in section 5.3.2 of Binney & Tremaine (2008), this method uses basis functions to represent the gravitational response of the stellar system. The response matrix **M** at a given complex frequency $\omega$ is generically

$$ \mathbf{M}(\omega) = \sum_n \int_L d\mathbf{J} \frac{\mathbf{G}_n(\mathbf{J})}{n \cdot \mathbf{\Omega}(\mathbf{J}) - \omega} \quad (1) $$

This equation encodes the physics of (linear) stability in self-gravitating systems. Here, $\sum_n$ is the sum over all (allowed) resonance vectors **n**; $\int_L d\mathbf{J}$ is the scan over the populated orbital (action) space (with an appropriate Landau prescription); $n \cdot \mathbf{\Omega}(\mathbf{J}) - \omega$ is the resonant amplification at a given (complex) frequency $\omega$ and for orbital frequencies $\mathbf{\Omega}$; $\mathbf{G}_n$ are functions of the system's distribution functions (DFs) and are matrices by virtue of decomposing the Newtonian pairwise interaction on potential basis elements that encode the long-range nature of the gravitational interaction (Kalnajs 1976).

The system sustains a mode at a complex frequency $\omega$ when the condition

$$ \det[\boldsymbol{\epsilon}(\omega)] = 0, \quad \text{with } \boldsymbol{\epsilon}(\omega) = \mathbf{I} - \mathbf{M}(\omega), \quad (2) $$

is met, where **I** is the identity matrix, and $\boldsymbol{\epsilon}$ is called the (gravitational) dielectric matrix. We write frequencies satisfying this equality as $\omega_M$, which we further break down into $\omega_M = \Omega_M + i\gamma_M$, with the oscillation frequency, $\Omega_M$, and the growth rate, $\gamma_M$. Modes may appear at any point in the complex plane. Modes with $\gamma_M > 0$ are unstable, modes with $\gamma_M = 0$ are neutral, and modes with $\gamma_M < 0$ are damped.

Previous studies in linear response have probed modal structure in self-gravitating discs (e.g. Zang 1976; Toomre 1981; Vauterin & Dejonghe 1996; Pichon & Cannon 1997; Evans & Read 1998a; Jalali & Hunter 2005; Fouvry et al. 2015; De Rijcke & Voulis 2016) and spheres (e.g. Polyachenko & Shukhman 1981; Weinberg 1989; Saha 1991; Weinberg 1991; Bertin et al. 1994; Murali & Tremaine 1998; Rozier et al. 2019; Fouvry & Prunet 2022; Weinberg 2023). Most calculations have been interested in studying instabilities, i.e. modes with $\gamma_M > 0$ (for a detailed look at stability, see Palmer 1994). Traditionally, finding damped modes has been challenging because of the analytical continuation required by Landau's prescription. Zang (1976) presented a first prediction in razor-thin discs while relying heavily on the scale invariance of the Mestel potential. Weinberg (1994) broke through by using rational functions to (numerically) continue calculations made in the upper half frequency plane to the lower half, finding a lopsided damped mode in the King family of spherical models. Standard techniques in plasma physics offer other options for analytic continuation, which were applied to stellar systems by Fouvry & Prunet (2022). Damped modes are interesting for their capacity to dissipate energy through interactions of the mode and resonating stellar orbits (Nelson & Tremaine 1999). This phenomenon may have implications for real astrophysical systems,

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Society

<footnote n="1"><sup>1</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, UMR 7095, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>IPhT, DRF-INP, UMR 3680, CEA, L'orme des Merisiers, Bât 774, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="4"><sup>4</sup>Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote>* E-mail: Michael.Petersen@roe.ac.uk</footnote>
including accelerating the overall relaxation (see e.g. Hamilton et al. 2018; Heggie, Breen & Varri 2020), as well as the secular dynamics captured by the Balescu–Lenard equation (Heyvaerts 2010; Chavanis 2012), and the so-called mode-particle interactions (Hamilton & Heinemann 2020).

In this paper, building upon the method developed in Fouvry & Prunet (2022, hereafter FP22), we introduce software libraries in the JULIA language to perform linear response calculations. Details of the toolbox are given in Appendix A. Linear response calculations are intricate, requiring expensive integrals over the full action phase that must be performed with care. In the past, codes and methods for linear response have often not been publicly shared, necessitating complicated re-implementations. By publishing our code, we hope to create a framework where others can build on the methods, or add their own. By using JULIA, we are able to rapidly explore parameter space and convergence rates, while retaining a high level of readability and optional interactivity.

The paper is organized as follows. In Section 2, we sketch the linear response computations and present the JULIA libraries developed to tackle them. In Section 3, we compute the response for the spherical Plummer model, probing its modal spectrum across a range of DFs, and compare with N-body simulations. In Section 4, we validate our extension to the linear response of razor-thin discs on the well-studied constant circular velocity discs (Zang 1976). We wrap up and conclude in Section 5. Throughout the main text, technical details are kept to a minimum and deferred to appendices or to relevant references.

## 2 LINEAR RESPONSE OF SPHERES AND DISCS

The linear response of spheres and discs can be split in independent harmonics $ℓ$ (historically denoted $m$ for the discs). For these systems and within the appropriate ‘resonance’ coordinate system ($u, v$) introduced in equation (10) of FP22, each element of the response matrix of harmonic $ℓ$ reads

$$M_{pq}^{\ell}(\omega) = \sum_{n \in \mathbb{Z}^2} \int_{-1}^{1} du \frac{G_{pq}^{\ell n}(u)}{u - \omega_n(\omega)}, \quad (3)$$

where the functions<fnref n="1" /> $G_{pq}^{\ell n}(u)$ depend on the Fourier transform of bi-orthogonal basis elements (indexed by $p$ and $q$), the resonance number, $n = (n_1, n_2)$, and involve an integral over some other coordinate $dv$, as explained in Appendix A4.1. These functions, which depend on the dimensionality of the problem, are detailed in Appendix B. In equation (3), $\omega_n(\omega)$ stands for the rescaled (complex) frequency introduced in equation (11) of FP22. Importantly, the imaginary parts of the rescaled frequency $\omega_n$ and $\omega$ share the same sign. Computing the linear response of a stellar system mostly entails carrying out the integrals outlined in equation (3). It is crucial to focus on the resonant denominator, adhering to Landau’s prescription as detailed in, for instance, equation (2) in the work by FP22. The goal behind producing LinearResponse.jl is to extend the prescription used for the isochrone sphere in FP22 to any numerically given potential and DFs, as well as to discs. To be as versatile as possible, the software is decomposed in four libraries:

(i) OrbitalElements.jl (Appendix A1), which, given any central potential and its two first derivatives, performs various changes of

orbital coordinates to compute the actions ($J_r, L$), orbital frequencies ($\Omega_1, \Omega_2$) and resonance coordinates ($u, v$).<fnref n="2" />

(ii) AstroBasis.jl (Appendix A2) provides bi-orthogonal bases of potential-density pairs for spheres and discs. This library is constructed such that a user could readily supply their own additional bases using a straightforward template.

(iii) FiniteHilbertTransform.jl (Appendix A3) performs the (Landau’s prescription-compliant) finite Hilbert transform<fnref n="3" /> using Legendre polynomials, as detailed in appendix D of FP22.

(iv) LinearResponse.jl (Appendix A4) is the driver library. Relating on the three (independent) previous libraries, LinearResponse.jl implements the computation of the response matrix $M'(\omega)$. It performs the Fourier transform of the basis elements, then computes the functions $G_{pq}^{\ell n}(J)$, before finally computing the finite Hilbert transform.

We refer to the appendices for details on the implementation and to the respective GitHub repositories for installation and usage. The bulk of the following sections is devoted to exercising the library in a series of simple systems.

## 3 PLUMMER SPHERE: THEORY AND SIMULATIONS

In order to test and validate LinearResponse.jl, we first examine models representing spherically symmetric globular clusters. Following FP22, we consider the isochrone potential in Appendix C. We use its analytical expressions to check our numerical calculation of the orbital elements, and recover all of FP22’s results. This is reassuring. We also point out two caveats raised by our exploration, namely (i) diverging gradients of the anisotropic DFs at the edge of the domain leading to convergence issues in the $ℓ = 2$ case, and (ii) a (neutral) translation mode that becomes damped in the $ℓ = 1$ case.

In this section, we expand our linear response study to the Plummer (1911) model. In particular, we first study the ($ℓ = 2$) radial orbit instability (ROI) for radially anisotropic models (Section 3.1), before investigating isotropic DFs and $ℓ = 1$ modes (Section 3.2). We then simulate the Plummer model with EXP, a basis function expansion N-body code, to test our predictions (Section 3.3). We provide the relevant equations for the Plummer model in Appendix D.

### 3.1 $ℓ = 2$ modes – ROI

The ROI arises in models with some degree of radial anisotropy, i.e. the DFs contains an enhancement at low angular momentum values (for a review, see Maréchal & Perez 2011). Often, this is accomplished via the introduction of some anisotropy radius, $r_a$, outside of which the fraction of radial orbits is enhanced. We follow an Osipkov–Merritt form for the DFs (Binney & Tremaine 2008). The mechanisms for the instability are discussed in Palmer (1994) and Polyachenko & Shukhman (2015). The classic ROI results in a quadrupole ($ℓ = 2$) zero oscillation frequency mode ($\Omega_M = 0$) that grows exponentially in time ($\gamma_M > 0$). In practice, the ROI is generated by the inner Lindblad resonance (ILR), $n = (-1, 2)$, and its opposite. This simplicity makes it a particularly compelling instability to use as benchmark.

<footnote n="2">²At present, the mappings are restricted to cored potentials. Extending them to cuspy potentials will be the purpose of future work.</footnote>

<footnote n="3">³As discussed in more detail in Appendix A3, equation (A10) is the finite-interval version of the Hilbert transformation (Tricomi 1951), hence the name of this software.</footnote>

<footnote n="1">¹Note that $G_{pq}^{\ell n}(u)$ in equation (3) slightly differs from $G_n(J[u])$ in equation (1). See Appendix A4.1 for details.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left panel: Growth rate $\gamma_M = \text{Im}[\omega_M]$ of the ROI as a function of the anisotropy radius $r_a$ for the Plummer model. Right panel: Predicted potential fluctuation as a function of radius for the ROI-unstable modes in the Plummer model. Different curves correspond to different values of $r_a$, which are colour-coded. Even though the growth rate drops by three orders of magnitude, the shape of the mode changes very little. All calculations are made using the Fiducial settings.</figcaption></figure>

Various studies have sought to establish the boundary for stability in the radially anisotropic Plummer model. Dejonghe & Merritt (1988) estimated the threshold for stability in the Plummer model to be $r_a/b_c \simeq 0.9$ using the criteria of Barnes, Goodman & Hut (1986), with $b_c$ the Plummer's scale radius (equation D1). Using the mean-field code from Merritt (1987), Dejonghe & Merritt (1988) also reported that the stability threshold was somewhere around $r_a/b_c \simeq 1.1$. Relying on direct N-body simulations, Breen, Varri & Heggie (2017) confirmed that $r_a/b_c = 0.75$ was unstable, but, in contrast to Dejonghe & Merritt (1988), Breen et al. (2017) found that $r_a/b_c = 1.0$ was stable.

We test our linear response machinery on a range of $r_a/b_c$ values. Fig. 1 shows the results of searching for modes using a 'Fiducial' set of control parameters.<fnref n="4" /> We find a threshold for instability of $r_a/b_c = 1.035$, as highlighted by the vertical dashed line in Fig. 1. Let us finally note that, here, the response matrix was computed using the basis provided by Clutton-Brock (1973) tailored for the Plummer potential. Hence, the basis matches the underlying mean profile at large radii: this drastically reduces numerical noise. We report results for varying control parameters in Table C1, finding that the uncertainty in growth rate from control parameters is a few per cent. This leads to no alteration in the threshold for instability.

Since individual $ℓ$ harmonics are decoupled in linear response calculations, we investigate $ℓ = (1, 3, 4)$ for similar instability behaviour. We do not find evidence for any modes for $ℓ \neq 2$.

## 3.2 A search for $ℓ = 1$ damped modes

We next search for evidence of damped $ℓ = 1$ modes, akin to the findings of Weinberg (1994) and Heggie et al. (2020) for King models. Using an isotropic DFs (see Appendix D), we perform a search in the complex plane. We find zeros in the determinant of the dielectric matrix $e$, but these appear to be spurious non-converged

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Predicted density mode shape vs. radius for $ℓ = 1$ damped modes, using several combinations of configuration parameters. Predictions from different configurations are indistinguishable in mode shape, despite a range predictions for $\omega_M$ (see Table C2). Predictions for the Plummer model are shown in black. The dashed curve is $d\rho/dr$ for the Plummer model. The predicted mode shapes are indistinguishable from $d\rho/dr$.</figcaption></figure>

realizations of the neutral translation mode. The right column of Table C2 lists the results for different configurations. As the number of resonances considered increases, the predicted mode's frequency shifts towards the origin of the complex frequency plane, as in the isochrone model (Appendix C2). In addition, the 'false' mode clearly resembles the density derivative $d\rho/dr$, i.e. the modal shape associated with an infinitesimal translation of the cluster. As these unconverged poles are the only one we found, we do not predict any robust $ℓ = 1$ modes. For completeness, we also search other harmonics $ℓ = (2, 3, 4)$, finding no evidence for damped modes.

Fig. 2 shows the predicted mode shape with radius for the Plummer model (black curve), as compared to $d\rho/dr$. As expected given the results of FP22, the predicted mode appears to be the neutral $\omega_M = 0 + 0i$ mode, pushed slightly away from the origin of the complex frequency plane due to incomplete convergence. We show the mode shape for the full set of configurations (Fig. 2), but all curves lie on top of one another. Surprisingly, even though the frequency of the mode is not converged, its shape appears to be. Here, we recover again that reconstructing a large-scale translation mode from orbital space

<footnote n="4">⁴The 'Fiducial' control parameters for spherical calculations take the software defaults for OrbitalElements.jl (Table A1) and LinearResponse.jl (Table A3), along with a basis scale $r_b = 5.0$ and 100 radial basis functions. For $ℓ = 2$ calculations, this results in 62 resonances to calculate. For $ℓ = 1$ calculations, this results in 42 resonances to calculate. For the Plummer sphere, the fiducial value of $\Omega_0$ is 2, cf. equation D2.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left panel: Growth of the $\ell = 2$ potential fluctuations energy ($\tilde{E}_\ell$, equation E3), through time in N-body simulations of the radially-anisotropic Plummer sphere for four different anisotropy radii, $r_a/b_c$. Simulations with $r_a/b_c = 0.95$ and $r_a/b_c = 1.00$ are not shown, as instabilities were not reliably measured. Thin curves are individual simulation realizations; thick curves are the ensemble average for each value of $r_a/b_c$. Curves are aligned on the time of the first peak of the potential perturbation, denoted $t = 0$, where the shape of each mode is measured. The exponential growth predicted by linear theory is clearly visible before the modes' saturation. Right panel: Shape of the $\ell = 2$ potential fluctuations measured from the N-body simulations (solid curves). The dashed curves are the corresponding predictions from LinearResponse.jl. The linear predictions and numerical measurements are in pleasant agreement.</figcaption></figure>

integrals and resonant interactions is, by design, a challenging task.

## 3.3 N-body simulation

In this section, we report the results of N-body simulations of the Plummer model. We test two different DFs (see Appendix D): one radially anisotropic, and one isotropic. We realise the N-body simulations using initial conditions from PlummerPlus (Breen et al. 2017). We run the models using EXP (Weinberg 1999; Petersen, Weinberg & Katz 2022), a basis function expansion N-body code. Further numerical details for both the initial conditions and the N-body integration are given in Appendix E.

### 3.3.1 Anisotropic simulations, ROI

To attempt to recover the ROI in the Plummer clusters, we select six different values of $r_a/b_c$ for which to run N-body simulations: [0.75, 0.80, 0.85, 0.90, 0.95, 1.0]. These values span the entire range of models predicted to be unstable. For each $r_a/b_c$, we run an ensemble of six simulations of $2^{18}$ particles each.

Fig. 3 shows the results for the $r_a/b_c = [0.75, 0.80, 0.85, 0.9]$ runs using EXP (see Appendix E). In general, the runs for any value of $r_a/b_c$ show similar behaviour, though the dispersion between realizations increases with $r_a/b_c$ (i.e. dispersion between realizations increases as the model approaches marginal stability). The solid curves in the right panel of Fig. 3 show the reconstruction of the potential fluctuations for the series of $r_a/b_c$ models; the qualitative agreement with the predicted mode shapes and trends (shown as dashed curves) with $r_a/b_c$ are remarkably high, though the N-body peak is consistently 10 per cent lower than the predicted peak. Measurements of the mode shape during the growth are consistent with the same shape. In the runs with $r_a/b_c = 0.95$, only three of the six runs performed appear to show any instability. Thus, $r_a/b_c = 0.95$ appears to be a weakly unstable model. We do not measure an instability for any of the $r_a/b_c = 1.0$ runs, in N-body simulations.

Table 1 compares the measured growth rates from the N-body simulations to the predictions from LinearResponse.jl. The agreement is quantitatively good, in that all measured growth rates are

<table><caption>**Table 1.** Comparison between the measured and predicted growth rates for the ROI in the Plummer sphere for various normalized anisotropy radii, $r_a/b_c$. The scale frequency, $\Omega_0$, is defined for the Plummer cluster (equation D2). Uncertainties on the N-body runs arise from variations among different realizations (i.e. initial conditions). Predictions are for the Fiducial configuration, with uncertainties measured from different configurations (Table C1). Measurements and predictions satisfyingly agree and retrieve that the smaller the anisotropy radius, the stronger the instability.</caption><thead><tr><th>r<sub>a</sub>/b<sub>c</sub></th><th>N-body (γ)/Ω<sub>0</sub></th><th>Linear response (γ)/Ω<sub>0</sub></th></tr></thead><tbody><tr><td>0.75</td><td>0.046 ± 0.002</td><td>0.048 ± 0.001</td></tr><tr><td>0.80</td><td>0.038 ± 0.002</td><td>0.038 ± 0.001</td></tr><tr><td>0.85</td><td>0.033 ± 0.003</td><td>0.029 ± 0.001</td></tr><tr><td>0.90</td><td>0.023 ± 0.003</td><td>0.020 ± 0.001</td></tr></tbody></table>

within $3\sigma$ of the predicted growth rate, and both approaches recover the qualitative trend of decreasing $\gamma_M$ with increasing $r_a/b_c$. We do not find evidence for a second unstable mode in the simulations. The numerical values from the Fiducial run, except with a basis scale length of $r_b = 1$, match the basis used to run the simulation parameters (Appendix E). This ensures a straightforward comparison with the N-body simulations. We find that the predicted eigenvectors (equation A20) qualitatively match the basis function coefficients in the N-body simulations.

In conclusion, the N-body simulations and the predicted linear response growth rate and mode shape match to a high degree, validating both our linear response approach, and the ability of EXP to make intricate dynamical measurements. Future work can increase the number of particles in the runs to validate the predicted growth rates, and attempt to recover lower growth rates near the stability boundary.

### 3.3.2 Isotropic simulations, $\ell = 1$ measurements

By analogy with the investigation of Section 3.2, we perform an ensemble of six simulations of Plummer spheres with an isotropic DFs (see Appendix E). The left panel of Fig. 4 shows the run with time of the (rescaled) energy, $\tilde{E}_1 \propto E_{\ell=1}$ (equation E3), for each simulation. We do not find any evidence for an $\ell = 1$ mode in
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left panel: Amplitude of the $\ell = 1$ potential fluctuations energy ($\tilde{E}_\ell$, equation E3) versus time in the isotropic N-body realizations of Plummer spheres. Surprisingly, this time evolution could be compatible with a wave kinetic equation, as in equation (4). The best-fitting solution is shown as a dashed grey curve. Right panel: Radial shape of $\ell = 1$ density fluctuations measured from the N-body simulations at three different times (denoted by vertical coloured dashed lines in the left panel). The dashed grey curve in the right panel shows $d\rho/dr$ versus radius for the Plummer cluster. As time increases, the measured dipolar perturbations get more and more similar to a simple translation of the cluster as a whole.</figcaption></figure>

isotropic Plummer models; the run of $\tilde{E}_1$ with time does not exhibit an exponential increase. Instead, $\tilde{E}_1$ increases and slowly approaches a saturated level that differs between realizations. Upon inspection, the amplitude $\tilde{E}_1$ is dominated by the lowest order function in each realization, i.e. these fluctuations are large scale. None of the basis elements show clear evidence of a coherent oscillation frequency.

The right panel of Fig. 4 may be compared with Fig. 2 for the predictions of modal shapes. While the linear predictions are consistent with the only apparent mode being the neutral mode with shape $d\rho/dr$, the results from the N-body simulation are generally not consistent with $d\rho/dr$. Instead, the peak density deviation is confined to smaller radii. Over time in the simulation, the lowest order function dominates further, such that the shape of the $\ell = 1$ density fluctuations further resembles $d\rho/dr$. However, over the duration of the simulations here, the cluster does not fully offset. Rather, the inner regions of the cluster offset, while the outer regions appear to remain roughly fixed. This likely owes to the long dynamical times in the cluster’s outskirts.

Interestingly, in Fig. 4, the time evolution of the total energy in the $\ell = 1$ fluctuations shows a clear sign of saturation. Rogister & Oberman (1968, equation 21 therein) put forward a wave kinetic equation (see also equation 22 in Hamilton & Heinemann 2020) that describes the (linear) saturation of potential fluctuations in self-gravitating systems. It reads

$$ \frac{d\tilde{E}_\ell}{dt} = 2\gamma_M (\tilde{E}_\ell - \tilde{E}_{th}), \quad (4) $$

where $\tilde{E}_\ell$ is the (rescaled) energy in the $\ell = 1$ gravitational fluctuations and $\tilde{E}_{th}$ the (rescaled) asymptotic level of dressed fluctuations. Equation (4) states that in a self-gravitating system which sustains a (weakly) damped mode with damping frequency $\gamma_M$, transients only fade away after a (few) $1/\gamma_M$, at which stage fluctuations are fully dressed by self-gravity (see e.g. Hamilton & Heinemann 2023). Solutions of equation (4) are readily obtained. Given the run of $\tilde{E}_{\ell=1}$ in Fig. 4, we can use the ensemble average $\ell = 1$ energy (the thick black curve in the left panel of Fig. 4) to find best-fitting values for $\tilde{E}_{th}$, $\gamma_M$ and the initial energy fluctuations, $\tilde{E}_{init} = \tilde{E}_\ell(t=0)$. We find $\tilde{E}_{init} = 1.0 \times 10^{-6}$, $\tilde{E}_{th} = 8.0 \times 10^{-3}$, and $\gamma_M/\Omega_0 = -0.001$. The time evolution of the numerical simulations

appears to be reasonably compatible with the wave kinetic equation in the isotropic case. However, the fitted solution for $\gamma_M$ is outside of the region that is numerically accessible with the current methods used in LinearResponse.jl, despite implying a relatively long time for the $\ell = 1$ energy to reach a thermal state. Additionally, the apparent scatter of values of $\tilde{E}_{th}$ for individual runs is somewhat unexpected. We defer further investigation to future work.

## 4 MESTEL RAZOR-THIN DISCS: VALIDATION

To emphasize the generality of LinearResponse.jl, we now consider razor-thin [two-dimensional (2D)] discs. In that case, stars still interact through the usual Newtonian interaction potential $U(r, r') = -G|r - r'|$, with G the gravitational constant, but are confined to a plane. Assuming an axisymmetric mean DFs, the stars’ mean-field orbits can be described by the two action variables ($J_r, L$), just like in the case of spheres. Hence, all the previous numerical methods from FP22 naturally apply to razor-thin discs. The only required modifications are (i) marginally modifying the linear response integrand, $G^{\text{in}}_{pq}$, in equation (3), given in Appendix B2, and (ii) implementing 2D bi-orthogonal basis elements (e.g. Clutton-Brock 1972; Kalnajs 1976), given in Appendix F1.

In order to validate our implementation, we aim to recover well-documented unstable modes in razor-thin discs. In practice, following the work of Zang (1976), we consider Mestel discs (Mestel 1963), whose scale-invariance allows for various analytical simplifications. None the less, in what follows, we do not use these simplifications and rather use the generic scheme from LinearResponse.jl. Following Evans & Read (1998a), the DFs of the stars is tapered in the central region. The presence of an inner cut-off mimicks an unresponsive central bulge, hence introducing a reflexive boundary<fnref n="5" />: the sharper this inner cut-out, the stronger the instability (Zang 1976). The disc’s outskirts are also tapered, though this does not impact the disc’s stability, provided that this external cut-out is sharp enough

<footnote n="5">⁵The Mestel potential diverges for $r \to 0$, making exactly radial orbits singular. As a result, we also softened the mean gravitational potential in the central region (Appendix F2). This softening is on a much smaller scale than the DFs’s inner cut-out, and does not impact the linear response predictions.</footnote>
<table><caption>**Table 2.** Complex frequency, $\omega_M/\Omega_0$, of the most unstable $\ell = 2$ mode of truncated Zang discs ($q = 6$), as one varies the index of the inner taper, $N$, with $\Omega_0$ the disc's scale frequency (see Appendix F2). Here, we compare the numerical prediction from *LinearResponse.jl* with the values from Evans & Read (1998b). We refer to Appendix F2 for detailed parameters. Both predictions satisfactorily agree within 1 per cent precision. The results prove robust to variations of the control parameters (Table F2 for the $N = 4$ disc).</caption><thead><tr><th>Cut-in index <i>N</i></th><th>This work</th><th>Evans & Read (1998b)</th></tr></thead><tbody><tr><td>4</td><td>0.878 + 0.126i</td><td>0.879 + 0.127i</td></tr><tr><td>6</td><td>0.901 + 0.221i</td><td>0.902 + 0.222i</td></tr><tr><td>8</td><td>0.922 + 0.265i</td><td>0.922 + 0.266i</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Isocontours of the determinant of the $\ell = 2$ (gravitational) dielectric matrix from equation (2) for the $N = 4$ Zang disc. The dominant mode obtained by Zang (1976) is highlighted with a yellow cross and is recovered within 1 per cent precision.</figcaption></figure>

and far enough (Evans & Read 1998b). We refer to Appendix F2 for details and notably the definition of the scaling frequency $\Omega_0$. Note that the method based on finite Hilbert transform implemented in *LinearResponse.jl* is not perfectly suited here. Indeed, because of its central divergence, the frequency range in Mestel discs is not finite. We deal with this particular issue by limiting the orbital domain probed by the code, constrained by the inner DFs taper (see Appendix A1.3 for details). It would be rewarding to extend FP22's methodology to the case of semi-finite frequency supports.

Once this model is set up, we perform stability analysis for two-armed modes, i.e. $\ell = 2$ modes, as one varies the properties of the inner taper.<fnref n="6" /> As reported in Table 2, we find a satisfying agreement between the (semi-analytical) predictions of Zang (1976) and Evans & Read (1998b), and the present linear predictions for the growth rate and oscillation frequency of the most unstable mode. These predictions have already been confirmed using numerical simulations by Sellwood & Evans (2001).

In Fig. 5, we present a typical map of the complex frequency plane one can obtain from *LinearResponse.jl*. More precisely, Fig. 5 illustrates the determinant of the (gravitational) dielectric matrix, $\boldsymbol{\varepsilon}(\omega)$ from equation (2), through its level contours in the complex-frequency plane. This determinant vanishes at the frequency $\omega_M/\Omega_0 = 0.878 + 0.126i$, i.e. the system supports a growing mode. The shape of this unstable mode is reported in Fig. 6, where we compare it with the result from Zang (1976). We find a quantitative match between both approaches. In Fig. 5, the saturation and ringing for damped frequencies are to be expected. They are a direct consequence of analytical continuation being an ill-conditioned numerical problem. Further explanation is given in Appendix A3.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Shape of the $\ell = 2$ dominant mode of the $N = 4$ Zang disc as predicted by *LinearResponse.jl* (coloured dashed lines) overlaid with the shape obtained by Zang (1976) (in black, fig. 9 therein). For both shapes, the contours denote the 10, 20, 40, 60, and 80 per cent of the peak density, with only the overdensity being represented. The dotted circles show the corotation (CR) and ILR radii of the mode. Both linear predictions are in very good agreement.</figcaption></figure>

As a concluding remark, let us note that here we used the generic basis from Clutton-Brock (1972), which is not tailored to asymptotically match the disc's underlying potential. Interestingly, this did not impact our ability to recover precisely the underlying modes. Expanding *LinearResponse.jl* to accommodate for more generic basis elements will nonetheless be the topic of future work.

## 5 SUMMARY AND CONCLUSIONS

In this paper, we presented *LinearResponse.jl*, an efficient, readable software written in JULIA to perform linear response calculations for self-gravitating stellar spheres and discs. Given a model's potential and DFs, it can predict the system's modal response, as well as the shapes of individual modes. We validate the four underlying libraries against known results for the isochrone sphere, before investigating the Plummer model.

In the case of the isochrone model, we find quantitative agreement with previous work. Studying the $\ell = 1$ response, we confirm that the results of FP22 correspond to an unconverged version of the neutral translation mode. We demonstrate that higher numerical fidelity reduces the drift of the neutral mode away from the origin of the complex frequency plane. Such explorations are straightforward with *LinearResponse.jl*.

In the case of the Plummer model, we perform our own comparison with $N$-body simulations performed using EXP. For the ROI, we find quantitative agreement between predictions and measurements in simulations. This agreement opens avenues for validating $N$-body simulations with linear response theory, as well as for studying instabilities in systems where linear response may be complicated. In the case of an isotropic DFs, we find the same unconverged neutral mode as in the isochrone model. We find no evidence for other modes. Finally, the $N$-body simulations of the isotropic model suggest that neutral translation modes may not be easily excited in a cluster's lifetime.

The extension of FP22's work to the stability analysis of razor-thin discs is convincingly validated against known semi-analytical results in Mestel discs (Zang 1976; Evans & Read 1998b). In particular, we satisfyingly retrieve both the frequency and shape of the $\ell = 2$

<footnote n="6">⁶The azimuthal harmonic number for discs is historically denoted $m$. Adapting here from the spherical case, we none the less denote it $\ell$.</footnote>
growing mode from Zang (1976). Accurately computing the linear response of these dynamically cold systems is key to appropriately studying their long-term relaxation.

All these promising results pave the way for a more systematic and rigorous exploration of linear response for a wider class of clusters and discs, including flattened or indeed triaxial systems. For example, radial orbit instabilities have long been suspected as being important for constraining dynamical models of elliptical galaxies. With novel data products such as integral field unit spectroscopy (e.g. the MaNGA survey described in Bundy et al. 2015), new constraints may be available for the orbital content of elliptical galaxies, using modal analysis with *LinearResponse.jl*.

The calculations performed in this paper may straightforwardly be extended to various other cored models and DFs: this is one of the main gains from this generic toolbox. Similarly, the software may also be used for other generic computations, such as computing orbital elements in various potentials. We hope that these libraries will serve as a base for further calculations and extensions, with the means to approach the speed of fully compiled codes and the readability of scripting languages.

Ultimately, one could eventually build the following other extensions of the software: (i) extending the mapping to cuspy profiles for which the range of frequencies is unbounded; (ii) generalizing the mapping to integrable rotating spherical clusters or flattened ones; (iii) accounting for the contributions from branch cuts; (iv) implementing analytic continuation via rational functions, following Weinberg (1994) (see appendix F in FP22, for a review); (v) building the self-gravitating amplification kernel in the time domain rather than the frequency domain (see e.g. Dootson & Magorrian 2022; Rozier et al. 2022, and references therein).

As for long-term relaxation, *LinearResponse.jl* may also be used on various fronts, for example to estimate torques that external perturbations apply on stellar systems. Indeed, the celebrated LBK formula (Lynden-Bell & Kalnajs 1972) integrates over resonance conditions. Phrased differently, it only requires one to account for the resonance part of Landau's prescription for neutral frequencies: the calculation precisely implemented here. Finally, *LinearResponse.jl* can also easily be adapted to compute the secular response of discs and spheres, through the drift and diffusion coefficients the Balescu-Lenard equation (see e.g. Hamilton et al. 2018), or to explore their adiabatic response as one slowly varies the mean field or the stellar cluster's DFs (see e.g. Reddish et al. 2022).

## ACKNOWLEDGEMENTS

We thank S. Rozier, M. Weinberg, C. Hamilton, S. de Rijcke, J. Peñarrubia, and A. Naik for insightful conversations. We thank the referee, James Binney, for constructive comments that strengthened the work. We also thank the Kavli Institute for Theoretical Physics for hosting the workshop 'Connecting Galaxies to Cosmology at High and Low Redshift' and the Higgs Centre for Theoretical Physics for hosting the workshop 'Secular evolution of self-gravitating systems'. This work is partially supported by the French Agence Nationale de la Recherche under Grant SEGAL ANR-19-CE31-0017 and by the National Science Foundation under grant no. NSF PHY-1748958. MSP acknowledges support from a UKRI Stephen Hawking Fellowship. We thank Stéphane Rouberol for the smooth running of the Infinity cluster, where the majority of the simulations were performed. This project is primarily built in the JULIA language (Bezanson et al. 2017). We used SCIPY (Virtanen et al. 2020) to fit the wave kinetic equation solution. For visualizations, we used

IPYTHON (Pérez & Granger 2007), NUMPY (Harris et al. 2020), and MATPLOTLIB (Hunter 2007).

## DATA AVAILABILITY

The JULIA libraries developed for this paper, *LinearResponse.jl*, *OrbitalElements.jl*, *AstroBasis.jl* and *FiniteHilbertTransform.jl*, are freely available on GitHub. The EXP simulations used for comparison are available upon reasonable request to the corresponding author.

## REFERENCES

* Barnes J., Goodman J., Hut P., 1986, ApJ, 300, 112

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

* Bezanson J., Edelman A., Karpinski S., Shah V. B., 2017, SIAM Review, 59, 65

* Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, NJ

* Bovy J., 2015, ApJS, 216, 29

* Breen P. G., Varri A. L., Heggie D. C., 2017, MNRAS, 471, 2778

* Bundy K. et al., 2015, ApJ, 798, 7

* Chavanis P.-H., 2012, Physica A, 391, 3680

* Clutton-Brock M., 1972, Ap&SS, 16, 101

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

* De Rijcke S., Voulis I., 2016, MNRAS, 456, 2024

* Dejonghe H., Merritt D., 1988, ApJ, 328, 93

* Dootson D., Magorrian J., 2022, preprint (arXiv:2205.15725)

* Evans N. W., Read J. C. A., 1998a, MNRAS, 300, 83

* Evans N. W., Read J. C. A., 1998b, MNRAS, 300, 106

* Fouvry J.-B., Prunet S., 2022, MNRAS, 509, 2443 (FP22)

* Fouvry J.-B., Pichon C., Magorrian J., Chavanis P.-H., 2015, A&A, 584, A129

* Fouvry J.-B., Hamilton C., Rozier S., Pichon C., 2021, MNRAS, 508, 2210

* Fridman A. M., Poliachenko V. L., 1984, Physics of Gravitating Systems. I. Springer, New York

* Hamilton C., Heinemann T., 2020, preprint (arXiv:2011.14812)

* Hamilton C., Heinemann T., 2023, MNRAS, 525, 4161

* Hamilton C., Fouvry J.-B., Binney J., Pichon C., 2018, MNRAS, 481, 2041

* Harris C. R. et al., 2020, Nature, 585, 357

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

* Hénon M., 1959, Ann. d'Astrophys., 22, 126

* Hénon M., 1971, Ap&SS, 14, 151

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

* Heyvaerts J., 2010, MNRAS, 407, 355

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

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

* Kalnajs A. J., 1971, ApJ, 166, 275

* Kalnajs A. J., 1976, ApJ, 205, 745

* King I. R., 1966, AJ, 71, 64

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

* Lynden-Bell D., 2015, MNRAS, 447, 1962

* Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1

* Maréchal L., Perez J., 2011, Transport Theor. Stat., 40, 425

* May A., Binney J., 1986, MNRAS, 221, 13

* Merritt D., 1985, AJ, 90, 1027

* Merritt D., 1987, ApJ, 319, 55

* Merritt D., 1988, AJ, 95, 496

* Mestel L., 1963, MNRAS, 126, 553

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

* Murali C., Tremaine S., 1998, MNRAS, 296, 749

* Nelson R. W., Tremaine S., 1999, MNRAS, 306, 1

* Osipkov L. P., 1979, Soviet Ast. Lett., 5, 42

* Palmer P. L., 1994, Stability of Collisionless Stellar Systems. Kluwer Academics Publ.Amsterdam

* Pérez F., Granger B. E., 2007, Comput. Sci. Eng., 9, 21

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

* Pichon C., Cannon R. C., 1997, MNRAS, 291, 616
Plummer H. C., 1911, MNRAS, 71, 460

Polyachenko V. L., Shukhman I. G., 1981, Soviet Ast., 25, 533

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

Press W. et al., 2007, Numerical Recipes, 3rd edn. Cambridge Univ. Press, Cambridge

Price-Whelan A. M., 2017, JOSS, 2, 388

Reddish J. et al., 2022, MNRAS, 512, 160

Rogister A., Oberman C., 1968, J. Plasma Phys., 2, 33

Rozier S., Fouvry J.-B., Breen P. G., Varri A. L., Pichon C., Heggie D. C., 2019, MNRAS, 487, 711

Rozier S., Famaey B., Siebert A., Monari G., Pichon C., Ibata R., 2022, ApJ, 933, 113

Saha P., 1991, MNRAS, 248, 494

Sellwood J. A., Evans N. W., 2001, ApJ, 546, 176

Toomre A., 1981, Structure and Evolution of Normal Galaxies. Cambridge Univ. Press, Cambridge, p. 111

Trefethen L. N., 2019, Approximation Theory and Approximation Practice. SIAM, Philadelphia, PA

Trefethen L. N., 2020, BIT Numer. Math., 60, 901

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

Tricomi F. G., 1951, Q. J. Math., 2, 199

Vasiliev E., 2019, MNRAS, 482, 1525

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

Virtanen P. et al., 2020, Nature Methods, 17, 261

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

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

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

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

Weinberg M. D., 2023, MNRAS, 525, 4962

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

## APPENDIX A: LINEAR RESPONSE SOFTWARE

This Appendix presents four distinct libraries, each of which controls a different unique aspect of our linear response calculation. The libraries are organized under a GitHub organization, JuliaStellarDynamics. The libraries are written in pure JULIA, a language that combines the speed of a compiled language and the readability of a scripting language. The ease of reading the code makes the libraries highly extendable without a significant barrier to entry, while the optimized compiler passes produce low-level code that can be efficiently evaluated. As many readers may be new to JULIA, we have created documentation and example notebooks for each library, available on GitHub. We refer the interested reader to the relevant documentation for direct examples of the software. In particular, a webpage for the GitHub organization lays out the linear response project fully: juliastellardynamics.github.io.

### A1 OrbitalElements.jl: orbit frequency computation

The JULIA library OrbitalElements.jl is optimized to calculate high-precision quantities for orbits in a static central potential, $\psi = \psi(r)$. For convenience, we do not use more generic libraries GALPY (Bovy 2015), GALA (Price-Whelan 2017), or AGAMA (Vasiliev 2019), but rather design our own specific library tailored for our needs. At the heart of OrbitalElements.jl is the ability, given a central potential and its two first derivatives, to convert nearly seamlessly between different orbital elements, i.e. different constants of motion, namely

(i) the pericentre and apocentre radii ($r_{per}$, $r_{apo}$);

(ii) the effective semimajor axis and eccentricity

$$a = \frac{r_{per} + r_{apo}}{2}, \quad e = \frac{r_{apo} - r_{per}}{r_{apo} + r_{per}}; \qquad (A1)$$

(iii) the energy and angular momentum (equation 1.3 in Lynden-Bell 2015)

$$E = \frac{r_{apo}^2 \psi(r_{apo}) - r_{per}^2 \psi(r_{per})}{r_{apo}^2 - r_{per}^2}, \quad L = \left[ \frac{2(\psi(r_{apo}) - \psi(r_{per}))}{r_{per}^2 - r_{apo}^2} \right]^{1/2}; \qquad (A2)$$

(iv) the actions ($J_r$, $L$), where

$$J_r = \frac{1}{\pi} \int_{r_{per}}^{r_{apo}} dr \, v_r, \qquad (A3)$$

with $J_r$ the radial action and $v_r = \sqrt{2(E - \psi(r)) - L^2/r^2}$ the instantaneous radial velocity;

(v) the orbital frequencies ($\Omega_1$, $\Omega_2$) associated respectively with the radial and azimuthal oscillations;

(vi) the frequency ratios ($\alpha$, $\beta$) from which the frequencies are computed (see equations 28 and 29 in Tremaine & Weinberg 1984)

$$\frac{1}{\alpha} = \frac{\Omega_0}{\Omega_1} = \frac{1}{\pi} \int_{r_{per}}^{r_{apo}} \frac{dr}{v_r}, \quad \beta = \frac{\Omega_2}{\Omega_1} = \frac{L}{\pi} \int_{r_{per}}^{r_{apo}} \frac{dr}{r^2 v_r}, \qquad (A4)$$

with $\Omega_0$ some given frequency scale (typically the central radial frequency in cored profiles);

(vii) the resonance-specific (i.e. dependent on $n$) coordinates ($u$, $v$) from FP22 (appendix B therein).

In practice, OrbitalElements.jl is centred around the effective semimajor axis and eccentricity ($a$, $e$), but straightforward conversions between different orbital labels exist as simple function calls. These changes of coordinates are a requirement for the linear response of self-gravitating systems. Indeed, by construction, it involves scanning the full orbital space and dealing appropriately with resonant denominators, as visible in equation (1). In FP22, these conversions were performed analytically for the isochrone model. With this library, we provide a generic computation of orbital elements for any central potential.<fnref n="7" />

#### A1.1 Forward mappings

The direct calculations compute orbital elements (frequency ratios, radial action, etc.) through integrals of the form $\int_{r_{per}}^{r_{apo}} dr...$. However, due to the diverging integrand $1/v_r$ in equations (A4), the radius $r$ is not an appropriate integration variable. To improve numerical accuracy, we cure divergences using a mapping from equation 51 in Hénon (1971). It reads

$$r(w) = a[1 + ef(w)]; \quad f(w) = w(\frac{3}{2} - \frac{1}{2}w^2), \qquad (A5)$$

with $w \in [-1, 1]$ the 'Hénon anomaly', and $w = -1$ (resp. $w = 1$) corresponding to the pericentre (resp. apocentre). Using this variable, the integrand

$$\Theta(w) = \frac{dr/dw}{v_r(w)}, \qquad (A6)$$

is no longer divergent and equations (A4) simply read

$$\frac{1}{\alpha} = \frac{1}{\pi} \int_{-1}^{1} dw \Theta(w), \qquad (7a)$$

$$\beta = \frac{L}{\pi} \int_{-1}^{1} dw \frac{\Theta(w)}{r^2(w)}. \qquad (7b)$$

<footnote n="7">⁷In practice, the library is currently limited to cored potentials. Its extension to cuspy potentials is left for future work.</footnote>
In practice, the evaluation of $\Theta$ can become numerically unstable close to the boundaries $w = \pm 1$. To cure this, we perform a second-order expansion for values of $w$ closer than some parameter EDGE. Values at the border are obtained from straightforward expansions (see e.g. equations 53–54 in Hénon 1971). The radial action, $J_r$, is readily and easily computed by integrating equation (A3) over $w$. In Fig. A1, we illustrate the function $w \mapsto \Theta(w)$ for the isochrone potential. For most orbits, $\Theta(w)$ is a smooth function that can be integrated using low-order schemes. We use the Simpson’s 1/3 rule with a parameter NINT to control the number of integration nodes.<fnref n="8" />

For near-circular orbits, i.e. $e \to 0$, and/or very small semimajor axes, i.e. $a \to 0$, the evaluation of $\Theta(w)$ can break down. To tackle this issue, we use analytic expressions for circular ($e=0$) and central ($a=0$) orbits, and then perform a second-order expansion for values close to the borders. For the sake of numerical stability, similar expansions are performed close to radial orbits, i.e. $e \to 1$. In that case, the values on the border are computed through usual integrations (there exist no analytical expressions) before being interpolated. The same interpolations are used to compute the orbit’s energy and actions. In practice, the regions of expansions are set by the parameters TOLECC and TOLA (rescaled by the characteristic radius parameter $rc$).

## A1.2 Backward mappings

We now have at our disposal forward transformations from semimajor axis and eccentricity to energy, actions and frequencies (ratios). These mappings are not analytically invertible. We therefore employ a Newton–Raphson descent to invert them, e.g. to construct the function $(\alpha, \beta) \to (a, e)$. It requires the knowledge of the forward mapping Jacobians, i.e. sets of derivatives. These derivatives are estimated via simple two-point finite differences, on the scales $\delta a$ and $\delta e$. Finally, the Newton–Raphson algorithm is controlled with a maximal number of iterations (ITERMAX) and a accuracy goal (inv$e$). The effective accuracy of the inversion is mainly determined by the orbital integration (NINT) and the finite differences ($\delta a$ and $\delta e$).

## A1.3 Resonances variables

Let us now discuss the computation of resonance-specific (i.e. dependent on $n$) coordinates $(u, v)$. Within these coordinates, resonances lines $(\omega_n = n \cdot \Omega = \text{cst.})$ become straight lines $(u = \text{cst.})$. Their construction from the frequency ratios, $(\alpha, \beta)$, is concisely presented in appendix B of FP22. We refer to this paper for definitions and notation.

In OrbitalElements.jl, we extend the resonance-specific mapping to situations where the potential and DFs are no fully self-consistent. In this situation, a significant portion of the frequency space might be unimportant for the system’s linear response. This is for example the case in tapered Zang discs (Appendix F2), where the frequency profile diverges in the ‘unpopulated’ centre.

We therefore truncate the domain in $(\alpha, \beta)$ effectively explored. This does not affect the computation of the frequencies, $\Omega$, but only the definition of the resonance variables $(u, v)$. This effective truncation is set by the two parameters $rmin$ and $rmax$. Then, the effective domain in $(\alpha, \beta)$ is restricted to the region between $\alpha_{\min} = \alpha_c(rmax)$ and $\alpha_{\max} = \alpha_c(rmin)$ with $\alpha_c(r) = \alpha(a = r, e = 0)$ the (outward decreasing) circular frequency ratio. For example, in

the case of a Mestel disc (Appendix F2), $rmin > 0$ ensures that the frequency support is finite, and that the $(u, v)$ domain is focused on orbital regions effectively populated.

Fortunately, compared to FP22, adding these truncations only amounts to slightly modifying the boundary values of the resonance frequency $\omega_n \in [\omega_n^{\min}, \omega_n^{\max}]$ and the variable $v \in [v_n^-, v_n^+]$. More precisely, the extrema of $\omega_n$ (see equation B7 in FP22) are now reached either in the four edges $(\alpha_{\min}, \frac{1}{2})$, $(\alpha_{\min}, \beta_c[\alpha_{\min}])$, $(\alpha_{\max}, \frac{1}{2})$, $(\alpha_{\max}, \beta_c[\alpha_{\max}])$, or along the curve $(\alpha, \beta_c[\alpha])$ with $\alpha_{\min} \le \alpha \le \alpha_{\max}$. For the $v$ variable, the first two constraints of equation (B10) in FP22 become $\alpha_{\min} \le v; v \le \alpha_{\max}$. Finally, the same procedure as in FP22 is used to obtain the boundary values.

For a fully self-gravitating system (e.g. a Plummer sphere as in Appendix D), one can stick to the default parameters values which do not introduce such domain restriction.

## A1.4 Parameters

For the user that wants straightforward simplicity, control parameters are largely hidden, and set to tested defaults. None the less, we also provide a documented interface to change them. Indeed, all exposed function calls allow the user to modify control parameters, through an optional final argument. The control parameters are summarized in Table A1, but are described in more detail throughout the previous subsections.

## A1.5 Tests

Figs A2 and A3 demonstrate the fidelity of our orbit frequency calculation and subsequent inversion. We use the default control parameters in OrbitalElements.jl for this accuracy benchmark.

First, in Fig. A2, we test the accuracy of $\alpha$ and $\beta$ calculations for the isochrone model (for which we know the analytic values of $\alpha$ and $\beta$). As a summary statistics of the accuracy, we define the distance between the numerically calculated values of $\alpha$ and $\beta$ (denoted $\tilde{\alpha}, \tilde{\beta}$) and the analytic values (simply denoted $\alpha, \beta$) through

$$ \Delta_{\alpha\beta} = \sqrt{(|\tilde{\alpha} - \alpha|/\alpha)^2 + (|\tilde{\beta} - \beta|/\beta)^2}. \quad (A8) $$

In practice, $\Delta_{\alpha\beta}$ is dominated by the accuracy contribution from $\tilde{\beta}$. To probe the full range of relevant semimajor axes, we sample $a$ using $1000$ log-spaced points $a/b_c \in [10^{-3}, 10^3]$, and $1000$ linear-spaced points $e \in [0, 1]$. The accuracy is better than 1 per cent at all points sampled, and in most cases is better than 0.001 per cent accurate.

In Fig. A3, we use the same $1000 \times 1000$ grid of $(a, e)$ to test the inversion from frequencies back to $(a, e)$. This time we use the Plummer potential to test both the numerical computation of $(a, e) \mapsto (\tilde{\alpha}, \tilde{\beta})$ and the subsequent numerical inversion $(\tilde{\alpha}, \tilde{\beta}) \mapsto (\tilde{a}, \tilde{e})$. The figure shows a combination of the relative difference $|\tilde{a} - a|/a$ and absolute difference $|\tilde{e} - e|$ for each point in the grid, given by

$$ \Delta_{ae} = \sqrt{(|\tilde{a} - a|/a)^2 + (|\tilde{e} - e|)^2}. \quad (A9) $$

In practice, we find that the inversions are quite accurate, with typical differences on the order of 0.001 per cent. However, in certain regions it becomes difficult to accurately recover the input $(a, e)$ values, namely orbits close to radial or circular, and orbits with small semimajor axes. Fortunately, these inaccuracies have little impact on later calculations.

In OrbitalElements.jl, we have taken care to optimize the calculation time and memory allocations. While individual frequency

<footnote n="8">⁸In practice, the smoothness of the function dictates the error in the resulting calculations, rather than the integration scheme.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Example of $\Theta(w)$ (equation A6) for different $(a, e)$ in the isochrone potential, normalized to $\Theta(w = 1) = 1$, and with $b_c$ the potential scale radius. The range of $w$ runs from -1 (pericentre) to 1 (apocentre). The curves are smooth and straightforward to integrate via low-order schemes. The line for $a/b_c = 0.001$ and $e = 0.1$ is not shown because it corresponds to an orbit for which interpolation is used (see Appendix A1.1).</figcaption></figure>

<table><caption>**Table A1.** Summary of the control parameters in OrbitalElements.jl. Further description is given throughout Appendix A1.</caption><thead><tr><th>Parameter</th><th>Data type</th><th>Default</th><th>Description</th></tr></thead><tbody><tr><td>EDGE</td><td>Float64</td><td>0.01</td><td>Tolerance in <i>w</i> before switching to pericentre and apocentre expansions in equation (A6)</td></tr><tr><td>NINT</td><td>Int64</td><td>32</td><td>Number of steps for the Simpson's 1/3 integration rule in <i>w</i> to compute equations (A7)</td></tr><tr><td>TOLECC</td><td>Float64</td><td>0.01</td><td>Tolerance in eccentricity before expanding the energy, actions and frequencies from circular or radial orbits (Section A1.1)</td></tr><tr><td>TOLA</td><td>Float64</td><td>0.01</td><td>Tolerance in semi-major axis before expanding the energy, actions and frequencies from <i>r</i> = 0 (Section A1.1)</td></tr><tr><td>ITERMAX</td><td>Int64</td><td>100</td><td>Number of Newton-Raphson descent steps for backward mappings (Section A1.2)</td></tr><tr><td>inv<sub>e</sub></td><td>Float64</td><td>1e-12</td><td>Target accuracy on Newton-Raphson inversions (Section A1.2)</td></tr><tr><td>Ω<sub>0</sub></td><td>Float64</td><td>1.0</td><td>Frequency normalization scale</td></tr><tr><td>da</td><td>Float64</td><td>1.e-4</td><td>Step size for numerical derivative w.r.t. the semi-major axis (Section A1.2)</td></tr><tr><td>de</td><td>Float64</td><td>1.e-4</td><td>Step size for numerical derivative w.r.t. the eccentricity (Section A1.2)</td></tr><tr><td>rc</td><td>Float64</td><td>1.0</td><td>Characteristic scale radius</td></tr><tr><td>rmin</td><td>Float64</td><td>0.0</td><td>Inner truncation radius for the resonance variables (<i>u</i>, <i>v</i>) (Section A1.3)</td></tr><tr><td>rmax</td><td>Float64</td><td>Inf</td><td>Outer truncation radius for the resonance variables (<i>u</i>, <i>v</i>) (Section A1.3)</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure A2.** Accuracy of the numerical frequency calculation for the (analytical) isochrone model, using the distance $\Delta_{\alpha\beta}$ (equation A8). Extremely radial orbits are the most difficult ones to deal with.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure A3.** Accuracy of the numerical inversion from $(\alpha, \beta) \to (a, e)$ for the Plummer model. For each grid point in $(a, e)$, we first compute the frequencies $(\alpha, \beta)$, then 'invert' these frequencies back to calculate $(\tilde{a}, \tilde{e})$. We define a relative accuracy metric $\Delta_{ae}$ in equation (A9), which is the colourmap. The effective inversion is stable and shows a very satisfactory overall performance. The accuracy owes primarily to the recovery of eccentricity; semimajor axis is generally recovered to a higher precision.</figcaption></figure>

calculations are completed in NINT integration steps, the inversion to $(a, e)$ can take up to ITERMAX longer. As a testimony of performance, the grid of $10^6$ $(a, e)$ values and their inversions in Fig. A3 took ~140s on a single core of an M1 Macbook Air.

## A2 AstroBasis.jl: basis computation

Fundamental to the matrix method are the chosen basis functions. The AstroBasis.jl library is an implementation of several different radial

basis functions, $r \mapsto U_p^\ell(r)$, with a straightforward interface.<fnref n="9" /> At present, AstroBasis.jl supports the bases from Clutton-Brock (1972),

<footnote n="9">⁹Generically, the radial functions, $U_p^\ell(r)$, can always be taken to be real (see section 3.2 in Heyvaerts 2010). Conveniently, this makes $G(u)$ a purely real function, as in equation (A10).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure A4.** Lowest order $ℓ = 2$ radial basis function for the Clutton-Brock (1973) basis as a function of radius. Three different values of rb are shown, $r_b/b_c = [1, \sqrt{3}, 3]$. The dashed grey curve is the ROI mode for the $r_a/b_c = 0.75$ Plummer model (see Fig. 1), which peaks at the same location as the $r_b/b_c = \sqrt{3}$ basis. The dotted grey curve is the ROI mode for the $r_a/b_c = 1.0$ isochrone model, which peaks at the same location as the $r_b/b_c = 3$ basis. Using a basis tailored for the underlying mean potential, as is the case here, greatly improves the convergence of the linear response prediction.</figcaption></figure>

<table><caption>**Table A2.** Summary of the control parameters in AstroBasis.jl. Further description is given in the text.</caption><thead><tr><th>Parameter</th><th>Data type</th><th>Default</th><th>Description</th></tr></thead><tbody><tr><td>rb</td><td>Float64</td><td>1.0</td><td>Radial scale of the basis functions</td></tr><tr><td>G</td><td>Float64</td><td>1.0</td><td>Gravitational constant in Poisson equation</td></tr><tr><td>nradial</td><td>Int64</td><td>none</td><td>Maximum radial order of the expansion</td></tr></tbody></table>

Clutton-Brock (1973), Kalnajs (1976), Fridman & Poliachenko (1984)/Weinberg (1989) (Bessel), and Hernquist & Ostriker (1992). The spherical calculations in this paper use the basis from Clutton-Brock (1973), which is re-described in appendix B of Fouvry et al. (2021). The associated lowest order function matches the Plummer potential exactly. The disc calculations in this paper use the basis from Clutton-Brock (1972).

The user must choose parameters to define the basis, including the gravitational constant $G$ ($G$), the scaling radius for the basis $r_b$ (rb), the maximum harmonic number $ℓ_{max}$ ($l_{max}$), and the maximum radial number $n_{max}$ (nradial). Basic parameters are listed in Table A2. Some bases may require additional parameters.

When defining a basis, the prefactors are computed in advance and tabulated so that calls to evaluate the basis are rapid. The functions at a given radius are evaluated on the fly. For bases computed by recursion, all radial orders are computed in one step to optimize memory usage.

Fig. A4 shows an example of basis functions from AstroBasis.jl. In this case, we compare the shape of the first $ℓ = 2$ element of Clutton-Brock (1973)'s basis to the shape of the ROI modes from the isochrone cluster (Appendix C1) and the Plummer cluster (Section 3.1). To give a sense of the similarity of the basis function to the predicted ROI mode shape, we chose the rb value to match the peak of the predicted shape for the mode. In general, for linear response calculations, one decides the maximum order for the expansion to resolve the expected structure of the mode. In many cases, the user will not have an a priori guess for the shape of the mode. This entails experimentation with rb and nradial.

In addition to analytic bases for a select few models, recent basis function implementations have created an opportunity for tailored

basis functions (Petersen et al. 2022; Lilley & van de Ven 2023). The applicability of these tailored basis functions to matrix method calculations should be explored further in future work.

## A3 FiniteHilbertTransform.jl: specialized Legendre integration

FP22 describes a method to adapt Landau's integral prescription to the case of self-gravitating systems in order to approximate the integral

$$M(\omega) = \int_{-1}^{1} du \frac{G(u)}{u - \omega}, \qquad (A10)$$

at all points in the complex plane, $\omega \in \mathbb{C}$, using Legendre polynomials, $P_k(u)$. We note that we must compute integrals of this form for each resonance n and basis element pairs (p, q), as highlighted in equation (3).

To compute equation (A10), we first approximate its numerator via

$$G(u) \approx \sum_{k=0}^{K_u-1} a_k P_k(u), \qquad (A11)$$

where, importantly, $K_u$ (defined as $K_u$ in FiniteHilbertTransform.jl) controls the quality of the approximation. In practice, the coefficients $a_k$ entering equation (A11) are estimated through a Gauss-Legendre quadrature with $K_u$ nodes (see equation D3 in FP22).<fnref n="10" />

Then, the integral from equation (A10) is computed through

$$M(\omega) = \sum_{k=0}^{K_u-1} a_k D_k(\omega), \qquad (A12)$$

where

$$D_k(\omega) = \int_{-1}^{1} du \frac{P_k(u)}{u - \omega}. \qquad (A13)$$

The function $D_k(\omega)$ explicitly implements Landau's prescription and is defined in appendix D of FP22. In brief, $D_k(\omega)$ is proportional to Legendre functions of the second kind when $Im[\omega] > 0$, with an added contribution when $Im[\omega] \le 0$ (cf. section 5.2.4 of Binney & Tremaine 2008). The library FiniteHilbertTransform.jl implements the computation of (i) the $a_k$ coefficients for a given function $G$, and (ii) the $\omega \to D_k(\omega)$ functions from equation (A12) at all points in the complex plane.

The present method has mainly been introduced to probe linear response for damped frequencies, i.e. in the lower half of the complex plane. In this regime, Landau's prescription requires to give a meaning to $G(\omega)$ with $\omega \in \mathbb{C}$. Phrased differently, one has to perform an analytical continuation of these $G$ functions, which involve intricate, non-analytically known functions in the present self-gravitating case. Given that analytical continuation is intrinsically a (severely) ill-conditioned numerical problem (Trefethen 2020), for damped frequencies, i.e. negative $Im[\omega]$, the effective numerical precision plays an important role in setting the floor for accuracy.

Indeed, the complex Legendre polynomials, $P_k(\omega)$, diverge in the lower-half of the complex frequency plane. In practice, this divergence only gets approximately cancelled out by the decaying

<footnote n="10">¹⁰For smooth enough $G(u)$, this is a robust approximation scheme (Trefethen 2019).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure A5.** Convergence of finite Hilbert Transform integral calculation for the $ℓ = 2$ radially-anisotropic Plummer cluster from Fig. 1, fixing $r_a/b_c = 1.0$, $n = (-1, 2)$, $p = q = 1$, and using the Fiducial configuration settings. Here, the response matrix is evaluated for $\omega = 0 + i\gamma$, with $\gamma < 0$ (i.e. damped frequencies; solid curves) and zero oscillation frequency. Dashed curves show the corresponding calculation for $-\gamma$, i.e. unstable frequencies. Panel (a) shows that the coefficients $a_k$. They decay for small $k$, but they eventually saturate, as a result of finite numerical accuracy. Panel (b) shows that for all $k > 20$, the quantity $D_k$ increases exponentially as a function of $k$. Panel (c) shows the individual components of $a_k D_k$ as a function of $k$. For $k$ large enough, this product does not decay anymore. Finally, panel (d) shows the cumulative sum in equation (A12) as a function of $k$. For $\gamma$ sufficiently close to the real frequency line, the sum is convergent. However, as Im[$\omega$] becomes more negative, the sum begins to diverge for smaller and smaller $k$. In contrast to the solid lines which show the damped frequencies (i.e. $\gamma < 0$), the dashed lines—which show the unstable frequencies (i.e. $\gamma > 0$)—are always convergent.</figcaption></figure>

coefficients $k \mapsto |a_k|$. Due to the limitations of finite numerical precision, this results in ringing and numerical saturation in the lower half of the complex frequency plane (see e.g. Fig. 5). Fig. A5 demonstrates why this numerical saturation takes place. For values of Im[$\omega$] becoming increasingly more negative, the sum in equation (A12) diverges more and more (panel d), owing to the combination of ultimately non-decreasing $|a_k|$ (panel a), and exponentially increasing $|D_k(\omega)|$ (panel b). We emphasize that numerical saturation is eventually unavoidable given that the $G$ functions are not analytically known, whatever the chosen method.

Implementation of strategies to mitigate this numerical divergence, i.e. to push the sum from equation (A12) further down in the complex frequency plane, are currently underway. To this extent, delaying the saturation of the decomposition coefficients, $a_k$, would prove useful. It requires an enhanced accuracy on the $u \mapsto G(u)$ functions, i.e. on (i)

the OrbitalElements.jl mappings and associated gradients, and (ii) the Fourier transform of basis elements. This would expectantly increase precision in both the upper-half and lower-half planes. Though, these series of coefficients will inevitably saturate at machine precision (at best). Therefore, appropriate, a posteriori, regularization procedures should also be designed to push significantly further down in the complex plane, e.g., by reducing the effective range of summation over $k$ in equation (A12).

## A4 LinearResponse.jl: linear response computation

The computation of the response matrix and associated by-products is performed by LinearResponse.jl. It mainly requires the user to provide (i) the considered gravitational potential (and two derivatives), (ii) the DFs (through its directional derivatives $n \cdot \partial F / \partial \mathbf{J}$), and (iii) a bi-orthogonal basis. Some of these are available via OrbitalElements.jl and AstroBasis.jl, but the user can easily supply their own.

For a given harmonic $\ell$ and for each resonance $n$ and each matrix element $(p, q)$, the calculations proceed in three phases. The first two aim at computing the $u \mapsto G(u)$ functions involved in equation (3), namely (i) by computing the Fourier transform of basis elements, and (ii) by performing an integral along the resonance line, i.e. over the resonance variable $v$. The third phase is to decompose these functions over Legendre polynomials, through the computation of the $a_k$ coefficients from equation (A11) using FiniteHilbertTransforrm.jl. Once these coefficients are known, the response matrix can be efficiently computed for any given complex frequency $\omega$.

### A4.1 Computing the $G$ functions

equation (1) reduces to equation (3) when taking

$$G(u) = \int_{0}^{1} dv' \frac{2}{\Omega_{0}(\omega_{n}^{\text{max}} - \omega_{n}^{\text{min}})} \left| \frac{\partial \mathbf{J}}{\partial(u, v')} \right| G(\mathbf{J}[u, v]), \quad (\text{A14})$$

where $\omega_n^{\text{min}}$ and $\omega_n^{\text{max}}$ are introduced in Appendix A1.3 and the $\mathbf{J} \mapsto G(\mathbf{J})$ functions can be found in Appendix B. In both cases, they involve the Fourier transform of basis elements (Tremaine & Weinberg 1984). They read

$$W_{p}^{\ell n}(\mathbf{J}) = \int_{-1}^{1} dw \frac{d\theta_{1}}{dw} U_{p}^{\ell}(r) \cos(n_{1}\theta_{1} + n_{2}[\theta_{2} - \phi]), \quad (\text{A15})$$

where the radius $r$ and the angles $\theta_1$ and $\theta_2 - \phi$ are implicit functions of the orbit, $\mathbf{J}$, and the Hénon anomaly, $w$, introduced in equation (A5). Following appendix B of Rozier et al. (2019), we perform the nested integration over $w$ in equation (A15) simultaneously, pushing the angles with

$$\frac{d}{dw}(\theta_1, \theta_2 - \phi) = (\Omega_1, \Omega_2 - L/r^2) \Theta(w), \quad (\text{A16})$$

where $\Theta$ follows from equation (A6). In practice, we use the RK4 scheme (see e.g. Press et al. 2007) with $Kw$ steps. Importantly, we perform the integration backward from apocentre $(w, \theta_1, \theta_2 - \phi) = (1, \pi, 0)$ to pericentre to mitigate errors.

For a fixed accuracy goal, more eccentric orbits ($e \to 0$) typically require more integration steps. We therefore add the ability to choose the number of steps adaptively, following

$$Kw \leftarrow \lceil Kw / (0.1 + (1 - e)) \rceil. \quad (\text{A17})$$

This option is set by the boolean parameter ADAPTIVEKW.

Computing the Fourier transforms at each point $(u, v')$ and for each resonance $n$ is the most intensive part of the linear response
<table><caption>**Table A3.** Summary of the control parameters in *LinearResponse.jl*. Further description is given throughout Appendix A4.</caption><thead><tr><th>Parameter</th><th>Data type</th><th>Default</th><th>Description</th></tr></thead><tbody><tr><td>Ku</td><td>Int64</td><td>200</td><td>Number of points to sample per resonance in <i>u</i> coordinate</td></tr><tr><td>Kv</td><td>Int64</td><td>200</td><td>Number of points to sample per resonance in <i>v</i> coordinate</td></tr><tr><td>Kw</td><td>Int64</td><td>200</td><td>Number of points to sample per orbit to compute Fourier transform of basis elements (equation A15)</td></tr><tr><td>ADAPTIVEKW</td><td>Bool</td><td>false</td><td>If true, automatically scale number of Kw points to optimize calculation (equation A17)</td></tr><tr><td>VMAPN</td><td>Int64</td><td>1</td><td>Exponent in the <i>v</i>→<i>v</i>' mapping from equation (A18). Larger values will sample more finely near <i>v</i><sub>n</sub><sup>-</sup></td></tr><tr><td>lharmonic</td><td>Int64</td><td>2</td><td>Harmonic <i>l</i> to be considered</td></tr><tr><td>nlmax</td><td>Int64</td><td>10</td><td>Maximum radial resonance number <i>n</i><sub>1</sub> to consider in equation (3) (azimuthal number <i>n</i><sub>2</sub> constrained by lharmonic)</td></tr></tbody></table>

computations. For the Fiducial configuration settings, each resonance takes approximately 16 seconds on an M1 Macbook Air. The calculation of the resonances can be parallelized. Helpfully, these quantities do not depend on the DFs but only on the potential. They are therefore computed once and stored in HDF5 files for future use.

The integral along the resonance line in equation (A14) is performed on a rescaled variable $v'$, compared to the choice $v = \alpha$ (for $n_2 \neq 0$) used in FP22. Indeed, in practice, we use

$$v' = \left[ \frac{v - v_{\text{n}}^{-}}{v_{\text{n}}^{+} - v_{\text{n}}^{-}} \right]^{\text{VMAPN}}, \quad (A18)$$

with $v \in [v_{\text{n}}^{-}, v_{\text{n}}^{+}]$ and VMAPN some integer parameter which allows for the variable to spread more evenly along the resonance line. This proves particularly useful for the ILR. Ultimately, the integral over $v'$ is performed using a simple mid-point rule with Kv points.

In practice, each function $u \mapsto G_{pq}^{\ell n}(u)$ is evaluated in Ku points, $\{u_k\}_k$, chosen according to the Gauss-Legendre quadrature. For the Fourier transform of basis elements, for each resonance **n**, we end up with $Ku \times Kv$ vectors $\{W_p^{\ell n}(u, v)\}_p$ of length nradial. Computing these vectors is the most memory- and time-consuming operation. We therefore cache files for each considered resonance **n**, in arrays of size $Ku \times Kv \times nradial$. For the $G$ functions, for each resonance **n**, we end up with $nradial \times nradial$ vectors $\{G_{pq}^{\ell n}(u_k)\}_k$ of length Ku. These are also cached.

## A4.2 Computing the response matrix

Once the $G$ functions tabulated, their Legendre decomposition coefficients from equation (A11) are computed using FiniteHilbertTransform.jl. For each resonance **n**, this results in Ku matrices of size $nradial \times nradial$.<fnref n="11" />

Then, for a given complex frequency $\omega$, the response matrix from equation (3) is computed using equation (A12). We write

$$M^{\ell}(\omega) = \sum_{\mathbf{n}} \sum_{k} D_k(\omega_{\mathbf{n}}) A_{k}^{\ell n}, \quad (A19)$$

where $\omega \mapsto \omega_{\mathbf{n}}(\omega)$ is given by equation (11) in FP22, and the $D_k$ functions are computed by via FiniteHilbertTransform.jl.

On top of this computation of the response matrix, *LinearResponse.jl* also supports a range of analysis steps, namely:

(i) Mapping the complex plane by computing the determinant of the (gravitational) dielectric matrix from equation (2) at a list of $\omega$ values, as illustrated in Fig. 5. When this determinant vanishes, the system supports a mode.

(ii) Identifying the frequency of a mode, $\omega_M$, using a gradient descent algorithm.¹²This is the method used in the main text to report on the modes' frequencies.

(iii) Identifying the shape of a mode, $\psi_M$, given its frequency. We write

$$\psi_M(\mathbf{r}) = \mathrm{Re} \left[ \sum_p X_p \psi_p(\mathbf{r}) \right], \quad (A20)$$

where $\{X_p\}_p$ is the kernel eigenvector of the dielectric matrix, and $\psi_p(\mathbf{r}) \propto U_p^{\ell}(r)$ (with an angular dependence set by the considered geometry) runs over the potential basis functions used in the response matrix computation.

## A4.3 Linear response parameters

As for *OrbitalElements.jl*, control parameters for *LinearResponse.jl* are set in a JULIA structure taken as an optional argument of all exposed functions. In addition, all control parameters from *OrbitalElements.jl* are passed to *LinearResponse.jl* as a structure, owing to the computations offloaded to *OrbitalElements.jl*. In general we find that the default settings result in high enough accuracy for our problems of interest. The relevant control parameters are summarized in Table A3.

The accuracy of *LinearResponse.jl* is not straightforward to predict. Let us nonetheless report on some heuristics. First, increasing the number nradial of basis elements leads to more oscillating radial functions, hence more difficult reconstructions. This requires therefore a more precise quadrature. The same with nlmax, the maximum radial resonance number. As it increases, the resonance structures become increasingly complicated and require additional care when estimating integrals.

To conclude this section, let us comment on performance. The generic computation of the self-gravitating linear response involves numerous nested integrals, from frequency computations to finite Hilbert transforms. For these intense computations, JULIA provides an ideal language to optimize the running time while keeping a high readability. As an example, for a given $(\ell, \mathbf{n})$, a typical computation of all the needed values of $\{W_p^{\ell n}(u, v)\}_p$, with 'Fiducial' parameters and nradial = 100, takes ~40 s on a single core of an M1 Macbook Air. Fortunately, parallelizing over resonances is straightforward in JULIA. Several resonances can be then computed at once, depending on the thread count of the hardware.

<footnote n="12">¹²This involves two new control parameters, POLETOL and NRSTEP, and requires the user to define a starting point. In practice, the pole finding only weakly depends on these control parameters, and reasonable defaults have been set.</footnote>

<footnote n="11">¹¹We use the symmetry $A_{pq} = A_{qp}$ to reduce the computational cost.</footnote>
# APPENDIX B: LINEAR RESPONSE INTEGRANDS

## B1 For spheres

As in FP22, in three-dimensional (3D) the full expression of $G^{\ell n}(\mathbf{J})$ in equation (1) is given by

$$G_{pq}^{\ell n}(\mathbf{J}) = -\frac{2(2\pi)^3}{2\ell + 1} |y_{\ell}^{n_2}|^2 L \left(\mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}}\right) W_p^{\ell n}(\mathbf{J}) W_q^{\ell n}(\mathbf{J}), \quad (B1)$$

where we introduced $y_{\ell}^{n} = Y_{\ell}^{n}(\frac{\pi}{2}, 0)$, with $Y_{\ell}^{m}(\vartheta, \phi)$ the spherical harmonics normalized so that $\int d\vartheta d\phi \sin(\vartheta)|Y_{\ell}^{m}(\vartheta, \phi)|^2 = 1$. Further details may be found in FP22, including the Jacobian for the conversion to $G^{\ell n}(u, v)$, needed in practice by LinearResponse.jl. In equation (B1), the Fourier-transformed basis elements follow from equation (A15).

## B2 For discs

It is straightforward to extend the linear response computations of FP22 to the case of razor-thin discs with axial symmetry. In this case, the definition of the response matrix (see section 5.3 in Binney & Tremaine 2008) immediately applies, without any integral over a third dummy action as in the 3D spherical case (see section 4 in Hamilton et al. 2018). Then, simply using

$$G_{pq}^{\ell n}(\mathbf{J}) = -(2\pi)^2 \delta_{\ell}^{n_2} \left(\mathbf{n} \cdot \frac{\partial F}{\partial \mathbf{J}}\right) W_p^{\ell n}(\mathbf{J}) W_q^{\ell n}(\mathbf{J}), \quad (B2)$$

instead of equation (B1), one can compute the linear response of razor-thin discs using all the numerical tools from LinearResponse.jl. Equation (B2) imposes the constraint $n_2 = \ell$ and the resonance $\mathbf{n} = (0, 0)$ does not contribute (to the $\ell = 0$ response matrix) as in the spherical case. It also involves the Fourier transform of 2D bi-orthogonal basis elements, following equation (A15), which are described in Appendix F1.

# APPENDIX C: VALIDATION WITH THE ISOCHRONE SPHERE

In this section, we revisit the calculations from FP22 for the spherical isochrone potential (Hénon 1959). Because the isochrone has analytic expressions for the frequencies, one can readily evaluate integrals to high numerical precision, as used in FP22. However, here we exercise the empirical frequency computation from OrbitalElements.jl (Appendix A1) for our tests of the isochrone potential. In what follows in this section, we use our numerical calculations of the frequencies to reproduce earlier results, testing the ROI and $\ell = 1$ damped mode calculations. The equations for the isochrone model potential and DFs can be found in FP22.

## C1 $\ell = 2$ modes – ROI

The ROI in the isochrone model with $r_a/b_c = 1.0$, with $b_c$ the isochrone scale radius, has become a regular test for different aspects of the matrix method implementation (Hamilton et al. 2018; Fouvry et al. 2021; Rozier et al. 2022; Lilley & van de Ven 2023). Using the radially-anisotropic DFs given in equation (G12) of Fouvry et al. (2021), we compute the growth rate for the fiducial $r_a/b_c = 1.0$ case. We find a fastest growing mode with $\gamma_M/\Omega_0[r_a/b_c = 1.0] = 0.0236 \pm 0.0015$. This matches the result of FP22 at the 2 per cent level, and Saha (1991) at the 4 per cent level. Additionally, the shape

of the recovered isochrone ROI mode matches that of FP22, which itself matches the mode shape of Saha (1991).

Given the efficiency of LinearResponse.jl, we are able to test a range of configuration parameters to determine how results may vary with control parameters. This allows us to establish some measure of uncertainty, as quoted in the above paragraph. In Table C1, we report the growth rate of the fastest growing mode with various control parameters.

Previous studies have investigated the dependence of the ROI w.r.t. the normalized anisotropy radius, $r_a/b_c$. A naive application of our libraries and the DFs given in Fouvry et al. (2021) suggests that the isochrone stability boundary is $r_a/b_c > 2$. However, in examining the results of our methodology, we discovered a few numerical pitfalls that could alter results. First, the integration domain should be appropriately limited: as $r_a/b_c$ decreases, regions of the $(J_r, L)$ plane become unpopulated (see e.g. fig. 4 in Hamilton et al. 2018). Second, the gradient of the DFs becomes extremely large at the $Q = 0$ boundary, such that numerical estimates readily break down. To circumvent this, Hamilton et al. (2018) applied a smoothing to the radially-anisotropic form of the DFs hence considering the smoothed DFs $\hat{f} = \exp(-\zeta/Q)f$. We find here that different values of $\zeta$, the smoothing length, result in significantly different estimates for the growth rate, at fixed $r_a/b_c$. We do not address either of these shortcomings in the present work. As a result, obtaining a clear and reliable stability transition radius in $r_a/b_c$ remains out of reach. We note that even though the growth frequency may be mis-estimated, the shape of the predicted modes is converged.

For completeness, we finally report literature values of the stability boundary in $r_a/b_c$. May & Binney (1986) estimated an instability threshold at the half-mass radius of $r_a/b_c = 1.67$. Merritt (1988) then identified an instability threshold of $r_a/b_c = 1.9$, from simulations using a harmonic code (Merritt 1987). Saha (1991) performed new simulations using the same harmonic code from Merritt (1987), producing qualitative, but not quantitative, agreement in the trend of growth rates. The linear response calculations in Saha (1991), using a tailored basis, suggested that unstable modes may be found for $r_a/b_c \le 4$. Future work may return to the question of stability with improved numerical treatment.

## C2 $\ell = 1$ modes

While searches for instabilities have been fairly commonplace in dynamical modelling, searches for damped modes in stellar clusters have not equally flourished. Most of this owes to the challenge of calculating the linear response of a system in the lower-half plane. Weinberg (1994) identified $\ell = 1$ damped modes in King (1966) models using linear response theory. Weinberg (1994) also compared this prediction with perturbed N-body simulations, finding qualitative agreement. Heggie et al. (2020) measured the same $\ell = 1$ mode in N-body simulations of a King model.

Extending the tools developed for the isochrone cluster, FP22 measured a damped mode for the isotropic isochrone model at $\omega_M/\Omega_0 = 0.0143 - 0.00142i$. As a first test, we reproduced this prediction using all the machinery from LinearResponse.jl. Table C2 lists results for a set of validation runs used to search for the $\ell = 1$ damped mode. There, we denote the set of parameters from FP22, for which we recover $\omega_M/\Omega_0 = 0.0143 - 0.00143i$, in tight agreement with FP22.

However, we noted that varying the control parameters of LinearResponse.jl results in a drift of the mode's frequency. The middle column of Table C2 lists the results for different runs. As the number of resonances considered increases (controlled by $n1max$),
<table><caption>**Table C1.** Predictions for the ROI $\ell = 2$ growth rates in the isochrone ($r_a/b_c = 1.0$) and Plummer ($r_a/b_c = 0.75$) clusters. The configurations with increasing $n1max$ are the most computationally expensive: $n1max \times 2$ requires 122 resonances, while $n1max \times 4$ requires 242 resonances.</caption><thead><tr><th>Run</th><th>Isochrone γ<sub>M</sub>/Ω<sub>0</sub></th><th>Plummer γ<sub>M</sub>/Ω<sub>0</sub></th></tr></thead><tbody><tr><td>Fiducial</td><td>0.0236</td><td>0.0477</td></tr><tr><td>Ku×2</td><td>0.0247</td><td>0.0477</td></tr><tr><td>Kv×2</td><td>0.0243</td><td>0.0477</td></tr><tr><td>Kw×2</td><td>0.0238</td><td>0.0477</td></tr><tr><td>nradial/2</td><td>0.0238</td><td>0.0470</td></tr><tr><td>n1max=1</td><td>0.0238</td><td>0.0470</td></tr><tr><td>n1max×2</td><td>0.0237</td><td>0.0478</td></tr><tr><td>n1max×4</td><td>0.0238</td><td>0.0479</td></tr><tr><td>rb/5</td><td>0.0236</td><td>0.0476</td></tr><tr><td>rb×4 (FP22)</td><td>0.0238</td><td>0.0477</td></tr><tr><td>VMAPN = 2</td><td>0.0251</td><td>0.0477</td></tr></tbody></table>

<table><caption>**Table C2.** Predictions for $\ell = 1$ mode locations in the complex plane. The search space is Re($\omega$)/$\Omega_0 \in [0.0, 0.1]$ and Im[$\omega$]/$\Omega_0 \in [-0.01, 0.05]$. In general, the more resonances are considered, the more the complex frequency of the recovered (damped) $\ell = 1$ mode converges to the frequency of the (neutral) translation mode $\omega_M = 0 + 0i$.</caption><thead><tr><th>Run</th><th>Isochrone ω<sub>M</sub>/Ω<sub>0</sub></th><th>Plummer ω<sub>M</sub>/Ω<sub>0</sub></th></tr></thead><tbody><tr><td>Fiducial</td><td>0.0143 – 0.0014i</td><td>0.0247 – 0.0004i</td></tr><tr><td>Ku×2</td><td>0.0150 – 0.0022i</td><td>0.0247 – 0.0004i</td></tr><tr><td>Kv×2</td><td>0.0143 – 0.0014i</td><td>0.0247 – 0.0004i</td></tr><tr><td>Kw×2</td><td>0.0143 – 0.0014i</td><td>0.0247 – 0.0004i</td></tr><tr><td>nradial/2</td><td>0.0143 – 0.0014i</td><td>0.0246 – 0.0004i</td></tr><tr><td>n1max×2</td><td>0.0084 – 0.0005i</td><td>0.0124 – 0.0001i</td></tr><tr><td>n1max×4</td><td>0.0048 – 0.0002i</td><td>0.0060 – 0.0001i</td></tr><tr><td>rb/5</td><td>0.0143 – 0.0014i</td><td>0.0247 – 0.0004i</td></tr><tr><td>rb×4 (FP22)</td><td>0.0143 – 0.0014i</td><td>0.0247 – 0.0004i</td></tr><tr><td>VMAPN=2</td><td>0.0143 – 0.0014i</td><td>0.0246 – 0.0004i</td></tr></tbody></table>

the predicted complex frequency of the $\ell = 1$ mode shifts towards the origin of the complex plane. On the bright side though, when holding other parameters fixed and varying only the scale length of the basis (rb), we find that the result is essentially converged: the dilatation of the basis can be absorbed in this problem.

For the calculations considered here, the cluster's pure neutral translation mode (the $\ell = 1$ mode at $\omega_M = 0 + 0i$, see Weinberg 1989; Murali 1999) is not converged, and appears as a 'false' damped mode. This seems to be the mode we are recovering here. We do not find evidence for any other $\ell = 1$ modes in the isotropic isochrone model.

False modes are more likely to manifest themselves in the lower half-plane due to the Landau prescription employed in computing the integral presented in equation (A13). Indeed, the divergence of $P_k(\omega)$ (equation (A13)) only contributes in the lower-half of the complex frequency plane.

Regrettably, due to numerical noise in the calculation, we cannot explore sufficiently low regions in the complex plane to identify authentic damped modes that might exist below our imposed limits. Hopefully, upcoming numerical enhancements should enable us to delve deeper into the complex frequency plane. Finally, let us stress that the presence of very damped modes are unlikely to be relevant in Nature. Indeed, in a $N = 10^5$ globular cluster, the evolution of fluctuations can be taken to be collisionless only for large enough amplitude, hence making globular clusters unable to 'resolve' very damped linear modes.

## APPENDIX D: THE PLUMMER SPHERE

Unlike the isochrone sphere, the Plummer sphere does not have analytical expressions for its orbital frequencies. We therefore use the tools in OrbitalElements.jl to compute the frequencies from the potential (and its two derivatives), as well as the resonantly-aligned $(u, v)$ coordinates. Here, the versatility of LinearResponse.jl shines: the user only needs to define the DFs.

The Plummer (1911) potential is given by

$$ \psi(r) = -GM/\sqrt{b_c^2 + r^2}, \qquad (D1) $$

with $b_c$ the scale radius. The frequency scale for the Plummer sphere, $\Omega_0$, is given by

$$ \Omega_0 = \Omega_1(r \to 0) = 2\sqrt{GM/b_c}. \qquad (D2) $$

The isotropic DFs is

$$ F(E) = \frac{24\sqrt{2}}{7\pi^3} \left( \frac{E}{E_0} \right)^{7/2}, \qquad (D3) $$

with the energy scale $E_0 = -GM/b_c$.

To introduce radial anisotropy, we use the Osipkov (1979); Merritt (1985) method to find a self-consistent DFs. The dimensionless quantity

$$ Q = \frac{1}{E_0} \left( E + \frac{L^2}{2r_a^2} \right) \qquad (D4) $$

is the single parameter in the DFs, which reads<fnref n="13" />

$$ F(Q) = \frac{24\sqrt{2}}{7\pi^3} Q^{7/2} \left[ 1 - \frac{1}{r_a^2} + \frac{7}{16r_a^2 Q^2} \right]. \qquad (D5) $$

## APPENDIX E: N-BODY TECHNIQUE

To compare our linear response theory predictions with N-body simulations, we generate realizations of the Plummer cluster and evolve them for several hundred dynamical times. The results are discussed in Section 3.3. This appendix gives a brief description of the N-body method, and the parameters of the N-body runs (see Table E1).<fnref n="14" />

The initial conditions are generated using PlummerPlus (Breen et al. 2017). We use a lightly modified version of PLUMMERPLUS.PY, specifically to produce outputs that interface directly with EXP (Petersen et al. 2022), a basis function expansion N-body code.<fnref n="15" /> The initial conditions are centred by the average positions and velocities in each dimension after sampling, resulting in a zero-momentum initial setup.

Simulations are run using EXP. For the basis, we use the Sturm-Liouville solution from EXP. In the case of the Plummer model, the Sturm-Liouville solver in EXP results in the same basis as the analytical Clutton-Brock (1973) basis (cf. Appendix A2).

Conveniently, EXP uses virial units, making direct comparison between the simulation and linear response calculations straightforward. We run simulations for 400 time units, with a maximum possible timestep of 0.2. In practice, EXP employs a block-step technique for subdividing timesteps in powers of 2. We allow a

<footnote n="13">¹³Equation (5) in Breen et al. (2017) does not include a square on the middle term in brackets. This appears to be a typesetting error, as the PLUMMERPLUS.PY software (see Appendix E) includes this square.</footnote>

<footnote n="14">¹⁴Any unspecified parameters are taken as software defaults.</footnote>

<footnote n="15">¹⁵The code we use is available at https://github.com/michael-petersen/PlummerPlus, forked from the original.</footnote>
<table><caption>**Table E1.** Summary of the control parameters for then EXP N-body runs.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>lmax</td><td>6</td><td>Number of azimuthal harmonics</td></tr><tr><td>nmax</td><td>24</td><td>Number of radial basis functions</td></tr><tr><td>N</td><td>262 144</td><td>Number of equal-mass particles</td></tr><tr><td>multistep</td><td>7</td><td>Maximum subdivision of the time-step</td></tr><tr><td>dt</td><td>0.2 [virial]</td><td>With multistep, minimum time-step is 0.2/2<sup>7</sup></td></tr><tr><td>nsteps</td><td>2 000</td><td>With dt, total run time is 400 virial units</td></tr></tbody></table>

maximum of 7 divisions, corresponding to a minimum timestep of ~1.5 × 10⁻³. The median timestep for particles is 0.2/2³ = 0.025.

As is the default behaviour in EXP, the centre of the expansion is pinned at the inertial centre, leaving to the basis to resolve any excursions. As the maximum excursion of the peak density is well within a scale radius of the simulated clusters, we do not expect any bias (cf. Petersen et al. 2022).

For computational expediency, we draw a relatively small number of particles per cluster,<fnref n="16" /> but realise six clusters per DFs to estimate the average response. For the ROI measurement in Fig. 3, we shift time such that the peak potential fluctuations line up.

To efficiently estimate the magnitude of the potential fluctuations, Δψ, we first compute the self-gravity in the non-axisymmetric harmonics. Given that the square of the EXP coefficients is the self-gravity (cf. section 1 of Petersen et al. 2022), and the magnitude of the amplitudes encodes Δψ, we may efficiently measure the potential fluctuations in either self-gravity or Δψ, simply using the output products of the simulation

$$ \Delta\psi^\ell = \sum_p a_p^\ell \psi_p^\ell, \qquad (E1) $$

where $a_p^\ell$ is the coefficient of the given basis function $\psi_p^\ell$. This is used in Figs 3 and 4.

To compute the energy in a given harmonic, we compute the sum square of the amplitudes of all radial orders in a given harmonic,

$$ E_{\ell} = \sum_{p} (a_{p}^{\ell})^{2}. \qquad (E2) $$

To provide a natural scaling, harmonics with $\ell > 0$ are normalized by the monopole to define

$$ \tilde{E}_{\ell} = E_{\ell} / E_{\ell=0}. \qquad (E3) $$

In Fig. 3, we measure the slope of ln($\tilde{E}_2$) to obtain the growth rate for the ROI simulations. The oscillation frequency of a mode, determined from a given coefficient for a harmonic $\ell$, is computed as

$$ \Omega_M = \arctan \left( \frac{\mathrm{Im}[a_p^\ell]}{\mathrm{Re}[a_p^\ell]} \right). \qquad (E4) $$

Measurements based on the coefficients are not perfect, in that both truncation noise and particle noise may bias the coefficients. Fortunately, these errors are generally small. In the context of the present simulations, reasonable estimates for the error on Δψ^ℓ/ψ^0 (where ψ^0 is the equilibrium model potential) and $\tilde{E}_\ell$ are of order 1 per cent.

# APPENDIX F: RAZOR-THIN DISCS

## F1 2D bi-orthogonal bases

Due to axial symmetry, it is natural to contemplate basis elements that are harmonically decoupled. Potential-density pairs for a disc then take the form $[U_p^\ell(r)e^{i\ell\phi}, D_p^\ell(r)e^{i\ell\phi}]$ where $(r, \phi)$ are the usual polar coordinates and $p$ (resp. $\ell$) stands for the radial (resp. azimuthal) number. As in the spherical case, a basis is then fully prescribed by its radial potential functions, $U_p^\ell$. In AstroBasis.jl, we implemented the 2D basis from Clutton-Brock (1972) (infinite radial support) and Kalnajs (1976) (finite radial support), with a careful treatment of the normalization. Naturally, it would be of interest to implement tailored basis functions (Petersen et al. 2022; Lilley & van de Ven 2023). This is left for future work.

## F2 Zang discs

The stability analysis of tapered Mestel discs was pioneered by Zang (1976) and subsequently generalized by Evans & Read (1998b). We refer to section 4.1 of Fouvry et al. (2015) for a concise presentation of these discs, their DFs and potentials.

We place ourselves in the same unit system as in Fouvry et al. (2015), hence setting $V_0 = G = R_i = 1$. Within this unit system, a natural scale frequency is $\Omega_0 = V_0/R_i = 1$. Here, $\Omega_0$ does not stand for the (infinite) central radial frequency. As in Zang (1976), we consider ‘fully’ self-gravitating discs (active fraction $\xi = 1$) with

<table><caption>**Table F1.** Summary of the control parameters for the disc’s linear response computations. Unspecified parameters are taken as defaults (Tables A1–A3).</caption><thead><tr><th>Library</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>OrbitalElements.jl</td><td>rmin</td><td>0.2</td></tr><tr><td rowspan="3">AstroBasis.jl</td><td>Basis</td><td>Clutton-Brock (1972)</td></tr><tr><td>rb</td><td>5.0</td></tr><tr><td>nradial</td><td>100</td></tr><tr><td>LinearResponse.jl</td><td>VMAPN</td><td>2</td></tr></tbody></table>

<table><caption>**Table F2.** Predictions for the N = 4 Zang disc most unstable mode: oscillation frequencies and growth rates when varying the control parameters. When considering enough resonances (nlmax sufficiently large), the predictions are strongly consistent. Spreading the point more evenly along the resonance line [via setting VMAPN = 2 in equation (A18)] prove particularly useful in reducing the number of integration points Kv needed to reach convergence, hence saving computational time and reducing memory allocation.</caption><thead><tr><th>Run</th><th>Ω<sub>M</sub>/Ω<sub>0</sub></th><th>γ<sub>M</sub>/Ω<sub>0</sub></th></tr></thead><tbody><tr><td>Fiducial</td><td>0.878</td><td>0.126</td></tr><tr><td>Ku×2</td><td>0.878</td><td>0.126</td></tr><tr><td>Kv×2</td><td>0.878</td><td>0.126</td></tr><tr><td>Kw×2</td><td>0.878</td><td>0.126</td></tr><tr><td>nradial/2</td><td>0.878</td><td>0.126</td></tr><tr><td>nlmax=1</td><td>0.842</td><td>0.028</td></tr><tr><td>nlmax×2</td><td>0.879</td><td>0.127</td></tr><tr><td>nlmax×4</td><td>0.879</td><td>0.127</td></tr><tr><td>rb/5</td><td>0.878</td><td>0.126</td></tr><tr><td>rb×4</td><td>0.878</td><td>0.128</td></tr><tr><td>VMAPN = 1</td><td>0.885</td><td>0.135</td></tr><tr><td>VMAPN = 1, Kv×2</td><td>0.881</td><td>0.128</td></tr><tr><td>VMAPN = 1, Kv×4</td><td>0.879</td><td>0.127</td></tr></tbody></table>

<footnote n="16">¹⁶No tests w.r.t. N were performed, though we expect that the measurement uncertainty will decrease like $\sqrt{N}$.</footnote>
radial velocity dispersion such that $q = 6$, while varying the inner taper exponent $v_t(=N)$. The outer taper is taken sufficiently far ($R_0 = 11.5$) and sharp ($\mu_t = M = 5$) not to affect the populated regions, hence not affecting the disc’s response.

For numerical purpose in OrbitalElements.jl, we consider a soft-ened Mestel potential

$$ \psi_M(r) = \frac{1}{2} V_0^2 \ln[\varepsilon^2 + (r/R_0)^2], \qquad (F1) $$

to prevent singularities for exactly radial orbits. In practice, $\varepsilon$ is taken sufficiently small ($\varepsilon = 10^{-5}$; $R_0 = 1$) compared to the DFs’s

inner taper radius $R_i = 1$. We checked that this softening had no impact on the predicted linear response. Table F1 summarizes the parameters used in LinearResponse.jl for the disc ‘Fiducial’ numerical applications. In Table F2, we report the location of the fastest growing mode for the $N = 4$ Zang disc with various control parameters. The ‘Fiducial’ run is well converged.

This paper has been typeset from a LaTeX/Latex file prepared by the author.