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

We present LinearResponse.jl, an efficient, versatile public library written in julia to compute the linear response of self-gravitating (3D spherically symmetric) stellar spheres and (2D 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.


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 is generically This equation encodes the physics of (linear) stability in selfgravitating systems.Here, n is the sum over all (allowed) resonance vectors n; ∫ L J is the scan over the populated orbital (action) space (with an appropriate Landau prescription); n• (J) − is the resonant amplification at a given (complex) frequency and for orbital frequencies ; G n are functions of the system's distribution function (DF) 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 when the condition det[ ( )] = 0, with ( ) is met, where I is the identity matrix and is called the (gravitational) dielectric matrix.We write frequencies satisfying this equality as M , which we further break down into M = Ω M +i M , with the oscillation frequency, Ω M , and the growth rate, M .Modes may appear at any point in the complex plane.Modes with M > 0 are unstable, modes with M = 0 are neutral, and modes with 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 Rĳcke & Voulis 2016) and spheres (e.g., Polyachenko & Shukhman 1981;Weinberg 1989Weinberg , 1991;;Saha 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 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, including accelerating the overall relaxation (see, e.g., Hamilton et al. 2018;Heggie et al. 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 organised 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 -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.

LINEAR RESPONSE OF SPHERES AND DISCS
The linear response of spheres and discs can be split in independent harmonics ℓ (historically denoted for the discs).For these systems and within the appropriate "resonance" coordinate system ( , ) introduced in equation (10) of FP22, each element of the response matrix of harmonic ℓ reads where the functions1 ℓn ( ) depend on the Fourier transform of bi-orthogonal basis elements (indexed by and ), the resonance number, n = ( 1 , 2 ), and involve an integral over some other coordinate d , as explained in Appendix A4.1.These functions, which depend on the dimensionality of the problem, are detailed in Appendix B. In equation (3), n ( ) stands for the rescaled (complex) frequency introduced in equation ( 11) of FP22.Importantly, the imaginary parts of the rescaled frequency n and 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 • OrbitalElements.jl(Appendix A1), which, given any central potential and its two first derivatives, performs various changes of orbital coordinates to compute the actions ( r , ), orbital frequencies (Ω 1 , Ω 2 ) and resonance coordinates ( , ).2 • 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.
• LinearResponse.jl(Appendix A4) is the driver library.Relying on the three (independent) previous libraries, LinearResponse.jlimplements the computation of the response matrix M ℓ ( ).It performs the Fourier transform of the basis elements, then computes the functions ℓ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.

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 DF 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 -body code, to test our predictions (Section 3.3).We provide the relevant equations for the Plummer model in Appendix D.

ℓ = 2 modes -Radial Orbit Instability
The ROI arises in models with some degree of radial anisotropy, i.e., the DF 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, a , outside of which the fraction of radial orbits is enhanced.The mechanisms for the instability are discussed in Palmer (1994); Polyachenko & Shukhman (2015).The classic ROI results in a quadrupole (ℓ = 2) zero oscillation frequency mode (Ω M = 0) that grows exponentially in time ( 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.
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 a / c ≃ 0.9 using the criteria of Barnes et al. (1986), with 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 a / c ≃ 1.1.Relying on direct -body simulations, Breen et al. (2017) confirmed that a / c = 0.75 was unstable, but, in contrast to Dejonghe & Merritt (1988), Breen et al. (2017) found that a / c = 1.0 was stable.
We test our linear response machinery on a range of a / c values.Fig. 1 shows the results of searching for modes using a "Fiducial" set of control parameters. 3We find a threshold for instability of a / 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 percent.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 ℓ ≠ 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 DF (see Appendix D), we perform a search in the complex plane.We find zeros in the determinant of the dielectric matrix , but these appear to be spurious non-converged realisations 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 /d , 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 /d .As expected given the results of FP22, the predicted mode appears to be the neutral 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 integrals and resonant interactions is, by design, a challenging task.

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

Anisotropic simulations, Radial Orbit Instability
To attempt to recover the ROI in the Plummer clusters, we select six different values of a / c for which to run -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 a / c , we run an ensemble of six simulations of 2 18 particles each.Fig. 3 shows the results for the a / c = [0.75,0.80, 0.85, 0.9] runs using EXP (see Appendix E).In general, the runs for any value of a / c show similar behaviour, though the dispersion between realisations increases with a / c (i.e., dispersion between realisations 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 a / c models; the qualitative agreement with the predicted mode shapes and trends (shown as dashed curves) with a / c are remarkably high, though the -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 a / c = 0.95, only three of the six runs performed appear to show any instability.Thus, a / c = 0.95 appears to be a weakly unstable model.We do not measure an instability for any of the a / c = 1.0 runs, in -body simulations.
Table 1 compares the measured growth rates from thebody simulations to the predictions from LinearResponse.jl.The agreement is quantitatively good, in that all measured growth rates are within 3 of the predicted growth rate, and both approaches recover the qualitative trend of decreasing M with increasing a / 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 scalelength of b = 1, match the basis used to run the simulation parameters (Appendix E).This ensures a straightforward comparison with the -body simulations.We find that the predicted eigenvectors (equation A20) qualitatively match the basis function coefficients in the -body simulations.
In conclusion, the -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.

Isotropic simulations, ℓ = 1 measurements
By analogy with the investigation of Section 3.2, we perform an ensemble of six simulations of Plummer spheres with an isotropic DF (see Appendix E).The left panel of Fig. 4 shows the run with time  of the (rescaled) energy, ˜ 1 ∝ ℓ=1 (Eq.E3), for each simulation.Interestingly, ˜ 1 increases and slowly approaches a saturated level that differs between realisations.Upon inspection, the amplitude ˜ 1 is dominated by the lowest-order function in each realisation, 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 /d , the results from the -body simulation are generally not consistent with d /d .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 ℓ = 1 density fluctuations further resembles d /d .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.Table 1.Comparison between the measured and predicted growth rates for the radial orbit instability in the Plummer sphere for various normalised anisotropy radii, a / c .The scale frequency, Ω 0 , is defined for the Plummer cluster (equation D2).Uncertainties on the -body runs arise from variations among different realisations (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.
Interestingly, in Fig. 4, the time evolution of the total energy in the ℓ = 1 fluctuations shows a clear sign of saturation.Rogister & Oberman (1968), equation ( 21 where ˜ ℓ is the (rescaled) energy in the ℓ = 1 gravitational fluctuations and ˜ 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 M , transients only fade away after a (few) 1/ 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 ˜ ℓ=1 in Fig. 4, we can use the ensemble average ℓ = 1 energy (the thick black curve in the left panel of Fig. 4) to find best-fitting values for ˜ th , M and the initial energy fluctuations, ˜ init = ˜ ℓ ( = 0).We find ˜ init = 1.0×10 −6 , ˜ th = 8.0×10 −3 , and M /Ω 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 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 ℓ = 1 energy to reach a thermal state.Additionally, the apparent scatter of values of ˜ th for individual runs is somewhat unexpected.
We defer further investigation to future work.

MESTEL RAZOR-THIN DISCS: VALIDATION
To emphasize the generality of LinearResponse.jl,we now consider razor-thin (2 ) discs.In that case, stars still interact through the usual Newtonian interaction potential (r, r ′ ) = − /|r − r ′ |, with the gravitational constant, but are confined to a plane.Assuming an axisymmetric mean DF, the stars' mean-field orbits can be described by the two action variables ( r , ), 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, ℓn , in equation ( 3), given in Appendix B2, and (ii) implementing 2 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.Nonetheless, in what follows, we do not use these simplifications and rather use the generic scheme from LinearResponse.jl.Following Evans & Read (1998a), the DF of the stars is tapered in the central region.The presence of an inner cut-off mimicks an unre- Complex frequency, M /Ω 0 , of the most unstable ℓ = 2 mode of truncated Zang discs ( = 6), as one varies the index of the inner taper, , with Ω 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% precision.The results prove robust to variations of the control parameters (Table F2 for the = 4 disc).
sponsive central bulge, hence introducing a reflexive boundary4 : 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 and far enough (Evans & Read 1998b).We refer to Appendix F2 for details and notably the definition of the scaling frequency Ω 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 DF 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., ℓ = 2 modes, as one varies the properties of the inner taper.5As 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, ( ) from equation ( 2), through its level contours in the complex-frequency plane.This determinant vanishes at the frequency M /Ω 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.
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.

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 DF, 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 ℓ = 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 -body simulations performed using EXP.For the radial orbit instability, we find quantitative agreement between predictions and measurements in simulations.This agreement opens avenues for validating -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 DF, we find the same unconverged neutral mode as in the isochrone model.We find no evidence for other modes.Finally, the -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 razorthin 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 ℓ = 2 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) generalising 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., Rozier et al. 2022;Dootson & Magorrian 2022, and references therein).
As for long-term relaxation, LinearResponse.jlmay 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.jlcan 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 DF (see, e.g., Reddish et al. 2022).

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.

APPENDIX A: LINEAR RESPONSE SOFTWARE A1 OrbitalElements.jl: Orbit Frequency Computation
The julia library OrbitalElements.jl is optimised to calculate high-precision quantities for orbits in a static central potential, = ( ).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 • the pericentre and apocentre radii ( per , apo ); • the effective semi-major axis and eccentricity with the radial action and = 2( − ( ))− 2 / 2 the instantaneous radial velocity; • the orbital frequencies (Ω 1 , Ω 2 ) associated respectively with the radial and azimuthal oscillations; • the frequency ratios ( , ) from which the frequencies are computed (see equations 28 and 29 in Tremaine & Weinberg 1984) with Ω 0 some given frequency scale (typically the central radial frequency in cored profiles); • the resonance-specific (i.e., dependent on n) coordinates ( , ) from FP22 (appendix B therein).
In practice, OrbitalElements.jl is centred around the effective semi-major axis and eccentricity ( , ), but straightforward conversions between different orbital labels exist as simple function calls.These change 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.6

A1.1 Forward mappings
The direct calculations compute orbital elements (frequency ratios, radial action, etc.) through integrals of the form ∫ apo per d ...However, due to the diverging integrand 1/ in equations (A4), the radius 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 with ∈ [−1, 1] the "Hénon anomaly", and = −1 (resp.= 1) corresponding to the pericentre (resp.apocentre).Using this variable, the integrand is no longer divergent and equations (A4) simply read In practice, the evaluation of Θ can become numerically unstable close to the boundaries = ±1.To cure this, we perform a second-order expansion for values of 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, , is readily and easily computed by integrating equation (A3) over .In Fig. A1, we illustrate the function ↦ → Θ( ) for the isochrone potential.For most orbits, Θ( ) 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. 7or near-circular orbits, i.e., → 0, and/or very small semimajor axes, i.e., → 0, the evaluation of Θ( ) can break down.To tackle this issue, we use analytic expressions for circular ( = 0) and central ( = 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., → 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 ( , ) ↦ → ( , ).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 da and de.Finally, the Newton-Raphson algorithm is controlled with a maximal number of iterations (ITERMAX) and a accuracy goal (inv ).The effective accuracy of the inversion is mainly determined by the orbital integration (NINT) and the finite differences (da and de).

A1.3 Resonances variables
Let us now discuss the computation of resonance-specific (i.e., dependent on n) coordinates ( , ).Within these coordinates, resonances lines ( n = n• = cst.)become straight lines ( = cst.).Their construction from the frequency ratios, ( , ), 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 DF are no fully selfconsistent.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" center.
We therefore truncate the domain in ( , ) effectively explored.This does not affect the computation of the frequencies, , but only the definition of the resonance variables ( , ).This effective truncation is set by the two parameters rmin and rmax.Then, the effective domain in ( , ) is restricted to the region between min = c (rmax) and max = c (rmin) with c ( ) = ( = , = 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 ( , ) domain is focused on orbital regions effectively populated.
Fortunately, compared to FP22, adding these truncations only  , 1 2 ), ( max , c [ max ]) or along the curve ( , c [ ]) with min ≤ ≤ max .For the variable, the first two constraints of equation (B10) in FP22 become min ≤ ; ≤ 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.Nonetheless, 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 summarised 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.jlfor this accuracy benchmark.
First, in Fig. A2, we test the accuracy of and calculations for the isochrone model (for which we know the analytic values of and ).As a summary statistics of the accuracy, we define the distance between the numerically-calculated values of and (denoted ˜ , ˜ ) and the analytic values (simply denoted , ) through In practice, Δ is dominated by the accuracy contribution from ˜ .
To probe the full range of relevant semi-major axes, we sample using 1000 log-spaced points / c ∈ [10 −3 , 10 3 ], and 1000 linearspaced points ∈ [0, 1].The accuracy is better than 1% at all points sampled, and in most cases is better than 0.001% accurate.In Fig. A3, we use the same 1000×1000 grid of ( , ) to test the inversion from frequencies back to ( , ).This time we use the Plummer potential, to test both the numerical computation of ( , ) ↦ → ( ˜ , ˜ ) and the subsequent numerical inversion In practice we find that the inversions are quite accurate, with typical differences on the order of 0.001%.However, in certain regions it becomes difficult to accurately recover the input ( , ) 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 optimise the calculation time and memory allocations.While individual frequency calculations are completed in NINT integration steps, the inversion to ( , ) can take up to ITERMAX longer.As a testimony of performance, the grid of 10 6 ( , ) 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, ↦ → ℓ ( ), with a straight- We define a relative accuracy metric Δ 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; semi-major axis is generally recovered to a higher precision.(1972).
The user must choose parameters to define the basis, including the gravitational constant (G), the scaling radius for the basis b (rb), the maximum harmonic number ℓ max (lmax), and the maximum radial number 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 optimise memory usage.Figure 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.self-gravitating case.Given that analytical continuation is intrinsically a (severely) ill-conditioned numerical problem (Trefethen 2020), for damped frequencies, i.e., negative Im[ ], the effective numerical precision plays an important role in setting the floor for accuracy.
Indeed, the complex Legendre polynomials, ( ), diverge in the lower-half of the complex frequency plane.In practice, this divergence only gets approximately cancelled out by the decaying coefficients ↦ → | |.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). Figure A5 demonstrates why this numerical saturation takes place.For values of Im[ ] becoming increasingly more negative, the sum in equation (A12) diverges more and more (panel d), owing to the combination of ultimately non-decreasing | | (panel a), and exponentially increasing | ( )| (panel b).We emphasize that numerical saturation is eventually unavoidable given that the 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, , would prove useful.It requires an enhanced accuracy on the ↦ → ( ) functions, i.e., on (i) the OrbitalElements.jlmappings 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, regularisation procedures should also be designed to push significantly further down in the complex plane, e.g., by reducing the effective range of summation over 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 DF (through its directional derivatives n• / J), and (iii) a bi-orthogonal basis.Some of these are available via OrbitalElements.jland AstroBasis.jl, but the user can easily supply their own.
For a given harmonic ℓ and for each resonance n and each matrix element ( , ), the calculations proceed in three phases.The first two aim at computing the ↦ → ( ) 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 .The third phase is to decompose these functions over Legendre polynomials, through the computation of the coefficients from equation (A11) using FiniteHilbertTransform.jl.Once these coefficients are known, the response matrix can be efficiently computed for any given complex frequency .

A4.1 Computing the G functions
Equation (1) reduces to equation (3) when taking where min n and max n are introduced in Appendix A1.3 and the J ↦ → (J) functions can be found in Appendix B. In both .In contrast to the solid lines which show the damped frequencies (i.e., < 0), the dashed lines -which show the unstable frequencies (i.e., > 0) -are always convergent.cases, they involve the Fourier transform of basis elements (Tremaine & Weinberg 1984).They read where the radius and the angles 1 and 2 − are implicit functions of the orbit, J, and the Hénon anomaly, , introduced in equation (A5).Following appendix B of Rozier et al. (2019), we perform the nested integration over in equation (A15) simultaneously, pushing the angles with Hilbert transforms.For these intense computations, julia provides an ideal language to optimise the running time while keeping a high readability.As an example, for a given (ℓ, n), a typical computation of all the needed values of { ℓn ( , )} , with "Fiducial" parameters and nradial = 100, takes ∼ 40 seconds on a single core of an M1 Macbook Air.Fortunately, parallelising over resonances is straightforward in julia.Several resonances can be then computed at once, depending on the thread count of the hardware.

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 3 spherical case (see section 4 in Hamilton et al. 2018).Then, simply using 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 2 = ℓ and the resonance n = (0, 0) does not contribute (to the ℓ = 0 response matrix) as in the spherical case.It also involves the Fourier transform of 2 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 ℓ = 1 damped mode calculations.The equations for the isochrone model potential and DFs can be found in FP22.

C1 ℓ = 2 modes -Radial Orbit Instability
The ROI in the isochrone model with a / c = 1.0, with 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 DF given in equation (G12) of Fouvry et al. (2021), we compute the growth rate for the fiducial a / c = 1.0 case.We find a fastest growing mode with M /Ω 0 [ a / c = 1.0] = 0.0236±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.jlwe 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 normalised anisotropy radius, a / c .A naive application of our libraries and the DF given in Fouvry et al. (2021) suggests that the isochrone stability boundary is a / 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 a / c decreases, regions of the ( , ) plane become unpopulated (see, e.g., figure 4 in Hamilton et al. 2018).Second, the gradient of the DF becomes extremely large at the = 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 DF hence considering the smoothed DF ˜ = exp(− / ) .We find here that different values of , the smoothing length, result in significantly different estimates for the growth rate, at fixed a / 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 a / 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 a / c .May & Binney (1986) estimated an instability threshold at the half-mass radius of a / c = 1.67.Merritt (1988) then identified an instability threshold of a / c = 1.9, from simulations using a harmonic code (Merritt 1987).Saha (1991) performed new simulations using the same harmonic code from MNRAS 000, 1-17 (0000) In general, the more resonances are considered, the more the complex frequency of the recovered (damped) ℓ = 1 mode converges to the frequency of the (neutral) translation mode M = 0+0i.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 a / c ≤ 4. Future work may return to the question of stability with improved numerical treatment.

C2 ℓ = 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 ℓ = 1 damped modes in King (1966) models using linear response theory.Weinberg (1994) also compared this prediction with perturbed -body simulations, finding qualitative agreement.Heggie et al. (2020) measured the same ℓ = 1 mode in -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 M /Ω 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 ℓ = 1 damped mode.There, we denote the set of parameters from FP22, for which we recover M /Ω 0 = 0.0143−0.00143i, in tight agreement with FP22.
However, we noted that varying the control parameters of LinearResponse.jlresults 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), the predicted complex frequency of the ℓ = 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 ℓ = 1 mode at 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 ℓ = 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 ( ) (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 = 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.

Figure 1 .
Figure 1.Left panel: Growth rate M = Im[ M ] of the radial orbit instability (ROI) as a function of the anisotropy radius 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 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.

Figure 2 .
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 M (see Table C2).Predictions for the Plummer model are shown in black.The dashed curve is d /d for the Plummer model.The predicted mode shapes are indistinguishable from d /d .
) 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

Figure 3 .
Figure 3. Left panel: Growth of the ℓ = 2 potential fluctuations energy ( ˜ ℓ , equation E3), through time in -body simulations of the radially-anisotropic Plummer sphere for four different anisotropy radii, a / c .Simulations with a / c = 0.95 and a / c = 1.00 are not shown, as instabilities were not reliably measured.Thin curves are individual simulation realisations; thick curves are the ensemble average for each value of a / c .Curves are aligned on the time of the first peak of the potential perturbation, denoted = 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 ℓ = 2 potential fluctuations measured from the -body simulations (solid curves).The dashed curves are the corresponding predictions from LinearResponse.jl.The linear predictions and numerical measurements are in pleasant agreement.

Figure 4 .
Figure 4. Left panel: Amplitude of the ℓ = 1 potential fluctuations energy ( ˜ ℓ , equation E3) vs. time in the isotropic -body realisations of Plummer spheres.Surprisingly, this time evolution could be compatible with a wave kinetic equation, as in equation (4).The best-fit solution is shown as a dashed gray curve.Right panel: Radial shape of ℓ = 1 density fluctuations measured from the -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 /d vs. 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.

Figure 5 .
Figure 5. Isocontours of the determinant of the ℓ = 2 (gravitational) dielectric matrix from equation (2) for the = 4 Zang disc.The dominant mode obtained by Zang (1976) is highlighted with a yellow cross and is recovered within 1% precision.

Figure 6 .
Figure 6.Shape of the ℓ = 2 dominant mode of the = 4 Zang disc as predicted by LinearResponse.jl(coloured dashed lines) overlaid with the shape obtained by Zang (1976) (in black, figure 9 therein).For both shapes, the contours denote the 10, 20, 40, 60 and 80% of the peak density, with only the overdensity being represented.The dotted circles show the corotation (CR) and inner Lindblad resonance (ILR) radii of the mode.Both linear predictions are in very good agreement.

•
the energy and angular momentum (equation 1.3 in Lynden-Bell 2015) = 2 apo ( apo )− 2 per ( per ) The figure shows a combination of the relative difference | ˜ − |/ and absolute difference | ˜ − | for each point in the grid, given by Δ

Figure A2 .
Figure A2.Accuracy of the numerical frequency calculation for the (analytical) isochrone model, using the distance Δ (equation A8).Extremely radial orbits are the most difficult ones to deal with.

Figure A3 .
Figure A3.Accuracy of the numerical inversion from ( , ) ↦ → ( , ) for the Plummer model.For each grid point in ( , ), we first compute the frequencies ( , ) , then "invert" these frequencies back to calculate ( ˜ , ˜ ).We define a relative accuracy metric Δ 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; semi-major axis is generally recovered to a higher precision.

Figure A5 .
Figure A5.Convergence of finite Hilbert Transform integral calculation for the ℓ = 2 radially-anisotropic Plummer cluster from Fig. 1, fixing a / c = 1.0, n = (−1, 2), = = 1, and using the Fiducial configuration settings.Here, the response matrix is evaluated for = 0+i , with < 0 (i.e., damped frequencies; solid curves) and zero oscillation frequency.Dashed curves show the corresponding calculation for − , i.e., unstable frequencies.Panel (a) shows that the coefficients .They decay for small , but they eventually saturate, as a result of finite numerical accuracy.Panel (b) shows that for all > 20, the quantity increases exponentially as a function of .Panel (c) shows the individual components of as a function of .For large enough, this product does not decay anymore.Finally, panel (d) shows the cumulative sum in equation (A12) as a function of .For sufficiently close to the real frequency line, the sum is convergent.However, as Im[ ] becomes more negative, the sum begins to diverge for smaller and smaller .In contrast to the solid lines which show the damped frequencies (i.e., < 0), the dashed lines -which show the unstable frequencies (i.e., > 0) -are always convergent.

Table 2 .
amounts to slightly modifying the boundary values of the resonance frequency n ∈ [ min n , max n ] and the variable ∈ [ − n , + n ].More precisely, the extrema of n (see equation B7 in FP22) are now reached either in the four edges ( min , 1 2 ), ( min , c [ min ]), ( max