Gravitation – Instabilities – Galaxies: kinematics and dynamics – Software: public release

1 Introduction↩︎

The response of self-gravitating systems to perturbations has been studied for more than half a century, mostly using the matrix method [1]. As discussed in detail in section 5.3.2 of [2], this method uses basis functions to represent the gravitational response of the stellar system. The response matrix \(\boldsymbol{\mathsf{M}}\) at a given complex frequency \(\omega\) is generically \[\label{eq:ResponseMatrix} \boldsymbol{\mathsf{M}}(\omega) = \sum_{\mathbf{n}}\! \int\limits_{\cal L} \!\!\mathrm{d}\mathbf{J}\frac{\boldsymbol{\mathsf{G}}_{\mathbf{n}}(\mathbf{J})}{\mathbf{n}\!\cdot\!\mathbf{\Omega}\left(\mathbf{J}\right)-\omega}.\tag{1}\] This equation encodes the physics of (linear) stability in self-gravitating systems. Here, \(\sum_{\boldsymbol{n}}\) is the sum over all (allowed) resonance vectors \({\mathbf{n}}\); \(\int_{\cal L} d{\boldsymbol{J}}\) is the scan over the populated orbital (action) space (with an appropriate Landau prescription); \({\mathbf{n}\!\cdot\!\mathbf{\Omega}\left(\mathbf{J}\right)\!-\!\omega}\) is the resonant amplification at a given (complex) frequency \(\omega\) and for orbital frequencies \(\mathbf{\Omega}\); \({\boldsymbol{\mathsf{G}}_{\mathbf{n}}}\) are functions of the system’s DFand are matrices by virtue of decomposing the Newtonian pairwise interaction on potential basis elements that encode the long-range nature of the gravitational interaction [3].

The system sustains a mode at a complex frequency \(\omega\) when the condition \[\label{eq:def95dielectric} \det[\boldsymbol{\mathsf{\varepsilon}}(\omega)]=0, \quad \mathrm{with} \; \boldsymbol{\mathsf{\varepsilon}}(\omega) = \boldsymbol{\mathsf{I}}- \boldsymbol{\mathsf{M}}(\omega),\tag{2}\] is met, where \(\boldsymbol{\mathsf{I}}\) is the identity matrix and \(\boldsymbol{\mathsf{\varepsilon}}\) is called the (gravitational) dielectric matrix. We write frequencies satisfying this equality as \(\omega_{\mathrm{M}}\), which we further break down into \({\omega_{\mathrm{M}}\!=\! \Omega_{\mathrm{M}}\!+\!\mathrm{i}\gamma_{\mathrm{M}}}\), with the oscillation frequency, \(\Omega_{\mathrm{M}}\), and the growth rate, \(\gamma_{\mathrm{M}}\). Modes may appear at any point in the complex plane. Modes with \({\gamma_{\mathrm{M}}\!>\!0}\) are unstable, modes with \({\gamma_{\mathrm{M}}\!=\!0}\) are neutral, and modes with \({\gamma_{\mathrm{M}}\!<\!0}\) are damped.

Previous studies in linear response have probed modal structure in self-gravitating discs [4]–[11] and spheres [12]–[20]. Most calculations have been interested in studying instabilities, i.e., modes with \({\gamma_{\mathrm{M}}\!>\! 0}\) [21]. Traditionally, finding damped modes has been challenging because of the analytical continuation required by Landau’s prescription. [4] presented a first prediction in razor-thin discs while relying heavily on the scale invariance of the Mestel potential. [22] 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 [19]. Damped modes are interesting for their capacity to dissipate energy through interactions of the mode and resonating stellar orbits [23]. This phenomenon may have implications for real astrophysical systems, including accelerating the overall relaxation [24], [25], as well as the secular dynamics captured by the Balescu–Lenard equation [26], [27], and the so-called mode-particle interactions [28].

In this paper, building upon the method developed in [19], hereafter , we introduce software libraries in the julialanguage to perform linear response calculations. Details of the toolbox are given in Appendix 6. 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 julialibraries 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 [4]. 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 \(\ell\) (historically denoted \(m\) for the discs). For these systems and within the appropriate “resonance” coordinate system \((u,v)\) introduced in equation (10) of , each element of the response matrix of harmonic \(\ell\) reads \[\label{eq:M95int95u} M^{\ell}_{pq}(\omega) = \sum_{\mathbf{n}\in\mathbb{Z}^2} {\int_{\substack{ \\[0.83ex] -1 \\\mathcal{L}}}^{1}}\!\! \mathrm{d}u \frac{G^{\ell\mathbf{n}}_{pq} (u)}{u - \varpi_{\mathbf{n}} (\omega)} ,\tag{3}\] where the functions¹ \(G^{\ell\mathbf{n}}_{pq}(u)\) depend on the Fourier transform of bi-orthogonal basis elements (indexed by \(p\) and \(q\)), the resonance number, \({\mathbf{n}\!=\!(n_1,n_2)}\), and involve an integral over some other coordinate \({ \mathrm{d}v }\), as explained in Appendix 6.4.1. These functions, which depend on the dimensionality of the problem, are detailed in Appendix 7. In equation 3 , \({\varpi_{\mathbf{n}} (\omega)}\) stands for the rescaled (complex) frequency introduced in equation (11) of . Importantly, the imaginary parts of the rescaled frequency \(\varpi_{\mathbf{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 . The goal behind producing LinearResponse.jlis to extend the prescription used for the isochrone sphere in  to any numerically-given potential and DF, as well as to discs. To be as versatile as possible, the software is decomposed in four libraries.

• OrbitalElements.jl(Appendix 6.1), which, given any central potential and its two first derivatives, performs various changes of orbital coordinates to compute the actions \((J_{\mathrm{r}},L)\), orbital frequencies \((\Omega_1,\Omega_2)\) and resonance coordinates \((u,v)\).²

• AstroBasis.jl(Appendix 6.2) 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.

• FiniteHilbertTransform.jl(Appendix 6.3) performs the (Landau’s prescription-compliant) finite Hilbert transform using Legendre polynomials, as detailed in appendix D of .

• LinearResponse.jl(Appendix 6.4) is the driver library. Relying on the three (independent) previous libraries, LinearResponse.jlimplements the computation of the response matrix \({\boldsymbol{\mathsf{M}}^{\ell}\!(\omega)}\). It performs the Fourier transform of the basis elements, then computes the functions \(G^{\ell\mathbf{n}}_{pq}(\mathbf{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 , we consider the isochrone potential in Appendix 8. We use its analytical expressions to check our numerical calculation of the orbital elements, and recover all of ’s results. This is reassuring. We also point out two caveats raised by our exploration, namely (i) diverging gradients of the anisotropic DFat the edge of the domain leading to convergence issues in the \({\ell\!=\!2}\) case, and (ii) a (neutral) translation mode that becomes damped in the \({\ell\!=\!1}\) case.

In this section, we expand our linear response study to the [29] model. In particular, we first study the (\({\ell\!=\!2}\)) ROI for radially anisotropic models (Section 3.1), before investigating isotropic DFsand \({\ell\!=\!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 9.

3.1 \({\ell\!=\!2}\) modes – Radial Orbit Instability↩︎

The ROI arises in models with some degree of radial anisotropy, i.e., the DFcontains an enhancement at low angular momentum values [30]. Often, this is accomplished via the introduction of some anisotropy radius, \(r_{\mathrm{a}}\), outside of which the fraction of radial orbits is enhanced. The mechanisms for the instability are discussed in [21], [31]. The classic ROI results in a quadrupole (\({\ell\!=\!2}\)) zero oscillation frequency mode (\({ \Omega_{\mathrm{M}}\!=\! 0 }\)) that grows exponentially in time (\({\gamma_{\mathrm{M}}\!>\!0}\)). In practice, the ROI is generated by the ILR, \(\mathbf{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. [32] estimated the threshold for stability in the Plummer model to be \({r_{\mathrm{a}}/b_{\mathrm{c}}\!\simeq\! 0.9}\) using the criteria of [33], with \(b_{\mathrm{c}}\) the Plummer’s scale radius (equation 28 ). Using the mean-field code from [34], [32] also reported that the stability threshold was somewhere around \({r_{\mathrm{a}}/b_{\mathrm{c}}\!\simeq\!1.1}\). Relying on direct \(N\)-body simulations, [35] confirmed that \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!0.75}\) was unstable, but, in contrast to [32], [35] found that \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!1.0}\) was stable.

We test our linear response machinery on a range of \(r_{\mathrm{a}}/b_{\mathrm{c}}\) values. Fig. 1 shows the results of searching for modes using a “Fiducial” set of control parameters.³ We find a threshold for instability of \({r_{\mathrm{a}}/b_{\mathrm{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 [36] 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 6, 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 \(\ell\) harmonics are decoupled in linear response calculations, we investigate \({\ell\!=\!(1,3,4)}\) for similar instability behaviour. We do not find evidence for any modes for \({\ell\!\ne\!2}\).

3.2 A search for \({\ell\!=\!1}\) damped modes↩︎

We next search for evidence of damped \({\ell\!=\!1}\) modes, akin to the findings of [22] and [25] for King models. Using an isotropic DF(see Appendix 9), we perform a search in the complex plane. We find zeros in the determinant of the dielectric matrix \(\boldsymbol{\mathsf{\varepsilon}}\), but these appear to be spurious non-converged realisations of the neutral translation mode. The right column of Table 7 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 8.2). In addition, the “false” mode clearly resembles the density derivative \({\mathrm{d}\rho/\mathrm{d}r}\), 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 \({\ell\!=\!1}\) modes. For completeness, we also search other harmonics \({\ell\!=\!(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 \({\mathrm{d}\rho/\mathrm{d}r}\). As expected given the results of , the predicted mode appears to be the neutral \({\omega_{\mathrm{M}}\!=\!0\!+\!0\mathrm{i}}\) 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.

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 9): one radially anisotropic, and one isotropic. We realise the \(N\)-body simulations using initial conditions from PlummerPlus [35]. We run the models using EXP[37], [38], a basis function expansion \(N\)-body code. Further numerical details for both the initial conditions and the \(N\)-body integration are given in Appendix 10.

3.3.1 Anisotropic simulations, Radial Orbit Instability↩︎

To attempt to recover the ROI in the Plummer clusters, we select six different values of \(r_{\mathrm{a}}/ b_{\mathrm{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_{\mathrm{a}}/b_{\mathrm{c}}\), we run an ensemble of six simulations of \(2^{18}\) particles each.

Fig. 3 shows the results for the \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\![0.75,0.80,0.85,0.9]}\) runs using EXP(see Appendix 10). In general, the runs for any value of \(r_{\mathrm{a}}/b_{\mathrm{c}}\) show similar behaviour, though the dispersion between realisations increases with \(r_{\mathrm{a}}/b_{\mathrm{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 \(r_{\mathrm{a}}/b_{\mathrm{c}}\) models; the qualitative agreement with the predicted mode shapes and trends (shown as dashed curves) with \(r_{\mathrm{a}}/b_{\mathrm{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_{\mathrm{a}}/b_{\mathrm{c}}\!=\!0.95}\), only three of the six runs performed appear to show any instability. Thus, \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!0.95}\) appears to be a weakly unstable model. We do not measure an instability for any of the \({r_{\mathrm{a}}/b_{\mathrm{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 within \({3\sigma}\) of the predicted growth rate, and both approaches recover the qualitative trend of decreasing \(\gamma_{\mathrm{M}}\) with increasing \(r_{\mathrm{a}}/b_{\mathrm{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 \({r_{\mathrm{b}}\!=\!1}\), match the basis used to run the simulation parameters (Appendix 10). This ensures a straightforward comparison with the \(N\)-body simulations. We find that the predicted eigenvectors (equation 25 ) 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 EXPto 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 DF(see Appendix 10). The left panel of Fig. 4 shows the run with time of the (rescaled) energy, \({\tilde{E}_{1}\!\propto\!E_{\ell=1} }\) (Eq. 31 ), for each simulation. Interestingly, \(\tilde{E}_{1}\) increases and slowly approaches a saturated level that differs between realisations. Upon inspection, the amplitude \(\tilde{E}_{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 \({\mathrm{d}\rho/\mathrm{d}r}\), the results from the \(N\)-body simulation are generally not consistent with \({\mathrm{d}\rho/\mathrm{d}r}\). 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 \({\mathrm{d}\rho/\mathrm{d}r}\). 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. [39], equation (21) therein, put forward a wave kinetic equation [28] that describes the (linear) saturation of potential fluctuations in self-gravitating systems. It reads \[\label{eq:wave95kinetic} \frac{\mathrm{d}\tilde{E}_{\ell}}{\mathrm{d}t} = 2\gamma_{\mathrm{M}}\left(\tilde{E}_{\ell} - \tilde{E}_{\mathrm{th}}\right),\tag{4}\] where \({\tilde{E}_{\ell}}\) is the (rescaled) energy in the \({ \ell \!=\! 1 }\) gravitational fluctuations and \(\tilde{E}_{\mathrm{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_{\mathrm{M}}\), transients only fade away after a (few) \({ 1 / \gamma_{\mathrm{M}}}\), at which stage fluctuations are fully dressed by self-gravity [40]. 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}_{\mathrm{th}}\), \(\gamma_{\mathrm{M}}\) and the initial energy fluctuations, \({\tilde{E}_{\mathrm{init}}\!=\! \tilde{E}_{\ell} (t\!=\!0)}\). We find \({\tilde{E}_{\mathrm{init}}\!=\!1.0\!\times\!10^{-6}}\), \({\tilde{E}_{\mathrm{th}}\!=\!8.0\!\times\!10^{-3}}\), and \({\gamma_{\mathrm{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_{\mathrm{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}_{\mathrm{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 (\({2D}\)) discs. In that case, stars still interact through the usual Newtonian interaction potential \({U(\mathbf{r},\mathbf{r^{\prime}})\!=\!-G/|\mathbf{r}-\mathbf{r^{\prime}}|}\), with \(G\) 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 \((J_{\mathrm{r}},L)\), just like in the case of spheres. Hence, all the previous numerical methods from naturally apply to razor-thin discs. The only required modifications are: (i) marginally modifying the linear response integrand, \({G^{\ell\mathbf{n}}_{pq}}\), in equation 3 , given in Appendix 7.2, and (ii) implementing \({2D}\) bi-orthogonal basis elements [3], [41], given in Appendix 11.1.

In order to validate our implementation, we aim to recover well-documented unstable modes in razor-thin discs. In practice, following the work of [4], we consider Mestel discs [42], 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 [8], the DFof 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⁴: the sharper this inner cut-out, the stronger the instability [4]. 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 [43]. We refer to Appendix 11.2 for details and notably the definition of the scaling frequency \(\Omega_0\). Note that the method based on finite Hilbert transform implemented in LinearResponse.jlis 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 DFtaper (see Appendix 6.1.3 for details). It would be rewarding to extend ’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.⁵ As reported in Table 2, we find a satisfying agreement between the (semi-analytical) predictions of [4] and [43], 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 [44].

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{\mathsf{\varepsilon}}(\omega) }\) from equation 2 , through its level contours in the complex-frequency plane. This determinant vanishes at the frequency \({\omega_{\mathrm{M}}/\Omega_0\!=\!0.878\!+\!0.126\mathrm{i}}\), 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 [4]. 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 6.3.

As a concluding remark, let us note that here we used the generic basis from [41], 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.jlto 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 juliato 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 \({\ell\!=\!1}\) response, we confirm that the results of  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 radial orbit instability, 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 DF, 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 ’s work to the stability analysis of razor-thin discs is convincingly validated against known semi-analytical results in Mestel discs [4], [43]. In particular, we satisfyingly retrieve both the frequency and shape of the \({\ell\!=\!2}\) growing mode from [4]. 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 [45], 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 [22] ; (v) building the self-gravitating amplification kernel in the time domain rather than the frequency domain [46], [47].

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 [48] 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 [24], or to explore their adiabatic response as one slowly varies the mean field or the stellar cluster’s DF [49].

Data Availability↩︎

The julialibraries developed for this paper, LinearResponse.jl, OrbitalElements.jl, AstroBasis.jland FiniteHilbertTransform.jl, are freely available on GitHub. The EXPsimulations used for comparison are available upon reasonable request to the corresponding author.

Acknowledgements↩︎

We thank S. Rozier, M. Weinberg, C. Hamilton, S. de Rijcke, J. Peñarrubia, and A. Naik for insightful conversations. 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. 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 julialanguage [50]. We used scipy [51] to fit the wave kinetic equation solution. For visualisations, we used iPython [52], numpy [53], and matplotlib [54].

6 Linear Response Software↩︎

6.1 OrbitalElements.jl: Orbit Frequency Computation↩︎

The julialibrary OrbitalElements.jlis optimised 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 [55], gala [56] or AGAMA [57], but rather design our own specific library tailored for our needs. At the heart of OrbitalElements.jlis 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 \((r_{\mathrm{per}},r_{\mathrm{apo}})\);

• the effective semi-major axis and eccentricity \[\label{eq:AE} a = \frac{r_{\mathrm{per}}+r_{\mathrm{apo}}}{2} , \quad e = \frac{r_{\mathrm{apo}}-r_{\mathrm{per}}}{r_{\mathrm{apo}}+r_{\mathrm{per}}};\tag{5}\]

• the energy and angular momentum [58] \[\label{eq:EL} E \!=\! \frac{r_{\mathrm{apo}}^2 \psi(r_{\mathrm{apo}}) \!-\! r_{\mathrm{per}}^2\psi(r_{\mathrm{per}})}{r_{\mathrm{apo}}^2 \!-\! r_{\mathrm{per}}^2} , \, L \!=\! \bigg[ \frac{2(\psi(r_{\mathrm{apo}}) \!-\! \psi(r_{\mathrm{per}}))}{r_{\mathrm{per}}^{-2} \!-\! r_{\mathrm{apo}}^{-2}} \bigg]^{1/2}\!\!\! ;\tag{6}\]

• the actions \((J_r,L)\), where \[\label{eq:Jr} J_r = \displaystyle\frac{1}{\pi}\int_{r_{\mathrm{per}}}^{r_{\mathrm{apo}}} \!\! \mathrm{d}r \, v_r,\tag{7}\] with \(J_{r}\) the radial action and \({v_r\!=\!\sqrt{2(E\!-\!\psi(r))\!-\!L^2/r^2}}\) the instantaneous radial velocity;

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

• the frequency ratios \((\alpha,\beta)\) from which the frequencies are computed [59] \[\label{eq:AlphaBeta} \frac{1}{\alpha} = \frac{\Omega_0}{\Omega_1} = \displaystyle \frac{1}{\pi}\int_{r_{\mathrm{per}}}^{r_{\mathrm{apo}}} \frac{\mathrm{d}r}{ v_r}, \quad \beta = \frac{\Omega_2}{\Omega_1} =\displaystyle \frac{L}{\pi}\int_{r_{\mathrm{per}}}^{r_{\mathrm{apo}}} \frac{\mathrm{d}r}{r^2 v_r},\tag{8}\] with \(\Omega_{0}\) some given frequency scale (typically the central radial frequency in cored profiles);

• the resonance-specific (i.e., dependent on \(\mathbf{n}\)) coordinates \((u,v)\) from (appendix B therein).

In practice, OrbitalElements.jlis centred around the effective semi-major axis and eccentricity \((a,e)\), 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 , 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.1.1 Forward mappings↩︎

The direct calculations compute orbital elements (frequency ratios, radial action, etc.) through integrals of the form \({ \int_{r_{\mathrm{per}}}^{r_{\mathrm{apo}}} \! \mathrm{d}r ... }\) However, due to the diverging integrand \({1/v_r}\) in equations 8 , the radius \(r\) is not an appropriate integration variable. To improve numerical accuracy, we cure divergences using a mapping from equation 51 in [60]. It reads \[\label{eq:Henon95mapping} r(w) = a\left[1+ef(w)\right]; \; f(w) = w\left(\tfrac{3}{2} - \tfrac{1}{2}w^2\right),\tag{9}\] 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{\mathrm{d}r / \mathrm{d}w}{v_r(w)}, \label{eq:Theta}\tag{10}\] is no longer divergent and equations 8 simply read \[\tag{11} \begin{align} \frac{1}{\alpha} &= \frac{1}{\pi}\int_{-1}^1 \mathrm{d}w \, \Theta(w), \tag{12} \\ \beta &= \frac{L}{\pi}\int_{-1}^1 \mathrm{d}w \, \frac{\Theta(w)}{r^2(w)}. \tag{13} \end{align}\]

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 [60]. The radial action, \(J_{r}\), is readily and easily computed by integrating equation 7 over \(w\). In Fig. 7, 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.⁷

For near-circular orbits, i.e., \({e\!\to\!0}\), and/or very small semi-major 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).

6.1.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)\!\mapsto\!(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 da and de. Finally, the Newton–Raphson algorithm is controlled with a maximal number of iterations (ITERMAX) and a accuracy goal (inv\(\varepsilon\)). The effective accuracy of the inversion is mainly determined by the orbital integration (NINT) and the finite differences (da and de).

6.1.3 Resonances variables↩︎

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

In OrbitalElements.jl, we extend the resonance-specific mapping to situations where the potential and DFare 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 11.2), where the frequency profile diverges in the “unpopulated” center.

We therefore truncate the domain in \({(\alpha,\beta)}\) effectively explored. This does not affect the computation of the frequencies, \(\mathbf{\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_{\mathrm{min}}\!=\!\alpha_{\mathrm{c}}(\texttt{rmax})}\) and \({\alpha_{\mathrm{max}}\!=\!\alpha_{\mathrm{c}}(\texttt{rmin})}\) with \({ \alpha_{\mathrm{c}}(r) \!=\! \alpha (a \!=\! r , e \!=\! 0) }\) the (outward decreasing) circular frequency ratio. For example, in the case of a Mestel disc (Appendix 11.2), \({\texttt{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 , adding these truncations only amounts to slightly modifying the boundary values of the resonance frequency \({\omega_{\mathbf{n}}\!\in\![\omega_{\mathbf{n}}^{\mathrm{min}},\omega_{\mathbf{n}}^{\mathrm{max}}]}\) and the variable \({v\!\in\![v_{\mathbf{n}}^{-},v_{\mathbf{n}}^{+}]}\). More precisely, the extrema of \({\omega_{\mathbf{n}}}\) are now reached either in the four edges \({(\alpha_{\mathrm{min}},\tfrac{1}{2})}\), \({(\alpha_{\mathrm{min}},\beta_{\mathrm{c}}[\alpha_{\mathrm{min}}])}\), \({(\alpha_{\mathrm{max}},\tfrac{1}{2})}\), \({(\alpha_{\mathrm{max}},\beta_{\mathrm{c}}[\alpha_{\mathrm{max}}])}\) or along the curve \({(\alpha,\beta_{\mathrm{c}}[\alpha])}\) with \({\alpha_{\mathrm{min}}\!\leq\!\alpha\!\leq\!\alpha_{\mathrm{max}}}\). For the \(v\) variable, the first two constraints of equation (B10) in become \({\alpha_{\mathrm{min}}\!\leq\!v; \, v\!\leq\!\alpha_{\mathrm{max}}}\). Finally, the same procedure as in is used to obtain the boundary values.

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

6.1.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 3, but are described in more detail throughout the previous subsections.

6.1.5 Tests↩︎

Figs. 8 and 9 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. 8, 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 \[\label{eq:def95Delta} \Delta_{\alpha\beta} = \sqrt{\left( |\tilde{\alpha}-\alpha| / \alpha \right)^2 + \left( |\tilde{\beta}-\beta | / \beta \right)^2}.\tag{14}\] In practice, \(\Delta_{\alpha\beta}\) is dominated by the accuracy contribution from \(\tilde{\beta}\). To probe the full range of relevant semi-major axes, we sample \(a\) using 1000 log-spaced points \({a/b_{\mathrm{c}}\!\in\![10^{-3},10^3]}\), and 1000 linear-spaced points \({e\!\in\![0,1]}\). The accuracy is better than 1% at all points sampled, and in most cases is better than 0.001% accurate.

In Fig. 9, 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{\left(|\tilde{a}-a| / a \right)^2 + \left(|\tilde{e}-e|\right)^2}. \label{eq:inversionaccuracyae}\tag{15}\] 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 \((a,e)\) values, namely orbits close to radial or circular, and orbits with small semi-major 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 \((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. 9 took \({\sim\!140\mathrm{s}}\) on a single core of an M1 Macbook Air.

6.2 AstroBasis.jl: Basis Computation↩︎

Fundamental to the matrix method are the chosen basis functions. The AstroBasis.jllibrary is an implementation of several different radial basis functions, \({r\!\mapsto\!U^{\ell}_p(r)}\), with a straightforward interface.⁸ At present, AstroBasis.jlsupports the bases from [41], [36], [3], [61]/[13] (Bessel), and [62]. The spherical calculations in this paper use the basis from [36], which is re-described in appendix B of [63]. The associated lowest-order function matches the Plummer potential exactly. The disc calculations in this paper use the basis from [41].

The user must choose parameters to define the basis, including the gravitational constant \(G\) (G), the scaling radius for the basis \(r_{\mathrm{b}}\) (rb), the maximum harmonic number \({\ell_{\mathrm{max}}}\) (lmax), and the maximum radial number \(n_{\rm max}\) (nradial). Basic parameters are listed in Table 4. Some bases may require additional parameters.

3pt

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 10 shows an example of basis functions from AstroBasis.jl. In this case, we compare the shape of the first \(\ell\!=\!2\) element of [36]’s basis to the shape of the ROI modes from the isochrone cluster (Appendix 8.1) 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 [38], [64]. The applicability of these tailored basis functions to matrix method calculations should be explored further in future work.

6.3 FiniteHilbertTransform.jl: Specialised Legendre Integration↩︎

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_{\substack{ \\[0.83ex] -1 \\\mathcal{L}}}^{1}}\mathrm{d}u \frac{G (u)}{u - \omega}, \label{eq:FHT}\tag{16}\] 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 \(\mathbf{n}\) and basis element pairs \({ (p,q) }\), as highlighted in equation 3 .

To compute equation 16 , we first approximate its numerator via \[G(u) \!\simeq\! \sum_{k = 0}^{K_{u}-1} a_{k} \, P_{k} (u) , \label{eq:approx95Gu}\tag{17}\] where, importantly, \(K_u\) (defined as Ku in FiniteHilbertTransform.jl) controls the quality of the approximation. In practice, the coefficients \(a_{k}\) entering equation 17 are estimated through a Gauss–Legendre quadrature with \(K_{u}\) nodes (see equation D3 in ).¹⁰

Then, the integral from equation 16 is computed through \[M (\omega) = \sum_{k = 0}^{K_u - 1} a_{k} D_{k} (\omega) , \label{eq:Msum}\tag{18}\] where \[D_{k} (\omega) = {\int_{\substack{ \\[0.83ex] -1 \\\mathcal{L}}}^{1}}\mathrm{d}u \frac{P_{k} (u)}{u - \omega}. \label{def95D}\tag{19}\] The function \(D_k(\omega)\) explicitly implements Landau’s prescription and is defined in appendix D of . In brief, \(D_k(\omega)\) is proportional to Legendre functions of the second kind when \(\mathrm{Im}[\omega]>0\), with an added contribution when \(\mathrm{Im}[\omega]\le0\) [2]. The library FiniteHilbertTransform.jlimplements the computation of (i) the \(a_k\) coefficients for a given function \(G\), and (ii) the \({\omega\!\mapsto\!D_k(\omega)}\) functions from equation 18 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 [67], for damped frequencies, i.e., negative \(\mathrm{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 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). Figure 11 demonstrates why this numerical saturation takes place. For values of \(\mathrm{Im}[\omega]\) becoming increasingly more negative, the sum in equation 18 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 18 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.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 \(k\) in equation 18 .

6.4 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 \({\mathbf{n}\!\cdot\!\partial F/\partial\mathbf{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 \(\ell\) and for each resonance \(\mathbf{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 17 using FiniteHilbertTransform.jl. Once these coefficients are known, the response matrix can be efficiently computed for any given complex frequency \(\omega\).

6.4.1 Computing the G functions↩︎

Equation 1 reduces to equation 3 when taking \[\label{eq:G95u95J} G(u) = \int_{0}^{1} \mathrm{d}v' \frac{2}{\Omega_0 (\omega_{\mathbf{n}}^{\mathrm{max}}- \omega_{\mathbf{n}}^{\mathrm{min}})}\left| \frac{\partial\mathbf{J}}{\partial(u,v')} \right| G(\mathbf{J}[u,v']),\tag{20}\] where \(\omega_{\mathbf{n}}^{\mathrm{min}}\) and \(\omega_{\mathbf{n}}^{\mathrm{max}}\) are introduced in Appendix 6.1.3 and the \({\mathbf{J}\!\mapsto\!G(\mathbf{J})}\) functions can be found in Appendix 7. In both cases, they involve the Fourier transform of basis elements [59]. They read \[\label{eq:Wln} W^{\ell\mathbf{n}}_p(\mathbf{J}) = \int_{-1}^{1} \mathrm{d}w \, \frac{\mathrm{d}\theta_1}{\mathrm{d}w}\, U^{\ell}_p (r) \, \cos\!\big(n_1 \theta_1 \!+\! n_2[\theta_2\!-\!\phi]\big),\tag{21}\] 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 9 . Following appendix B of [18], we perform the nested integration over \(w\) in equation 21 simultaneously, pushing the angles with \[\frac{\mathrm{d}}{\mathrm{d}w} (\theta_1, \theta_2\!-\!\phi) = \left(\Omega_1,\Omega_2\!-\! L / r^{2} \right)\Theta(w) , \label{grad95theta}\tag{22}\] where \(\Theta\) follows from equation 10 . In practice, we use the RK4 scheme [68] 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\!\rightarrow\!0\)) typically require more integration steps. We therefore add the ability to choose the number of steps adaptively, following \[\label{eq:ADAPTIVEKW} \texttt{Kw} \leftarrow \lceil \texttt{Kw}/(0.1+(1-e))\rceil.\tag{23}\] This option is set by the boolean parameter ADAPTIVEKW.

Computing the Fourier transforms at each point \((u,v')\) and for each resonance \(\mathbf{n}\) is the most intensive part of the linear response computations. For the Fiducial configuration settings, each resonance takes approximately 16 seconds on an M1 Macbook Air. The calculation of the resonances can be parallelised. Helpfully, these quantities do not depend on the DF 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 20 is performed on a rescaled variable \(v'\), compared to the choice \({v\!=\!\alpha}\) (for \({n_2\!\neq\!0}\)) used in . Indeed, in practice, we use \[v'=\left[\frac{v-v_{\mathbf{n}}^{-}}{v_{\mathbf{n}}^{+}-v_{\mathbf{n}}^{-}}\right]^\texttt{VMAPN}, \label{eq:VMAPN}\tag{24}\] with \({v\!\in\![v_{\mathbf{n}}^{-},v_{\mathbf{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 midpoint rule with Kv points.

In practice, each function \({u\!\mapsto\!G^{\ell\mathbf{n}}_{pq}(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 \(\mathbf{n}\), we end up with \({\texttt{Ku}\!\times\!\texttt{Kv}}\) vectors \(\{W^{\ell\mathbf{n}}_p(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 \(\mathbf{n}\), in arrays of size \({\texttt{Ku}\!\times\!\texttt{Kv}\!\times\!\texttt{nradial}}\). For the \(G\) functions, for each resonance \(\mathbf{n}\), we end up with \({\texttt{nradial}\!\times\!\texttt{nradial}}\) vectors \({\{{G}^{\ell \mathbf{n}}_{pq}(u_k)}\}_k\) of length Ku. These are also cached.

6.4.2 Computing the response matrix↩︎

Once the \(G\) functions tabulated, their Legendre decomposition coefficients from equation 17 are computed using FiniteHilbertTransform.jl. For each resonance \(\mathbf{n}\), this results in Ku matrices of size \({\texttt{nradial}\!\times\!\texttt{nradial}}\).¹¹

Then, for a given complex frequency \(\omega\), the response matrix from equation 3 is computed using equation (18 ). We write \[\boldsymbol{\mathsf{M}}^{\ell}(\omega) = \sum_{\mathbf{n}} \sum_{k} D_k(\varpi_{\mathbf{n}})\,\boldsymbol{\mathsf{A}}^{\ell\mathbf{n}}_k ,\] where \({\omega\!\mapsto\!\varpi_{\mathbf{n}}(\omega)}\) is given by equation (11) in , and the \(D_k\) functions are computed by via FiniteHilbertTransform.jl.

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

1. 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.

2. Identifying the frequency of a mode, \(\omega_{\mathrm{M}}\), using a gradient descent algorithm.¹² This is the method used in the main text to report on the modes’ frequencies.

3. Identifying the shape of a mode, \(\psi_\mathrm{M}\), given its frequency. We write \[\label{eq:def95Xp} \psi_{\mathrm{M}}(\mathbf{r}) = \mathrm{Re}\bigg[ \sum_{p} X_p \, \psi_p(\mathbf{r}) \bigg] ,\tag{25}\] where \({ \{ X_p \}_{p} }\) is the kernel eigenvector of the dielectric matrix, and \({ \psi_p (\mathbf{r})\!\propto\!U^{\ell}_{p} (r) }\) (with an angular dependence set by the considered geometry) runs over the potential basis functions used in the response matrix computation.

6.4.3 Linear response parameters↩︎

As for OrbitalElements.jl, control parameters for LinearResponse.jlare set in a juliastructure taken as an optional argument of all exposed functions. In addition, all control parameters from OrbitalElements.jlare passed to LinearResponse.jlas 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 summarised in Table 5.

The accuracy of LinearResponse.jl is not straightforward to predict. Let us nonetheless report on some heuristics. First, increasing the number \(\texttt{nradial}\) of basis elements leads to more oscillating radial functions, hence more difficult reconstructions. This requires therefore a more precise quadrature. The same with \(\texttt{n1max}\), 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, juliaprovides an ideal language to optimise 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^{\ell\mathbf{n}}_{p}(u,v)\}_{p}}\), with “Fiducial” parameters and \({\texttt{nradial}\!=\!100}\), takes \({\sim\!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.

7 Linear response integrands↩︎

7.1 For spheres↩︎

As in , in \({3D}\) the full expression of \(\boldsymbol{\mathsf{G}}^{\ell\mathbf{n}}(\mathbf{J})\) in equation 1 is given by \[\label{eq:def95Gpq} G^{\ell \mathbf{n}}_{pq} (\mathbf{J}) \!=\! - \frac{2 (2 \pi)^{3}}{2 \ell + 1} \big| y_{\ell}^{n_{2}} \big|^{2} L\, \left(\mathbf{n}\!\cdot\! \frac{\partial F}{\partial\mathbf{J}}\right) W_{p}^{\ell\mathbf{n}} (\mathbf{J}) W_{q}^{\ell\mathbf{n}} (\mathbf{J}) ,\tag{26}\] where we introduced \({y_{\ell}^{n}\!=\!Y_{\ell}^{n}(\tfrac{\pi}{2},0)}\), with \({Y_{\ell}^{m}(\vartheta,\phi)}\) the spherical harmonics normalised so that \({\!\int\! \mathrm{d}\vartheta \mathrm{d}\phi \sin\!(\vartheta) |Y_{\ell}^{m}(\vartheta,\phi)|^{2} \!=\! 1}\). Further details may be found in , including the Jacobian for the conversion to \(\boldsymbol{\mathsf{G}}^{\ell \mathbf{n}}(u,v)\), needed in practice by LinearResponse.jl. In equation 26 , the Fourier-transformed basis elements follow from equation 21 .

7.2 For discs↩︎

It is straightforward to extend the linear response computations of  to the case of razor-thin discs with axial symmetry. In this case, the definition of the response matrix [2] immediately applies, without any integral over a third dummy action as in the \({3D}\) spherical case [24]. Then, simply using \[\label{eq:def95Gpq95discs} G_{pq}^{\ell\mathbf{n}}(\mathbf{J}) = - (2 \pi)^2 \,\delta_{\ell}^{n_2}\, \left( \mathbf{n}\!\cdot\! \frac{\partial F}{\partial\mathbf{J}} \right) W_{p}^{\ell\mathbf{n}} (\mathbf{J}) \,W_{q}^{\ell\mathbf{n}} (\mathbf{J}),\tag{27}\] instead of equation 26 , one can compute the linear response of razor-thin discs using all the numerical tools from LinearResponse.jl. Equation 27 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 21 , which are described in Appendix 11.1.

8 Validation with the isochrone sphere↩︎

In this section, we revisit the calculations from  for the spherical isochrone potential [69]. Because the isochrone has analytic expressions for the frequencies, one can readily evaluate integrals to high numerical precision, as used in . However, here we exercise the empirical frequency computation from OrbitalElements.jl(Appendix 6.1) 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 DFscan be found in .

8.1 \({\ell\!=\!2}\) modes – Radial Orbit Instability↩︎

The ROI in the isochrone model with \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!1.0}\), with \(b_{\mathrm{c}}\) the isochrone scale radius, has become a regular test for different aspects of the matrix method implementation [24], [46], [63], [64]. Using the radially-anisotropic DFgiven in equation (G12) of [63], we compute the growth rate for the fiducial \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!1.0}\) case. We find a fastest growing mode with \({\gamma_{\mathrm{M}}/\Omega_{0}[r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!1.0]\!=\!0.0236\!\pm\!0.0015}\). This matches the result of  at the 2 per cent level, and [15] at the 4 per cent level. Additionally, the shape of the recovered isochrone ROI mode matches that of , which itself matches the mode shape of [15].

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 6, we report the growth rate of the fastest growing mode with various control parameters.

Previous studies have investigated the dependence of the ROIw.r.t.the normalised anisotropy radius, \(r_{\mathrm{a}}/b_{\mathrm{c}}\). A naive application of our libraries and the DFgiven in [63] suggests that the isochrone stability boundary is \({r_{\mathrm{a}}/b_{\mathrm{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_{\mathrm{a}}/b_{\mathrm{c}}\) decreases, regions of the \((J_r,L)\) plane become unpopulated [24]. Second, the gradient of the DFbecomes extremely large at the \({Q\!=\!0}\) boundary, such that numerical estimates readily break down. To circumvent this, [24] applied a smoothing to the radially-anisotropic form of the DF  hence considering the smoothed DF\({\tilde{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_{\mathrm{a}}/b_{\mathrm{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_{\mathrm{a}}/b_{\mathrm{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_{\mathrm{a}}/b_{\mathrm{c}}}\). [70] estimated an instability threshold at the half-mass radius of \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!1.67}\). [71] then identified an instability threshold of \({r_{\mathrm{a}}/b_{\mathrm{c}}\!=\!1.9}\), from simulations using a harmonic code [34]. [15] performed new simulations using the same harmonic code from [34], producing qualitative, but not quantitative, agreement in the trend of growth rates. The linear response calculations in [15], using a tailored basis, suggested that unstable modes may be found for \({r_{\mathrm{a}}/b_{\mathrm{c}}\!\le\!4}\). Future work may return to the question of stability with improved numerical treatment.

8.2 \({\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. [22] identified \({\ell\!=\!1}\) damped modes in [72] models using linear response theory. [22] also compared this prediction with perturbed \(N\)-body simulations, finding qualitative agreement. [25] measured the same \({\ell\!=\!1}\) mode in \(N\)-body simulations of a King model.

Extending the tools developed for the isochrone cluster, measured a damped mode for the isotropic isochrone model at \({\omega_{\mathrm{M}}/\Omega_{0}\!=\!0.0143\!-\!0.00142\mathrm{i}}\). As a first test, we reproduced this prediction using all the machinery from LinearResponse.jl. Table 7 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 , for which we recover \({ \omega_{\mathrm{M}}/ \Omega_0 \!=\! 0.0143 \!-\! 0.00143\mathrm{i}}\), in tight agreement with .

However, we noted that varying the control parameters of LinearResponse.jlresults in a drift of the mode’s frequency. The middle column of Table 7 lists the results for different runs. As the number of resonances considered increases (controlled by \(\texttt{n1max}\)), 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 [13], [73] 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 (19 ). Indeed, the divergence of \(P_{k} (\omega)\) (equation 19 ) 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.

9 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.jlto compute the frequencies from the potential (and its two derivatives), as well as the resonantly-aligned \({(u,v)}\) coordinates. Here, the versatility of LinearResponse.jlshines: the user only needs to define the DF.

The [29] potential is given by \[\psi(r) = - GM / \sqrt{b_{\mathrm{c}}^2+r^2} , \label{eq:Plummer95Pot}\tag{28}\] with \(b_{\mathrm{c}}\) the scale radius. The frequency scale for the Plummer sphere, \(\Omega_0\), is given by \[\Omega_0 = \Omega_1(r\to0) = 2\sqrt{ GM / b_{\mathrm{c}}}. \label{eq:plummeromega0}\tag{29}\] The isotropic DFis \[F(E) = \frac{24\sqrt{2}}{7\pi^3}\left(\frac{E}{E_0}\right)^{7 / 2},\] with the energy scale \(E_0=-GM/b_{\mathrm{c}}\).

To introduce radial anisotropy, we use the [74], [75] method to find a self-consistent DF. The dimensionless quantity \[Q = \frac{1}{E_0}\left(E+\frac{L^2}{2 r_{\mathrm{a}}^2}\right)\] is the single parameter in the DF, which reads¹³ \[F(Q) = \frac{24\sqrt{2}}{7\pi^3} \, Q^{7 / 2} \,\left[1-\frac{1}{r_{\mathrm{a}}^2}+\frac{7}{16r_{\mathrm{a}}^2Q^2}\right].\]

10 \(N\)-Body Technique↩︎

To compare our linear response theory predictions with \(N\)-body simulations, we generate realisations 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 8).¹⁴

3pt

The initial conditions are generated using PlummerPlus[35]. We use a lightly modified version of PlummerPlus.py, specifically to produce outputs that interface directly with EXP [38], a basis function expansion \(N\)-body code.¹⁵ The initial conditions are centred by the median positions and velocities in each dimension after sampling.

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 EXPresults in the same basis as the analytical [36] basis (cf.Appendix 6.2).

Conveniently, EXPuses 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, EXPemploys a block-step technique for subdividing timesteps in powers of 2. We allow a maximum of 7 divisions, corresponding to a minimum timestep of \({\sim\!1.5\!\times\!10^{-3}}\). The median timestep for particles is \({0.2/2^3\!=\!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 [38].

For computational expediency, we draw a relatively small number of particles per cluster,¹⁶ but realise six clusters per DFto 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, \({ \Delta \psi }\), we first compute the self-gravity in the non-axisymmetric harmonics. Given that the square of the EXPcoefficients is the self-gravity [38], and the magnitude of the amplitudes encodes \(\Delta \psi\), we may efficiently measure the potential fluctuations in either self-gravity or \(\Delta\psi\), simply using the output products of the simulation \[\Delta \psi^{\ell} = \sum_p a_p^{\ell} \, \psi_{p}^{\ell},\] 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} \big( a_p^{\ell} \big)^{2} . \label{eq:def95Al}\tag{30}\] To provide a natural scaling, harmonics with \({\ell\!>\!0}\) are normalised by the monopole to define \[\tilde{E}_{\ell} = E_{\ell} / E_{\ell = 0} . \label{eq:def95tE}\tag{31}\] In Fig. 3, we measure the slope of \({\ln\!\left(\tilde{E}_2\right)}\) 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_{\mathrm{M}}= \arctan\left( \mathrm{Im}[a_p^{\ell}] / \mathrm{Re}[a_p^{\ell}] \right) .\] 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 \({\Delta\psi^{\ell}/\psi_0}\) (where \(\psi_0\) is the equilibrium model potential) and \(\tilde{E}_{\ell}\) are of order 1%.

11 Razor-thin discs↩︎

11.1 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) \, \mathrm{e}^{\mathrm{i}\ell\phi},\, D_{p}^{\ell} (r) \, \mathrm{e}^{\mathrm{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 [41] (infinite radial support) and [3] (finite radial support), with a careful treatment of the normalisation. Naturally, it would be of interest to implement tailored basis functions [38], [64]. This is left for future work.

3pt

11.2 Zang discs↩︎

The stability analysis of tapered Mestel discs was pioneered by [4] and subsequently generalised by [43]. We refer to section 4.1 of [10] for a concise presentation of these discs, their DFsand potentials.

We place ourselves in the same unit system as in [10], hence setting \({V_0\!=\!G\!=\!R_{\mathrm{i}}\!=\!1}\). Within this unit system, a natural scale frequency is \({ \Omega_0\!=\!V_0/R_{\mathrm{i}}\!=\!1 }\). Here, \(\Omega_0\) does not stand for the (infinite) central radial frequency. As in [4], we consider “fully” self-gravitating discs (active fraction \({\xi\!=\!1}\)) with radial velocity dispersion such that \({q\!=\!6}\), while varying the inner taper exponent \(\nu_{\mathrm{t}}(=\!N)\). The outer taper is taken sufficiently far (\({R_{\mathrm{o}}\!=\!11.5}\)) and sharp (\({\mu_{\mathrm{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 softened Mestel potential \[\label{eq:truncated95Mestel95potential} \psi_\mathrm{M}(r) = \tfrac{1}{2}V_0^2 \ln [\varepsilon^2 + (r/R_0)^2],\tag{32}\] to prevent singularities for exactly radial orbits. In practice, \(\varepsilon\) is taken sufficiently small (\({\varepsilon\!=\!10^{-5}}\); \({R_0\!=\!1}\)) compared to the DF’s inner taper radius \({R_{\mathrm{i}}\!=\!1}\). We checked that this softening had no impact on the predicted linear response. Table 9 summarises the parameters used in LinearResponse.jl for the disc “Fiducial” numerical applications. In Table 10, 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.

References↩︎

---

1. Note that \(G^{\ell \mathbf{n}}_{pq} (u)\) in equation 3 slightly differs from \(\boldsymbol{\mathsf{G}}_{\mathbf{n}}(\mathbf{J}[u])\) in equation 1 . See Appendix 6.4.1 for details.↩︎

2. At present, the mappings are restricted to cored potentials. Extending them to cuspy potentials will be the purpose of future work.↩︎

3. The “Fiducial” control parameters for spherical calculations take the software defaults for OrbitalElements.jl(Table 3) and LinearResponse.jl(Table 5), along with a basis scale \({r_{\mathrm{b}}\!=\!5.0}\) and 100 radial basis functions. For \({\ell\!=\!2}\) calculations, this results in 62 resonances to calculate. For \({\ell\!=\!1}\) calculations, this results in 42 resonances to calculate. For the Plummer sphere, the fiducial value of \(\Omega_0\) is 2, cf. equation 29 .↩︎

4. 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 11.2). This softening is on a much smaller scale than the DF’s inner cut-out, and does not impact the linear response predictions.↩︎

5. The azimuthal harmonic number for discs is historically denoted \({m}\). Adapting here from the spherical case, we nonetheless denote it \(\ell\).↩︎

6. In practice, the library is currently limited to cored potentials. Its extension to cuspy potentials is left for future work.↩︎

7. In practice, the smoothness of the function dictates the error in the resulting calculations, rather than the integration scheme.↩︎

8. Generically, the radial functions, \(U_{p}^{\ell} (r)\), can always be taken to be real [26]. Conveniently, this makes \({ G(u) }\) a purely real function, as in equation 16 .↩︎

9. Equation 16 is the finite-interval version of the Hilbert transformation [65], hence the name of this software.↩︎

10. For smooth enough \(G(u)\), this is a robust approximation scheme [66].↩︎

11. We use the symmetry \({A_{pq}\!=\!A_{qp}}\) to reduce the computational cost.↩︎

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.↩︎

13. Equation (5) in [35] 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 10) includes this square.↩︎

14. Any unspecified parameters are taken as software defaults.↩︎

15. The code we use is available at https://github.com/michael-petersen/PlummerPlus, forked from the original.↩︎

16. No tests w.r.t.\(N\) were performed, though we expect that the measurement uncertainty will decrease like \(\sqrt{N}\).↩︎