© The Authors 2024

1. Introduction

Perturbative methods provide a powerful approach in stellar dynamics to elucidating observed features of galaxies and star clusters (Fridman & Polyachenko 1984; Binney & Tremaine 2008; Bertin 2014). Typically, the problem is transformed into solving an eigenvalue problem represented as a matrix equation, which remains a complex task (e.g. Kalnajs 1971, 1976, 1977; Polyachenko & Shukhman 1981; Polyachenko 2004, 2005). In such scenarios, it is highly beneficial to have a set of test analytical solutions. These can be used independently; for instance, to validate codes or address complex theoretical questions.

A well-known test perturbation involves shifting an entire spherical system. For instance, if a small shift, denoted by ζ, occurs along the z axis, where z = r cosθ and (r, θ, φ) are the coordinates of the spherical coordinate system, the resulting density and potential perturbations are $\rho (r,\theta )=-\zeta {\rho }_0^{\prime }(r)cos\theta$ and $\mathrm{\Phi }(r,\theta )=-\zeta {\mathrm{\Phi }}_0^{\prime }(r)\cos \theta$, respectively. This represents a dipole shift perturbation corresponding to the spherical harmonic P_l = 1(cosθ)=cosθ. It’s clear that the eigenfrequency, ω, for this perturbation is zero. This test has been frequently employed for code verification in stability studies (e.g., Tremaine 2005; Polyachenko & Shukhman 2015).

Another test solution, called “dilation mode”, recently proposed by Polyachenko & Shukhman (2023, hereafter PS23), applies to models with ergodic stable distribution functions, dF(E)/dE <  0, that contain a “single length parameter”. By fixing the total mass of the sphere, M, and altering the length parameter, we can achieve the same self-consistent equilibrium model but with a different length parameter. The simplicity of this solution allows for explicit expressions for the distribution function, potential, and density across all orders of perturbation theory.

This particularly helps to elucidate the concept of perturbation energy, a second-order quantity in amplitude, which is not calculable within linear theory. We’ve shown that the accurate expression for perturbation energy, constructed with second-order perturbations, doesn’t coincide with the well-known expression for perturbation energy, constructed as a bilinear form of the first-order perturbations, as is given in Sect. 5.4.2 of Binney & Tremaine (2008). These two forms of second-order energies – the pseudoenergy and the true energy – are both integrals of motion in the absence of external forces, but they can differ in both magnitude and sign.

In this study, we aim to derive the specified analytical solution using a well-established matrix method. However, we find that applying this method to stationary radial perturbations presents significant challenges.

At first glance, it may seem unnecessary to use matrix methods to investigate the spectra of eigenfrequencies in inherently stable models. After all, we know that in such cases the spectrum is continuous and consists solely of real frequencies of Van Kampen (1955) modes. However, the exploration of possible exponential decay (Landau damping) in these stable systems is linked to the superposition of van Kampen modes (Polyachenko et al. 2021; Lau & Binney 2021). This leads to the need to search for complex zeros with Im(ω)< 0 in the analytic continuation into the lower half-plane, ω, of the determinant, 𝒟(ω), which emerges in the matrix method. The introduction of our work, PS23, specifically discusses the rationale behind employing matrix methods in the analysis of stable models. It outlines the difficulties and challenges that arise when implementing this approach (for additional context and understanding, see also Weinberg 1994 and Polyachenko et al. 2021).

Section 2 provides a brief description of the dilation mode using the isochrone model (Hénon 1960) as an example. Section 3 delineates the problem and proposes a solution. Section 4 presents our conclusions. Finally, the appendices contain detailed derivations of some necessary relationships.

2. Models of spherical systems with a single length parameter and disturbances caused by its stretching or shrinking

2.1. The idea and an example

We consider an isotropic (ergodic) spherical model described by an equilibrium distribution function (DF) that contains a single characteristic length parameter, denoted as ℓ. For such models, the unperturbed potential and density take the form

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_0(r;\ell )=\frac{\mathit{MG}}{\ell }\varphi \left(\frac{r}{\ell }\right),\end{array}$$(1) $$\begin{array}{c}\hfill {\rho }_0(r;\ell )=\frac{M}{{\ell }^3}\varrho \left(\frac{r}{\ell }\right),\end{array}$$(2)

where ϕ(x) and 𝜚(x) are related by the Poisson equation,

$$\begin{array}{c}\hfill \frac{1}{x^2}\frac{\mathrm{d}}{\mathrm{d}x}\left[x^2\frac{\mathrm{d}\varphi (x)}{\mathrm{d}x}\right]=4\pi \varrho (x),\end{array}$$

and the distribution function is given by

$$\begin{array}{c}\hfill F_0(\mathcal{E};\ell )=\frac{1}{{(MG\ell )}^{3/2}}\mathcal{F}(\mathcal{E}),0\le \mathcal{E}\le \mathrm{\Psi }(0)\equiv -\varphi (0).\end{array}$$(3)

Here, minus dimensionless energy, ℰ, is also considered as a function of the velocity, v, radius, r, and length parameter, ℓ:

$$\begin{array}{cc}\hfill \mathcal{E}=\mathcal{E}(r,v;\ell )& =-\frac{\ell }{\mathit{MG}}E=-\frac{\ell }{\mathit{MG}}[\frac{1}{2}{v}^2+{\mathrm{\Phi }}_0(r,\ell )]\hfill \\ \hfill & =-\frac{\ell }{\mathit{MG}}[\frac{1}{2}{v}^2+\frac{\mathit{MG}}{\ell }\varphi \left(\frac{r}{\ell }\right)].\hfill \end{array}$$(4)

The dimensionless function, ℱ(ℰ), is normalized by the condition ∫ℱ d³r d³v = 1. It should be noted that not all spherical or even isotropic models fit the universal forms (Eqs. (1)–(3)). The King model (Michie & Bodenheimer 1963; King 1966), for instance, contains two length scales, which makes this model irrelevant.

It’s easy to understand that if we fix the total mass, M, but change ℓ, ℓ → ℓ + δℓ, we obtain the same self-consistent equilibrium model but with a different length parameter, ℓ. This means that the eigenfrequency, ω, of the mode corresponding to such a stretching or shrinking is zero. This fact can serve as a test for various codes when studying the dynamics of perturbations in spherical systems.

Concrete models with this form are known. For example, there is the isochrone model by Hénon (1960). For this model, the functions ϕ(x) and 𝜚(x) are

$$\begin{array}{c}\hfill \varphi (x)=-\frac{1}{1+a};\varrho (x)=\frac{1}{4\pi }\frac{1+2a}{{(1+a)}^2{a}^3},a=\sqrt{1+x^2},\end{array}$$(5)

where x = r/ℓ and the DF is given by

$$\begin{array}{cc}\hfill \mathcal{F}(\mathcal{E})& =\frac{1}{\sqrt{2}{(2\pi )}^3}\frac{\sqrt{\mathcal{E}}}{{[2(1-\mathcal{E})]}^4}\times [64{\mathcal{E}}^4-240{\mathcal{E}}^3+320{\mathcal{E}}^2-66\mathcal{E}+27\hfill \\ \hfill & +3(16{\mathcal{E}}^2+28\mathcal{E}-9)\frac{arcsin\sqrt{\mathcal{E}}}{\sqrt{\mathcal{E}(1-\mathcal{E})}}],\hfill \end{array}$$(6)

for $0\le \mathcal{E}\le \frac{1}{2}$. Several relevant models are discussed in the work by PS23, including those by Hernquist (1990) and Jaffe (1983). The polytrope model with the Plummer (1911) potential is also included, with explicit expressions provided for ϕ(x), 𝜚(x), and ℱ(ℰ). These models are also featured in the book by Binney & Tremaine (2008).

2.2. Perturbations induced by variations in the length parameter, ℓ

We consider the variations in the potential, density, and DF associated with the variation in the length parameter, ℓ. Assuming that ℓ = ℓ₀ + δℓ and expanding Eqs. (1)–(3) and (4) in Taylor series over ε, we obtain

$$\begin{array}{c}\hfill \mathrm{\Phi }(r;\ell )={\mathrm{\Phi }}_0(r;{\ell }_0)+\epsilon {\mathrm{\Phi }}_1(r;{\ell }_0)+{\epsilon }^2{\mathrm{\Phi }}_2(r;{\ell }_0)+\mathcal{O}({\epsilon }^3),\end{array}$$(7) $$\begin{array}{c}\hfill \rho (r;\ell )={\rho }_0(r;{\ell }_0)+\epsilon {\rho }_1(r;{\ell }_0)+{\epsilon }^2{\rho }_2(r;{\ell }_0)+\mathcal{O}({\epsilon }^3),\end{array}$$(8) $$\begin{array}{cc}\hfill F(\mathcal{E}(r,v;\ell ),r;\ell )& =F_0(\mathcal{E}(r,v;{\ell }_0);{\ell }_0)+\epsilon f_1(\mathcal{E}(r,v;{\ell }_0),r;{\ell }_0)\hfill \\ \hfill & +{\epsilon }^2{f}_2(\mathcal{E}(r,v;{\ell }_0),r;{\ell }_0)+\mathcal{O}({\epsilon }^3).\hfill \end{array}$$(9)

Here, ε = δℓ/ℓ₀, |ε|≪1 is the expansion parameter. Supposing further that G = M = ℓ₀ = 1, we have in the first order for the potential, density, and DF, respectively:

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_1(r)& =-[r\varphi (r)]\prime ,\hfill \end{array}$$(10) $$\begin{array}{cc}\hfill {\rho }_1(r)& =-[3\varrho (r)+r{\varrho }^{\prime }(r)],\hfill \end{array}$$(11) $$\begin{array}{cc}\hfill f_1(\mathcal{E},r)& =-\frac{3}{2}\mathcal{F}+{\mathcal{F}}^{\prime }(\mathcal{E})\{\mathcal{E}+{[r\varphi (r)]}^{\prime }\},\hfill \end{array}$$(12)

where now ℰ ≡ ℰ(r, v; ℓ₀) is minus unperturbed dimensionless energy. The prime symbol on functions denotes differentiation with respect to the corresponding argument. We note that since ℓ₀ = 1, we may not distinguish between the variables r and x ¹. Using the explicit form of perturbations, we can verify that the perturbation of mass is indeed zero:

$$\begin{array}{c}\hfill \int {\rho }_1(\mathbf{r}){\mathrm{d}}^3\mathbf{r}=\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}{f}_1(\mathbf{r},\mathbf{v})=0.\end{array}$$(13)

In what follows, we omit the index, “1”, of perturbations Φ₁, ρ₁, and f₁. Using Eq. (10), we can rewrite Eq. (12) as

$$\begin{array}{c}\hfill f(\mathcal{E},r)=-\frac{3}{2}\mathcal{F}(\mathcal{E})+{\mathcal{F}}^{\prime }(\mathcal{E})[\mathcal{E}-\mathrm{\Phi }(r)].\end{array}$$(14)

The perturbations of first-order Φ, ρ, and f represent a test perturbation, which corresponds to an eigenfrequency, ω = 0, in the eigenvalue problem.

If we study the dynamics of perturbations by solving the system of evolution equations for the Fourier amplitudes of the perturbed DF, f_n, specifying the initial DF in the form of Eq. (14) and the potential in the form of Eq. (10), we should obtain ∂f_n/∂t = 0. The test was conducted in PS23 using model (6). It confirmed the conservation of all harmonics, f_n, and consequently the entire DF: f(r, v, t)=f(r, v; 0).

Below, we demonstrate how the test is implemented in the traditional matrix method using the example of the isochrone model (6).

3. Using the dilation mode for testing in the matrix method

3.1. Clutton-Brock biorthonormal pairs for spherical systems

In order to apply the matrix method, we expanded the perturbed potential, Φ(r), and density, ρ(r), over the basis functions Φ^α(r) and ρ^α(r), where α = 1,  2,  3, …

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=\sum _{\alpha =1}^{\infty }{C}^{\alpha }{\mathrm{\Phi }}^{\alpha }(r),\rho (r)=\frac{1}{4\pi }\sum _{\alpha =1}^{\infty }{C}^{\alpha }{\rho }^{\alpha }(r),\end{array}$$(15)

where {Φ^α(r),ρ^α(r)} is a suitable set of potential-density basis pairs and C^α are expansion coefficients. These pairs should satisfy the equation

$$\begin{array}{c}\hfill (\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}+\frac{2}{r}\frac{\mathrm{d}}{\mathrm{d}r}){\mathrm{\Phi }}^{\alpha }(r)={\rho }^{\alpha }(r),\end{array}$$(16)

so that the Poisson equation is automatically satisfied, and the biorthogonality condition

$$\begin{array}{c}\hfill \underset{0}{\overset{\infty }{\int }}{\mathrm{\Phi }}^{\alpha }(r){\rho }^{\beta }(r)r^2\mathrm{d}r=-{\delta }^{\alpha \beta },\end{array}$$(17)

where δ^αβ is the Kronecker symbol². The coefficients, C^α, are found from the biorthogonality condition:

$$\begin{array}{c}\hfill C^{\alpha }=-\underset{0}{\overset{\infty }{\int }}\mathrm{\Phi }(r){\rho }^{\alpha }(r)r^2\mathrm{d}r=-4\pi \underset{0}{\overset{\infty }{\int }}{\mathrm{\Phi }}^{\alpha }(r)\rho (r)r^2\mathrm{d}r.\end{array}$$(18)

Since the isochrone model has an infinite radius and is regular at the center, it is convenient to choose the Clutton-Brock (1973) functions (hereafter CB) as biorthonormal potential-density basis pairs. In the general case when the potential and density depend on angular variables, ρ, Φ ∝ P_l(cosθ), the CB basis functions are expressed quite intricately in terms of Gegenbauer polynomials. However, in the case of radial perturbations (l = 0), they simplify to

$$\begin{array}{c}\hfill {\mathrm{\Phi }}^{\alpha }(r)=-\frac{2}{\sqrt{\pi (4{\alpha }^2-1)}}\frac{a(r)}{r}sin[2\alpha arcsin\frac{1}{a(r)}],\end{array}$$(19) $$\begin{array}{c}\hfill {\rho }^{\alpha }(r)=2\sqrt{\frac{4{\alpha }^2-1}{\pi }}\frac{1}{r{a}^3(r)}sin[2\alpha arcsin\frac{1}{a(r)}],\end{array}$$(20)

where $a(r)=\sqrt{1+r^2}$.

It is easy to verify that ρ^α and Φ^α are indeed related by the Eq. (16) and satisfy the biorthogonality condition (Eq. (17)). In Fig. 1, two basis functions of the potential Φ^α(r) are shown as an example, corresponding to α = 1 and α = 10.

Fig. 1. Basis functions of potential for α = 1 and α = 10 from the Clutton-Brock set. The number of nodes of Φα(r) is equal to α − 1.

In particular, for the isochrone model (6), considering that for it the perturbed potential, according to Eq. (10), is

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=\frac{1}{a(r)[1+a(r)]},\end{array}$$(21)

we have from Eq. (18)

$$\begin{array}{c}\hfill C^{\alpha }=-\underset{0}{\overset{\infty }{\int }}\mathrm{d}r{r}^2\frac{{\rho }^{\alpha }(r)}{a(1+a)},\end{array}$$(22)

or after substitution, r = cotψ, where $0\le \psi \le \frac{1}{2}\pi$,

$$\begin{array}{c}\hfill C^{\alpha }=-2\sqrt{\frac{4{\alpha }^2-1}{\pi }}\underset{0}{\overset{\pi /2}{\int }}\frac{cos\psi {sin}^2\psi sin(2\alpha \psi )}{1+sin\psi }\mathrm{d}\psi .\end{array}$$(23)

It should be noted that as r → ∞, the basis functions of the CB potential Φ^α(r) decrease as ∼1/r. Indeed, from Eq. (19) we have

$$\begin{array}{c}\hfill {\mathrm{\Phi }}^{\alpha }(r)\approx -\frac{q^{\alpha }}{r},q^{\alpha }=\frac{4\alpha }{\sqrt{\pi (4{\alpha }^2-1)}}.\end{array}$$(24)

As we will discuss further, this circumstance plays an important role in using the mode as a test. The conservation of the total mass of a spherical system during radial perturbation requires that the perturbed potential, Φ(r), and the density decrease faster than 1/r and 1/r³ at the periphery, respectively.

Therefore, when dealing with radial perturbations, it would be most optimal to use biorthogonal pairs of basis functions {Φ^α(r),ρ^α(r)} that each satisfy this condition. However, such basis biorthogonal pairs are not known in the literature, as is noted in the most comprehensive list of known basis pairs provided by Lilley et al. (2018). Hence, the condition for mass conservation should be expressed as nullifying the contribution, 𝒪(1/r), as r → ∞ in the total sum (Eq. (15)), describing the disturbance of potential. This can be derived from Eqs. (15) and (24)

$$\begin{array}{c}\hfill \sum _{\alpha =1}^{\infty }{q}^{\alpha }{C}^{\alpha }=0.\end{array}$$(25)

We emphasize that relation (25) expressing mass conservation is not an additional equation for the coefficients, C^α, but should be satisfied automatically if these coefficients are found by the matrix method correctly.

3.2. Matrix approach to radial disturbances

To describe quasiperiodic motions, it’s advantageous to use angles and actions instead of spherical coordinates and momentums. This is feasible for spherically symmetric potentials. The advantage of these variables is that angular variables, w, change uniformly over time with corresponding frequencies and are cyclic, while the unperturbed Hamiltonian, H₀ = −ℰ, depends only on the action variables J = (J_r, J_θ, J_φ) in the form H₀ = H₀(J_r, L)= − ℰ, where L = J_θ + |J_φ| is the angular momentum, and the radial action J_r is

$$\begin{array}{c}\hfill J_r=\frac{1}{2\pi }\oint v_r\mathrm{d}r.\end{array}$$(26)

The radial action and angular momentum are conserved along the unperturbed orbit. The radial frequency that corresponds to J_r is

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_r=\frac{\partial H_0(J_r,L)}{\partial J_r}.\end{array}$$(27)

It should be noted that in the isochrone model an analytic expression for Hamiltonian H₀(J_r, L) is available (see, e.g., Eq. (3.226a) in Binney & Tremaine 2008):

$$\begin{array}{c}\hfill H_0(J_r,L)=-\frac{1}{2}\frac{1}{[J_r+\frac{1}{2}(L+\sqrt{L^2+4}){]}^2},\end{array}$$(28)

so that the frequency, Ω_r, is solely dependent on energy and expressed as Ω_r(ℰ)=(2ℰ)^3/2.

The relationship between ℰ, J_r, and L allows us to use either the pair (ℰ, L) or (J_r, L) to specify the orbit. In the subsequent discussion, we use ϖ ≡ (ℰ, L) to designate orbits and the phase space. Once the orbit is specified, the radius of the particle becomes a function of the radial angle, w, which is conjugate to radial action J_r; that is, r = r(ϖ, w). From Eq. (14) we conclude that in case of the radial perturbations of interest the perturbed DF in action-angle variables is f = f(ϖ, w), while the unperturbed DF depends solely on energy, $\mathcal{F}=\mathcal{F}(\mathcal{E})$.

Now, it is helpful to recall what the matrix method for finding eigenvalues is. At this point, we will not consider perturbations of the general form, in other words, with potential and density depending also on angular variables θ and φ, as in the works by Polyachenko & Shukhman (1981), Weinberg (1991), Saha (1991), Bertin et al. (1994), but limit ourselves to the case of radial perturbations. The linearized collisionless Boltzmann equation in action-angle variables in this case has a simple form:

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}=-{\mathrm{\Omega }}_r\frac{\partial }{\partial w}(f+{\mathcal{F}}^{\prime }\mathrm{\Phi }).\end{array}$$(29)

We expanded the perturbed distribution function, f, potential, Φ, and density, ρ in a Fourier series in terms of the cyclic variable, w, and sought solutions in the form of e^−iωt:

$$\begin{array}{c}\hfill \omega f_n=n{\mathrm{\Omega }}_r(f_n+{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_n).\end{array}$$(30)

Here,

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_n(\mathit{\varpi })& =\oint \mathrm{d}w{\mathrm{e}}^{-\mathrm{i}nw}\mathrm{\Phi }(r(\mathit{\varpi },w)),\mathrm{\Phi }(r)=\frac{1}{2\pi }\sum _{n=-\infty }^{\infty }{\mathrm{\Phi }}_n{\mathrm{e}}^{\mathrm{i}nw},\hfill \end{array}$$(31) $$\begin{array}{cc}\hfill f_n(\mathit{\varpi })& =\oint \mathrm{d}w{\mathrm{e}}^{-\mathrm{i}nw}f(\mathcal{E},r(\mathit{\varpi },w)),f=\frac{1}{2\pi }\sum _{n=-\infty }^{\infty }{f}_n{\mathrm{e}}^{\mathrm{i}nw}.\hfill \end{array}$$(32)

Further, ϖ ≡ (ℰ,L) stands for the set of variables ℰ and L. The integration ∫d²ϖ (…) means integration over the accessible phase space. In the case of the isochrone model, this region is enclosed on the (ℰ,L)-plane between the lines $0\le \mathcal{E}\le \frac{1}{2}$, 0 <  L <  ∞, and the line of circular orbits $0\le L^2\le L_{\mathrm{circ}}^2(\mathcal{E})\equiv {(1-2\mathcal{E})}^2/(2\mathcal{E})$; that is,

$$\begin{array}{c}\hfill \int {\mathrm{d}}^2\mathit{\varpi }(\dots )={\int }_0^{1/2}\frac{\mathrm{d}\mathcal{E}}{{\mathrm{\Omega }}_r(\mathcal{E})}{\int }_0^{L_{\mathrm{circ}}^2(\mathcal{E})}\mathrm{d}(L^2)(\dots ).\end{array}$$(33)

Finally, for the isochrone potential, it is possible to analytically obtain a relation between the angular variable, w, and r (see also Ramond & Perez 2021). In the parametric form (parameter ξ) it reads:

$$\begin{array}{c}\hfill r(\mathit{\varpi },\xi )=\sqrt{{\left(\frac{1-pcos\xi }{2\mathcal{E}}\right)}^2-1},w(\mathit{\varpi },\xi )=\xi -psin\xi ,\end{array}$$(34)

where $p=\sqrt{{(1-2\mathcal{E})}^2-2\mathcal{E}{L}^2}$, − π ≤ ξ ≤ π, and − π ≤ w ≤ π. From Eq. (30), we have

$$\begin{array}{c}\hfill f_n=\frac{n{\mathrm{\Omega }}_r{\mathcal{F}}^{\prime }}{\omega -n{\mathrm{\Omega }}_r}{\mathrm{\Phi }}_n,\end{array}$$(35)

and for the density perturbation,

$$\begin{array}{c}\hfill \rho (r)=\frac{1}{2\pi }\int \mathrm{d}\mathbf{v}\sum _{n=-\infty }^{\infty }{f}_n{e}^{\mathrm{i}nw}=\frac{1}{2\pi }\int \mathrm{d}\mathbf{v}\sum _{n=-\infty }^{\infty }\frac{n{\mathrm{\Omega }}_r{\mathcal{F}}^{\prime }}{\omega -n{\mathrm{\Omega }}_r}{\mathrm{\Phi }}_n{e}^{\mathrm{i}nw}.\end{array}$$(36)

Multiplying both sides of Eq. (36) by Φ^α(r), integrating over the volume d³r = 4πr²dr, and taking into account that d³r d³v = d³J d³w = (2π)²d²ϖ dw, we find

$$\begin{array}{c}\hfill C^{\alpha }=-2\pi \int {\mathrm{d}}^2\mathit{\varpi }\sum _{n=-\infty }^{\infty }\frac{n{\mathrm{\Omega }}_r{\mathcal{F}}^{\prime }}{\omega -n{\mathrm{\Omega }}_r}{\mathrm{\Phi }}_n\oint \mathrm{d}w{\mathrm{\Phi }}^{\alpha }(r)e^{\mathrm{i}nw}.\end{array}$$(37)

Since r(−w) = r(w), we have

$$\begin{array}{c}\hfill \oint \mathrm{d}w{\mathrm{\Phi }}^{\alpha }(r)e^{\mathrm{i}nw}=\oint \mathrm{d}w{\mathrm{\Phi }}^{\alpha }(r)e^{-\mathrm{i}nw}={\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi }).\end{array}$$(38)

By substituting the expression for Φ_n with the expansion derived from Eq. (15), ${\mathrm{\Phi }}_n(\mathit{\varpi })={\displaystyle \sum _{\beta =1}^{\infty }}{C}^{\beta }{\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi })$, we get an infinite homogeneous system of linear equations³

$$\begin{array}{c}\hfill C^{\alpha }-\sum _{\beta =1}^{\infty }{\mathcal{M}}^{\alpha \beta }(\omega )C^{\beta }=0,\end{array}$$(39) $$\begin{array}{c}\hfill {\mathcal{M}}^{\alpha \beta }(\omega )=-2\pi \sum _{n=-\infty }^{\infty }\int {\mathrm{d}}^2\mathit{\varpi }{\mathrm{\Omega }}_r(\mathcal{E}){\mathcal{F}}^{\prime }(\mathcal{E})\frac{n{\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi }){\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi })}{\omega -n{\mathrm{\Omega }}_r(\mathcal{E})}.\end{array}$$(40)

From Eq. (39), we obtain the dispersion equation to find the eigenfrequency, ω:

$$\begin{array}{c}\hfill det\Vert {\delta }^{\alpha \beta }-{\mathcal{M}}^{\alpha \beta }(\omega )\Vert =0,\end{array}$$(41)

where δ^αβ is the Kronecker symbol. The matrix method for solving the eigenvalue problem consists of finding the zeros of the determinant in Eq. (41) and the eigenvector C^α.

It is important to emphasize that the equality of the sum in Eq. (25) to zero should be automatically satisfied if the coefficients, C^α, are found using the matrix method from the system of linear equations, provided the CB biorthogonal pairs Eqs. (19) and (20) are taken as basis functions.

3.3. Using the dilation mode to control the matrix method

At first glance, it may seem that a natural test for the correctness of the calculation within the matrix method using the stationary (ω = 0) dilation mode is to check if the determinant, ded‖D^αβ‖, is identically zero, where the matrix, D^αβ, is defined as (see Eq. (40))

$$\begin{array}{c}\hfill D^{\alpha \beta }={\delta }^{\alpha \beta }-M^{\alpha \beta },\end{array}$$(42)

with the matrix

$$\begin{array}{c}\hfill M^{\alpha \beta }={\mathcal{M}}^{\alpha \beta }(\omega =0)=2\pi \sum _{n\ne 0}\int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }(\mathcal{E}){\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi }){\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi }),\end{array}$$(43)

where the sum over n does not include the contribution from the zero harmonic, n = 0.

Of course, achieving an exact zero value for the determinant is impossible in principle, because in the numerical calculation, we are forced to truncate the infinite two-dimensional matrix, D^αβ, to finite dimensions N^α × N^α, which means keeping a finite number of N^α basis functions in the perturbed potential:

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=\sum _{\alpha =1}^{N^{\alpha }}{C}^{\alpha }{\mathrm{\Phi }}^{\alpha }(r).\end{array}$$(44)

However, the expectation that the determinant of the truncated matrix, D^αβ (Eq. (42)), would tend to zero, as N^α → ∞, was not confirmed. It turned out that this is not the case. This determinant tends not to zero, but to a finite limit, det ≈ −0.018, which is achieved already at N^α ≈ 40, and with a further increase in N^α to huge values, up to N^α = 5000, it remains approximately equal to the value mentioned above. It is also worth noting that the absence of convergence of the determinant to zero cannot be justified by the limited number of terms in the sum over the harmonics, n, in the matrix (43), which we usually have to resort to in the numerical procedure. The point is that for this mode, there is an elegant analytical way to calculate the matrix, M^αβ, by first taking into account all the harmonics, n, including the zero one, n = 0, (we call such a semi-finished matrix ${\stackrel{\sim }{M}}^{\alpha \beta }$), and then subtracting the contribution of the zero harmonic numerically:

$$\begin{array}{c}\hfill M^{\alpha \beta }={\tilde{M}}^{\alpha \beta }-M_0^{\alpha \beta },\end{array}$$(45)

where

$$\begin{array}{cc}\hfill {\tilde{M}}^{\alpha \beta }& =2\pi \sum _{n=-\infty }^{\infty }\int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }(\mathcal{E}){\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi }){\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi }),\hfill \end{array}$$(46) $$\begin{array}{cc}\hfill M_0^{\alpha \beta }& =2\pi \int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_0^{\alpha }(\mathit{\varpi }){\mathrm{\Phi }}_0^{\beta }(\mathit{\varpi }).\hfill \end{array}$$(47)

This means that in this way, we obtain the opportunity to include all harmonics in the matrix M^αβ, except for the zero one, and avoid the need to limit the number of included harmonics, n (which also, in principle, could be a reason for the determinant not converging to zero). This is what was actually done. The method of calculation and the explicit expression for the matrix ${\stackrel{\sim }{M}}^{\alpha \beta }$ are given in Appendix A.

Thus, it creates the impression that the initial idea of using the dilation mode as a test by checking the identity of the determinant of the matrix D^αβ (42) to zero seems inconsistent.

To understand what is going on, we need to go back to the original Eq. (30) for the harmonics. At ω = 0, it implies

$$\begin{array}{c}\hfill f_n=-{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_n,\end{array}$$(48)

for all n ≠ 0, which is the same equation that we obtain from Eq. (35) by setting ω = 0. However, at n = 0, the situation is completely different.

Indeed, from Eq. (30), it can be seen that at ω ≠ 0, the zero harmonic of the distribution function becomes zero, f_{n = 0}(ϖ) = 0, and is thus uniquely determined. On the other hand, at ω = 0, the Eq. (30) for n = 0 turns into the identity 0 = 0 for any f_{n = 0}. This means that at ω = 0, we cannot in principle obtain the zero harmonic of the distribution function. If we formally set n = 0 in Eq. (48), we obtain an incorrect expression for f_{n = 0}, namely $f_{n=0}=-{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_{n=0}$, while the correct expression for f_{n = 0}, as can be seen from Eq. (14) in the case of our stationary dilation mode, is

$$\begin{array}{c}\hfill f_{n=0}=2\pi [-\frac{3}{2}\mathcal{F}(\mathcal{E})+{\mathcal{F}}^{\prime }(\mathcal{E})]-{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_{n=0}(\mathit{\varpi }),\end{array}$$(49)

which differs from Eq. (48) by an additional term. In other words, it is not equal to zero, as for modes with ω = 0, but it is also not equal to the value $-{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_{n=0}(\mathit{\varpi })$ that formally follows from Eq. (48). In the case of the isochrone model, the calculation of the zero harmonic of the potential Φ_{n = 0}(ϖ) results in the final explicit form of f_{n = 0}:

$$\begin{array}{c}\hfill f_{n=0}(\mathit{\varpi })=2\pi \{[-\frac{3}{2}\mathcal{F}(\mathcal{E})+\mathcal{E}{\mathcal{F}}^{\prime }(\mathcal{E})]-{\mathcal{F}}^{\prime }(\mathcal{E})\frac{{(2\mathcal{E})}^{3/2}}{\sqrt{L^2+4}}\}.\end{array}$$(50)

The non-vanishing of the zero harmonic in the case of the radial stationary mode, as we will show below, leads to the fact that the system of linear equations for the coefficients, C^α, becomes inhomogeneous. This means that the core matrix Eq. (41) is incorrect for ω = 0. Therefore, we do not have the right to demand the determinant, ded‖D^αβ‖, to be identically zero as a test.

In simpler terms, we can formulate the reason for the incorrectness of the procedure as follows. We cannot first set ω = 0 in the right-hand side of Eq. (35), and then cancel the numerator and denominator by n, because for n = 0 this leads to an error in computing the zero harmonic. Although for all other harmonics f_{n ≠ 0}, the procedure is correct.

Nevertheless, it turns out that a slight modification of the system of linear equations, which leads to the dispersion Eq. (42) for ω = 0, makes it possible to use this mode as a test. We will demonstrate this.

To do this, first, we will obtain a correct system of equations for the coefficients, C^α, which is valid in the case of stationary perturbations. We will see that it is inhomogeneous, and this is precisely the reason why the determinant, ded‖D^αβ‖, differs from zero.

Going back to the expression (14) for the perturbed DF, $f=[-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }(\mathcal{E})]-{\mathcal{F}}^{\prime }(\mathcal{E})\mathrm{\Phi }(r)$, we can use it to derive a system of equations for the coefficients, C^α, of the potential Φ(r) expansion in terms of basis functions. Recalling that, according to Eq. (18),

$$\begin{array}{c}\hfill C^{\alpha }=-\int {\mathrm{d}}^3\mathbf{r}\rho (r){\mathrm{\Phi }}^{\alpha }(r)=-\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}f{\mathrm{\Phi }}^{\alpha }(r),\end{array}$$(51)

and that d³r d³v = d³J d³w = (2π)²d²ϖ dw, we get

$$\begin{array}{c}\hfill C^{\alpha }=-\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}f{\mathrm{\Phi }}^{\alpha }(r)=-{(2\pi )}^2\int {\mathrm{d}}^2\varpi \oint \mathrm{d}wf{\mathrm{\Phi }}^{\alpha }.\end{array}$$(52)

Expanding f in Eq. (52) into harmonics (see Eqs. (31) and (32)), we obtain

$$\begin{array}{c}\hfill C^{\alpha }=-2\pi \sum _n\int {\mathrm{d}}^2\mathit{\varpi }{f}_n(\mathit{\varpi }){\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi }).\end{array}$$(53)

Using Eq. (14) for the harmonics, f_n(ϖ), we have

$$\begin{array}{c}\hfill f_n(\mathit{\varpi })=-{\mathcal{F}}^{\prime }{\mathrm{\Phi }}_n(\mathit{\varpi })+2\pi [-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }(\mathcal{E})]{\delta }_{n,0},\end{array}$$(54)

or, taking into account that ${\mathrm{\Phi }}_n(\mathit{\varpi })={\displaystyle \sum C^{\beta }}{\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi })$,

$$\begin{array}{c}\hfill f_n(\mathit{\varpi })=-{\mathcal{F}}^{\prime }\sum _{\beta }{C}^{\beta }{\mathrm{\Phi }}_n^{\beta }+2\pi [-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }(\mathcal{E})]{\delta }_{n,0},\end{array}$$(55)

where δ_nm is the Kronecker delta. Substituting Eq. (55) into the right-hand side of Eq. (53), we obtain a system of inhomogeneous equations,

$$\begin{array}{c}\hfill C^{\alpha }-\sum _{\beta =1}^{\infty }{\tilde{M}}^{\alpha \beta }{C}^{\beta }={\tilde{A}}^{\alpha },\end{array}$$(56)

where ${\stackrel{\sim }{M}}^{\alpha \beta }$ is defined by Eq. (46), and

$$\begin{array}{c}\hfill {\tilde{A}}^{\alpha }=-{(2\pi )}^2\int {\mathrm{d}}^2\mathit{\varpi }{\mathrm{\Phi }}_{n=0}^{\alpha }(\mathit{\varpi })[-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }(\mathcal{E})].\end{array}$$(57)

Appendix A demonstrates how to obtain simple expressions for ${\stackrel{\sim }{M}}^{\alpha \beta }$ and ${\stackrel{\sim }{A}}^{\alpha }$ without resorting to action-angle variables.

In the double sum over β and n, denoted as ${\displaystyle \sum _{\beta }}{\stackrel{\sim }{M}}^{\alpha \beta }{C}^{\beta }$ in Eq. (56), we isolate the term where n = 0. Since ${\displaystyle \sum _{\beta }}{\mathrm{\Phi }}_{n=0}^{\beta }(\mathit{\varpi })C^{\beta }={\mathrm{\Phi }}_{n=0}(\mathit{\varpi })$, according to Eq. (47), the contribution of this term to the left-hand side of Eq. (56) is

$$\begin{array}{c}\hfill -\sum _{\beta }{M}_0^{\alpha \beta }{C}^{\beta }=-(2\pi )\int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }(\mathcal{E}){\mathrm{\Phi }}_{n=0}^{\alpha }(\mathit{\varpi }).\end{array}$$(58)

By shifting it to the right-hand side of Eq. (56), we obtain

$$\begin{array}{c}\hfill C^{\alpha }-\sum _{\beta =1}^{\infty }{M}^{\alpha \beta }{C}^{\beta }=A^{\alpha },\end{array}$$(59)

where

$$\begin{array}{c}\hfill A^{\alpha }={\tilde{A}}^{\alpha }+(2\pi )\int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }(\mathcal{E}){\mathrm{\Phi }}_{n=0}^{\alpha }(\mathit{\varpi }).\end{array}$$(60)

By utilizing the expression (Eq. (49)) for the zero harmonic f_{n = 0} and the expression (57) for ${\stackrel{\sim }{A}}^{\alpha }$, we can formulate the final system of inhomogeneous linear equations (Eq. (59)) for the coefficients, C^α, as follows:

$$\begin{array}{cc}& \sum _{\beta }{D}^{\alpha \beta }{C}^{\beta }=A^{\alpha },\hfill \end{array}$$(61) $$\begin{array}{cc}& D^{\alpha \beta }={\delta }^{\alpha \beta }-M^{\alpha \beta }\equiv {\delta }^{\alpha \beta }-2\pi \int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }\sum _{n\ne 0}{\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi }){\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi }),\hfill \end{array}$$(62) $$\begin{array}{cc}& A^{\alpha }=-(2\pi )\int {\mathrm{d}}^2\mathit{\varpi }{\mathrm{\Phi }}_{n=0}^{\alpha }(\mathit{\varpi })f_{n=0}(\mathit{\varpi }).\hfill \end{array}$$(63)

Once again, it’s important to emphasize that the inhomogeneity of the system of equations (Eq. (61)) (meaning the nonzero value of A^α) is entirely due to the zero harmonic of the perturbed distribution function, f_{n = 0}(ϖ), which, as we saw earlier, equals zero when ω = 0. This is why the inhomogeneous system of equations (Eq. (61)), which we have in the case of ω = 0, cannot be obtained by a simple limiting transition from the homogeneous system (Eq. (39)). In other words, the case of stationary radial perturbations is a special case.

The left-hand side of the inhomogeneous system (Eq. (61)) is described by the same matrix, D^αβ, the zero determinant of which we assumed in order to use it as a test. If this determinant were zero, it would imply that this system has no non-trivial solutions. However, since we know that the coefficients, C^α, are not zero (they are determined independently by the direct expansion of the perturbed potential, Φ(r) (see Eq. (23)), it becomes clear that the determinant, ded‖D^αβ‖ ≠ 0. As was noted earlier, this fact was already confirmed by direct numerical calculation.

Thus, we manage to understand the reason for the paradox we encountered earlier: on the one hand, we have a matrix Eq. (41) that is supposed to hold, and on the other hand, assuming ω = 0 in it, we don’t obtain the expected identity.

In Appendix B, it is shown that if we multiply each of the equations in the system (Eq. (61)) by the weight factor, $q^{\alpha }=4\alpha /\sqrt{\pi (4{\alpha }^2-1)}$, and then sum over α, we obtain the equation

$$\begin{array}{c}\hfill \sum _{\beta =1}^{\infty }{q}^{\beta }{C}^{\beta }=0.\end{array}$$(64)

This aligns perfectly with the mass conservation Eq. (25) that we previously discussed. This fact means that the solution of the system (Eq. (61)) automatically satisfies the mass conservation condition (Eq. (25)), as is required by a correct system. On the other hand, it provides a hint as to how to modify the inhomogeneous system of Eq. (61) in order to make the determinant of the modified system a test of the correctness of the matrix equation.

Indeed, given that the solution to the inhomogeneous system (Eq. (61)) inherently satisfies the condition (Eq. (25)), it’s clear that if we augment the system (Eq. (61)) with the Eq. (25), it will transform into a system of linearly dependent equations. Consequently, its determinant must be zero. The vanishing of the determinant of this augmented system will serve as the test.

The matter is complicated by the fact that, for a numerical check of the determinant’s equality to zero, we are forced to deal with a truncated system, limiting ourselves to a finite number of coefficients, C^α. The truncation process in this case requires caution. Let the number of coefficients, C^α, retained in the system be denoted as N^α. Denoting the finite sums as

$$\begin{array}{c}\hfill \sum _{\alpha =1}^{N^{\alpha }}{q}^{\alpha }{M}^{\alpha \beta }={\mathcal{M}}_{N^{\alpha }}^{\beta },\sum _{\alpha =1}^{N^{\alpha }}{q}^{\alpha }{A}^{\alpha }={\mathcal{A}}_{N^{\alpha }},\end{array}$$(65)

we find (see Eq. (B.2)) that the limit of these sums as N^α → ∞ is zero:

$$\begin{array}{c}\hfill \underset{N^{\alpha }\to \infty }{lim}{\mathcal{M}}_{N^{\alpha }}^{\beta }={\mathfrak{M}}^{\beta }=0,\underset{N^{\alpha }\to \infty }{lim}{\mathcal{A}}_{N^{\alpha }}=\mathfrak{A}=0.\end{array}$$(66)

We denote by $\mathfrak{D}(N^{\alpha })$ the determinant of size N^α × N^α of the truncated modified system of equations for finding N^α coefficients, C^α. The modification to the system of N^α equations is that now it contains N^α − 1 equations:

$$\begin{array}{c}\hfill \sum _{\beta =1}^{N^{\alpha }}{D}^{\alpha \beta }{C}^{\beta }\equiv \sum _{\beta =1}^{N^{\alpha }}[{\delta }^{\alpha \beta }-2\pi \int {\mathrm{d}}^2\mathit{\varpi }{\mathcal{F}}^{\prime }\sum _{n\ne 0}{\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi }){\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi })]C^{\beta }=A^{\alpha },\end{array}$$(67)

where α ranges from 0 to N^α − 1, and the last N^α-th equation is replaced by the additional equation⁴

$$\begin{array}{c}\hfill \sum _{\beta =1}^{N^{\alpha }-1}{q}^{\beta }{C}^{\beta }=0.\end{array}$$(68)

Therefore, a correct test for the matrix equation in the case of stationary perturbations, ω = 0, should consist of confirming the tendency of $\mathfrak{D}(N^{\alpha })$ to approach zero when N^α tends to infinity:

$$\begin{array}{c}\hfill \mathfrak{D}(N^{\alpha })=\left|\begin{array}{ccccc}{D}^{11}& D^{12}& \dots & D^{1,N^{\alpha }-1}& D^{1,N^{\alpha }}\\ & & & & \\ D^{21}& D^{22}& \dots & D^{2,N^{\alpha }-1}& D^{2,N^{\alpha }}\\ \dots & \dots & \dots & \dots & \dots \\ D^{N^{\alpha }-1,1}& D^{N^{\alpha }-1,2}& \dots & D^{N^{\alpha }-1,N^{\alpha }-1}& D^{N^{\alpha }-1,N^{\alpha }}\\ & & & & \\ q^1& q^2& \dots & q^{N^{\alpha }-1}& 0\end{array}\right|.\end{array}$$(69)

We performed this test for the isochrone model and confirmed that $\mathfrak{D}(N^{\alpha })$ indeed tends to zero with the growth of N^α as 1/(N^α)². Figure 2 shows the numerically found dependence of the determinant, $\mathfrak{D}$ (Eq. (69)), on N^α.

Fig. 2. Dependence on Nα of the determinant, 𝔇 (Eq. (69)), of the matrix of size Nα × Nα of the linear system of equations consisting of Nα − 1 equations of the inhomogeneous system (Eq. (67)) plus Eq. (68) – blue line. It is well approximated by fit −0.4/(Nα)2 – dotted orange line.

Another correctness check is the calculation of the coefficients, C^α, from the inhomogeneous system (Eq. (67)), where we set the coefficient, C^Nα = 0. Thus, we reduced the system of N^α linearly dependent Eqs. (67) + (68) to a system of N^α − 1 linearly independent equations and lowered the rank of the matrix. The coefficients, C^α, for α = 1,2,…,N^α − 1 (where N^α = 1000) have been determined and are displayed in Fig. 3. Upon comparison, it is evident that they are indistinguishable from those obtained by directly expanding the known perturbed potential, Φ(r), using Eq. (23), as expected.

Fig. 3. Expansion coefficients, Cα, of the perturbed potential Eq. (21) over Clutton-Brock (1973) basis functions (Eq. (19)) found from the inhomogeneous system (Eq. (67)). The first coefficient, Cα = 1 = −0.30339, is negative and the rest are positive. The coefficients are well described by the fit 0.54/α2.

4. Conclusion

The dilation mode, first introduced in PS23, provides explicit formulas for the distribution function, potential, and density across all orders of the perturbation theory. This allows us to explore various aspects of the perturbation theory as it applies to spherical stellar systems. This mode is present in ergodic systems with a single length parameter. When using matrix methods to solve the eigenvalue problem, this perturbation should yield an eigenfrequency of ω = 0 and known eigenfunctions.

We used the dilation mode in the isochrone model (Hénon 1960) as a test case to try and replicate it using the matrix method proposed by Polyachenko & Shukhman (1981). The method requires the determinant, det∥D^αβ(ω)∥, to be zero for the eigenfrequencies. The uniqueness of the radial case when using the matrix method was previously highlighted by Bertin et al. (1994) and Polyachenko & Shukhman (2015). We show that the matrix method indeed works for radial disturbances, yielding both the disturbance frequency and the potential’s eigenfunction, but after accounting for the specific features of the l = 0 case.

When we directly substituted ω = 0 into the expression for the matrix elements, D^αβ(ω), we obtained a nonzero determinant. The nonzero determinant was due to an unconventional relationship (Eq. (49)) between the perturbed harmonics, n = 0, of the distribution function’s expansion and the potential. Specifically, it was observed that f_n = 0 ≠ 0. This observation led to inhomogeneity in the equations for the potential decomposition coefficients, C^α.

The peculiarity of this situation can be explained as follows. For radial perturbations where ω ≠ 0 (i.e., ∂/∂t ≠ 0), the zero harmonic of the perturbed DF, f_n = 0, is zero. This is because the contribution of the zero harmonic of the perturbed potential does not generate a disturbing force since there is no gradient, ∂/∂w = ∂/∂r = 0. However, the stationary case where ω = 0 is different. In this case, the zero harmonic of the DF does not vanish: since the perturbed state is the same self-consistent model but with adjusted energy and scale length, its new DF is

$$F(\mathcal{E}(\ell ={\ell }_0+\delta \ell );\ell ={\ell }_0+\delta \ell )=F_0(\mathcal{E}({\ell }_0);{\ell }_0)+\epsilon f(\mathcal{E}({\ell }_0),r).$$

This results in the emergence of the zero harmonic f_n = 0 = ∮f(ℰ, r) dw.

At first glance, it might seem that the matrix method does not work. However, an additional constraint on the coefficients, C^α, is imposed by the need for total mass conservation, as is expressed by the relationship (Eq. (25)). Naturally, this should be automatically fulfilled if the system of equations for C^α is correct. In fact, Appendix B confirms that this is indeed the case. When this condition is added to the matrix of the system of equations, it renders the system linearly dependent. As a result, the determinant of the matrix becomes zero. This condition is tied to the fact that for l = 0, the potential decomposition functions from the biorthonormal set decreases at the periphery as r⁻¹. For the coefficients, C^α, that do not satisfy this relationship, this could imply that the disturbed potential of the system has an asymptotic behavior ∝r⁻¹, which would correspond to a change in the system’s mass.

For the case of non-radial disturbances l >  0, there are two key points to note. Firstly, there is no inhomogeneity in the equations for C^α. Secondly, although the conservation of mass also takes place, it is not expressed in the relationship of coefficients, C^α. This is because the potential decomposition functions from the biorthonormal set decreases at the periphery as r^{−1 − l}; in other words, faster than Coulomb’s. Therefore, they are not obligated to follow a relationship similar to Eq. (25). In other words, there is no correlation between the coefficients, C^α, unlike in the case of radial perturbations where l = 0.

Acknowledgments

The work was carried out with the financial support of the Program of the Presidium of the Russian Academy of Sciences No. 28 “Space: Studies of Fundamental Processes and their Interrelations” (subprogram II “Astrophysical Objects as Cosmic Laboratories”), as well as the Ministry of Science and Higher Education of the Russian Federation (Ilia Shukhman).

References

Appendix A: Derivation of matrix ${\stackrel{\sim }{M}}^{\alpha \beta }$

It is possible to obtain the matrix ${\stackrel{\sim }{M}}^{\alpha \beta }$, which includes all the Fourier harmonics, e^inw, quite simply, without passing to the action-angle variables and even without the procedure of expansion in Fourier harmonics itself. For the perturbed DF we have

$$\begin{array}{c}\hfill f=-{\mathcal{F}}^{\prime }\mathrm{\Phi }(r)+(-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }).\end{array}$$(A.1)

We multiply both sides of the equation by −Φ^α(r) and integrated over the entire phase volume. Taking into account that (see (52)) we obtain

$$\begin{array}{cc}\hfill C^{\alpha }& =\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}{\mathcal{F}}^{\prime }\mathrm{\Phi }(r){\mathrm{\Phi }}^{\alpha }(r)\hfill \\ \hfill & -\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}(-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }){\mathrm{\Phi }}^{\alpha }(r).\hfill \end{array}$$(A.2)

Substitution of Φ(r) in the form of series

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=\sum _{\beta }{C}^{\beta }{\mathrm{\Phi }}^{\beta }(r)\end{array}$$(A.3)

into the first term on the right-hand side of (A.2) yields

$$\begin{array}{cc}& C^{\alpha }=\sum _{\beta }\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}{\mathcal{F}}^{\prime }{\mathrm{\Phi }}^{\alpha }(r){\mathrm{\Phi }}^{\beta }(r)C^{\beta }\hfill \\ \hfill & -\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}(-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }){\mathrm{\Phi }}^{\alpha }(r),\hfill \end{array}$$(A.4)

or, in a compact form,

$$\begin{array}{c}\hfill C^{\alpha }=\sum _{\beta }{\tilde{M}}^{\alpha \beta }{C}^{\beta }+{\tilde{A}}^{\alpha },\end{array}$$(A.5)

where

$$\begin{array}{c}\hfill {\tilde{M}}^{\alpha \beta }=\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}{\mathcal{F}}^{\prime }{\mathrm{\Phi }}^{\alpha }(r){\mathrm{\Phi }}^{\beta }(r),\end{array}$$(A.6) $$\begin{array}{c}\hfill {\tilde{A}}^{\alpha }=-\int {\mathrm{d}}^3\mathbf{r}{\mathrm{d}}^3\mathbf{v}(-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime }){\mathrm{\Phi }}^{\alpha }(r).\end{array}$$(A.7)

Since the integrands (A.6) and (A.7) depend only on ℰ and r, we proceed to integration over these variables:

$$\begin{array}{c}\hfill {\tilde{M}}^{\alpha \beta }=4\pi \int \mathrm{d}r{r}^2{\mathrm{\Phi }}^{\alpha }(r){\mathrm{\Phi }}^{\beta }(r)I_1(r),\end{array}$$(A.8) $$\begin{array}{c}\hfill {\tilde{A}}^{\alpha }=-4\pi \int \mathrm{d}r{r}^2{\mathrm{\Phi }}^{\alpha }(r)I_2(r),\end{array}$$(A.9)

where

$$\begin{array}{c}\hfill I_1(r)=4\pi \underset{0}{\overset{\mathrm{\Psi }}{\int }}\mathrm{d}\mathcal{E}{\mathcal{F}}^{\prime }(\mathcal{E})\sqrt{2\mathrm{\Psi }-2\mathcal{E}},\end{array}$$(A.10) $$\begin{array}{c}\hfill I_2(r)=\underset{0}{\overset{\mathrm{\Psi }}{\int }}\mathrm{d}\mathcal{E}(-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime })\sqrt{2\mathrm{\Psi }-2\mathcal{E}}.\end{array}$$(A.11)

The integrals I_1, 2(r) are taken analytically. Below, we show how to take them. We give the result,

$$\begin{array}{c}\hfill I_1(r)=\frac{1}{4\pi }\frac{8a^2+9a+3}{a^4(1+a)},\end{array}$$(A.12) $$\begin{array}{c}\hfill I_2(r)=\frac{1}{4\pi }\frac{2a^2+6a+3}{a^4{(1+a)}^2},\end{array}$$(A.13)

where $a=\sqrt{1+r^2}$. Next, we go in (A.6) and (A.7) from variable r to variable ψ:

$$\begin{array}{c}\hfill r=cot\psi ,a=\frac{1}{sin\psi },{\mathrm{\Phi }}^{\alpha }(r(\psi ))=-\frac{2}{\sqrt{\pi (4{\alpha }^2-1)}}\frac{sin(2\alpha \psi )}{cos\psi },\end{array}$$(A.14)

and find the final expressions for ${\stackrel{\sim }{M}}^{\alpha \beta }$ and ${\stackrel{\sim }{A}}^{\alpha }$ in the form of simple single integrals

$$\begin{array}{cc}\hfill {\overline{M}}^{\alpha \beta }& =\frac{4}{\pi \sqrt{(4{\alpha }^2-1)(4{\beta }^2-1)}}\underset{0}{\overset{\pi /2}{\int }}\mathrm{d}\psi \frac{sin(2\alpha \psi )sin(2\beta \psi )}{sin\psi (1+sin\psi )}\hfill \\ \hfill & \times (8+9sin\psi +3{sin}^2\psi ),\hfill \end{array}$$(A.15) $$\begin{array}{cc}\hfill {\overline{A}}^{\alpha }& =\frac{2}{\sqrt{\pi (4{\alpha }^2-1)}}\underset{0}{\overset{\pi /2}{\int }}\mathrm{d}\psi \frac{cos\psi sin(2\alpha \psi )}{{(1+sin\psi )}^2}\hfill \\ \hfill & \times (2+8sin\psi +3{sin}^2\psi ).\hfill \end{array}$$(A.16)

Finally, we show how the integrals, I_1, 2(r), are calculated.

(i) Calculation of I₁(r). Performing integration by parts, we find that

$$\begin{array}{cc}\hfill I_1(r)& =4\pi \underset{0}{\overset{\mathrm{\Psi }}{\int }}\mathrm{d}\mathcal{E}{\mathcal{F}}^{\prime }\sqrt{2\mathrm{\Psi }-2\mathcal{E}}=4\pi \underset{0}{\overset{\mathrm{\Psi }}{\int }}\frac{\mathcal{F}\mathrm{d}\mathcal{E}}{\sqrt{2\mathrm{\Psi }-2\mathcal{E}}}\hfill \\ \hfill & =4\pi \frac{d}{\mathrm{d}\mathrm{\Psi }}\underset{0}{\overset{\mathrm{\Psi }}{\int }}\mathrm{d}\mathcal{E}\mathcal{F}\sqrt{2\mathrm{\Psi }-2\mathcal{E}}=\frac{\mathrm{d}\varrho }{\mathrm{d}\mathrm{\Psi }}=\frac{{\varrho }^{\prime }}{{\mathrm{\Psi }}^{\prime }}=-\frac{{\varrho }^{\prime }}{{\varphi }^{\prime }}.\hfill \end{array}$$(A.17)

Since for the isochrone model (6) the explicit expressions for the unperturbed density, 𝜚, and potential, ϕ = −Ψ, are (see (5))

$$\begin{array}{c}\hfill \varrho =\frac{1}{4\pi }\frac{1+2a}{{(1+a)}^2{a}^3},\varphi =-\mathrm{\Psi }=-\frac{1}{1+a},\end{array}$$(A.18)

we obtain for their derivatives

$$\begin{array}{c}\hfill {\varrho }^{\prime }=-\frac{1}{4\pi }\frac{8a^2+9a+3}{a^4{(1+a)}^3}{a}^{\prime },{\varphi }^{\prime }=\frac{a^{\prime }}{{(1+a)}^2},\end{array}$$(A.19)

and finally find I₁ in the form (A.12).

(ii) Calculation of I₂(r). Integrating by parts, we get

$$\begin{array}{c}\hfill I_2(r)=\underset{0}{\overset{\mathrm{\Psi }}{\int }}\mathrm{d}\mathcal{E}(-\frac{3}{2}\mathcal{F}+\mathcal{E}{\mathcal{F}}^{\prime })\sqrt{2\mathrm{\Psi }-2\mathcal{E}}=-3\varrho +\varphi \frac{{\varrho }^{\prime }}{{\varphi }^{\prime }}.\end{array}$$(A.20)

Again, using (A.18) and (A.19) yields I₂ in the form (A.13).

Appendix B: Derivation of the conservation of mass equation

We show in more detail how the condition for conservation of mass (25) automatically follows from an inhomogeneous system (61). We multiply in (61) the equation with number α by

$$q^{\alpha }=\frac{4\alpha }{\sqrt{\pi (4{\alpha }^2-1)}}$$

and sum over all α. We get

$$\begin{array}{c}\hfill \sum _{\beta =1}^{\infty }(q^{\beta }-{\mathfrak{M}}^{\beta })C^{\beta }=\mathfrak{A},\end{array}$$(B.1)

where

$${\mathfrak{M}}^{\beta }\equiv \sum _{\alpha =1}^{\infty }{q}^{\alpha }{M}^{\alpha \beta }=2\pi \int {\mathrm{d}}^2\varpi {\mathcal{F}}^{\prime }\sum _{n\ne 0}[\sum _{\alpha =1}^{\infty }{q}^{\alpha }{\mathrm{\Phi }}_n^{\alpha }(\mathit{\varpi })]{\mathrm{\Phi }}_n^{\beta }(\mathit{\varpi }),$$ $$\mathfrak{A}\equiv \sum _{\alpha =1}^{\infty }{q}^{\alpha }{A}^{\alpha }=-2\pi \int {\mathrm{d}}^2\mathit{\varpi }{f}_0(\mathit{\varpi })[\sum _{\alpha =1}^{\infty }{q}^{\alpha }{\mathrm{\Phi }}_0^{\alpha }(\mathit{\varpi })].$$