Open Access
Volume 684, April 2024
Article Number A58
Number of page(s) 10
Section Galactic structure, stellar clusters and populations
Published online 04 April 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 and , respectively. This represents a dipole shift perturbation corresponding to the spherical harmonic Pl = 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

(1) (2)

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

and the distribution function is given by


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


The dimensionless function, ℱ(ℰ), is normalized by the condition ∫ℱ d3r d3v = 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


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


for . 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 ℓ = ℓ0 + δℓ and expanding Eqs. (1)–(3) and (4) in Taylor series over ε, we obtain

(7) (8) (9)

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

(10) (11) (12)

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


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


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, fn, specifying the initial DF in the form of Eq. (14) and the potential in the form of Eq. (10), we should obtain ∂fn/∂t = 0. The test was conducted in PS23 using model (6). It confirmed the conservation of all harmonics, fn, 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, …


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


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


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


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, ρ, Φ ∝ Pl(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

(19) (20)

where .

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.

thumbnail 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


we have from Eq. (18)


or after substitution, r = cotψ, where ,


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


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/r3 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)


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, H0 = −ℰ, depends only on the action variables J = (Jr, Jθ, Jφ) in the form H0 = H0(Jr, L)= − ℰ, where L = Jθ + |Jφ| is the angular momentum, and the radial action Jr is


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


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


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

The relationship between ℰ, Jr, and L allows us to use either the pair (ℰ, L) or (Jr, 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 Jr; 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, .

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:


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:



(31) (32)

Further, ϖ ≡ (ℰ,L) stands for the set of variables ℰ and L. The integration ∫d2ϖ (…) 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 <  L <  ∞, and the line of circular orbits ; that is,


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:


where , − πξπ, and − πwπ. From Eq. (30), we have


and for the density perturbation,


Multiplying both sides of Eq. (36) by Φα(r), integrating over the volume d3r = 4πr2dr, and taking into account that d3r d3v = d3J d3w = (2π)2d2ϖ dw, we find


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


By substituting the expression for Φn with the expansion derived from Eq. (15), , we get an infinite homogeneous system of linear equations3

(39) (40)

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


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))


with the matrix


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:


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 ), and then subtracting the contribution of the zero harmonic numerically:



(46) (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 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


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, fn = 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 fn = 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 fn = 0, namely , while the correct expression for fn = 0, as can be seen from Eq. (14) in the case of our stationary dilation mode, is


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 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 fn = 0:


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 fn ≠ 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, , 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),


and that d3r d3v = d3J d3w = (2π)2d2ϖ dw, we get


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


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


or, taking into account that ,


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


where is defined by Eq. (46), and


Appendix A demonstrates how to obtain simple expressions for and without resorting to action-angle variables.

In the double sum over β and n, denoted as in Eq. (56), we isolate the term where n = 0. Since , according to Eq. (47), the contribution of this term to the left-hand side of Eq. (56) is


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




By utilizing the expression (Eq. (49)) for the zero harmonic fn = 0 and the expression (57) for , we can formulate the final system of inhomogeneous linear equations (Eq. (59)) for the coefficients, Cα, as follows:

(61) (62) (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, fn = 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, , and then sum over α, we obtain the equation


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


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


We denote by 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:


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


Therefore, a correct test for the matrix equation in the case of stationary perturbations, ω = 0, should consist of confirming the tendency of to approach zero when Nα tends to infinity:


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

thumbnail 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, CNα = 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.

thumbnail 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 fn = 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, fn = 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

This results in the emergence of the zero harmonic fn = 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−1. 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−1, 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.


In the context of this paper, we don’t need the second-order perturbations Φ2, ρ2, and f2. Their explicit expressions can be found in the cited work PS23. The Taylor expansion mentioned above applies to all relevant models including described in PS23, not just the isochrone model used here.


The purpose of introducing a biorthonormal basis is to avoid numerical integration of the radial part of the Poisson equation by decomposing the potential and density into basis functions (e.g., Clutton-Brock 1973; Polyachenko & Shukhman 1981; Lilley et al. 2018).


Here, lowercase latin and uppercase greek indices correspond to the Fourier harmonic numbers and the basis function numbers, respectively.


Formally, it is obtained by summing Nα − 1 equations of the system (Eq. (67)), with the weight, qα, and discarding the contribution of on the left-hand side and on the right-hand side of the resulting total equation. As can be seen from Eq. (66), these contributions are small when Nα ≫ 1. It’s easy to see that retaining these contributions in the corresponding parts of Eq. (68) would result in a linear dependence of the equations in this system for any Nα, and therefore the determinant would be zero for any Nα. However, such a trivial nullification of the determinant cannot be considered to be a valid test.


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


  1. Bertin, G. 2014, Dynamics of Galaxies, 2nd edn. (Cambridge University Press) [Google Scholar]
  2. Bertin, G., Pegoraro, F., Rubini, F., & Vesperini, E. 1994, ApJ, 434, 94 [Google Scholar]
  3. Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
  4. Clutton-Brock, M. 1973, Ap&SS, 23, 55 [Google Scholar]
  5. Fridman, A. M., & Polyachenko, V. L. 1984, Physics of Gravitating Systems. I - Equilibrium and stability (New York: Springer) [Google Scholar]
  6. Gradshteyn, I. S., & Ryzhik, I. M. 2015, in Table of Integrals, Series, and Products, 8th edn., eds. D. Zwillinger, & V. Moll (Amsterdam: Academic Press) [Google Scholar]
  7. Hénon, M. 1960, Ann. Astrophys., 23, 474 [NASA ADS] [Google Scholar]
  8. Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
  9. Jaffe, W. 1983, MNRAS, 202, 995 [NASA ADS] [Google Scholar]
  10. Kalnajs, A. J. 1971, ApJ, 166, 275 [CrossRef] [Google Scholar]
  11. Kalnajs, A. J. 1976, ApJ, 205, 751 [Google Scholar]
  12. Kalnajs, A. J. 1977, ApJ, 212, 637 [Google Scholar]
  13. King, I. R. 1966, AJ, 71, 64 [Google Scholar]
  14. Lau, J. Y., & Binney, J. 2021, MNRAS, 507, 2241 [Google Scholar]
  15. Lilley, E. J., Sanders, J. L., Evans, N. W., & Erkal, D. 2018, MNRAS, 476, 2092 [Google Scholar]
  16. Michie, R. W., & Bodenheimer, P. H. 1963, MNRAS, 126, 269 [Google Scholar]
  17. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  18. Polyachenko, E. V. 2004, MNRAS, 348, 345 [Google Scholar]
  19. Polyachenko, E. V. 2005, MNRAS, 357, 559 [Google Scholar]
  20. Polyachenko, V. L., & Shukhman, I. G. 1981, Sov. Ast., 25, 533 [Google Scholar]
  21. Polyachenko, E. V., & Shukhman, I. G. 2015, MNRAS, 451, 601 [Google Scholar]
  22. Polyachenko, E. V., & Shukhman, I. G. 2023, Astron. Rep., 67, 1156, [arXiv:2311.05551] (PS23) [Google Scholar]
  23. Polyachenko, E. V., Shukhman, I. G., & Borodina, O. I. 2021, MNRAS, 503, 660 [Google Scholar]
  24. Ramond, P., & Perez, J. 2021, J. Math. Phys., 62, 112704D [Google Scholar]
  25. Saha, P. 1991, MNRAS, 248, 494 [Google Scholar]
  26. Tremaine, S. 2005, ApJ, 625, 143 [Google Scholar]
  27. Van Kampen, N. G. 1955, Physica, 21, 949 [Google Scholar]
  28. Weinberg, M. D. 1991, ApJ, 368, 66 [Google Scholar]
  29. Weinberg, M. D. 1994, ApJ, 421, 481 [Google Scholar]

Appendix A: Derivation of matrix

It is possible to obtain the matrix , which includes all the Fourier harmonics, einw, 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


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


Substitution of Φ(r) in the form of series


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


or, in a compact form,



(A.6) (A.7)

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

(A.8) (A.9)


(A.10) (A.11)

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

(A.12) (A.13)

where . Next, we go in (A.6) and (A.7) from variable r to variable ψ:


and find the final expressions for and in the form of simple single integrals

(A.15) (A.16)

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

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


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


we obtain for their derivatives


and finally find I1 in the form (A.12).

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


Again, using (A.18) and (A.19) yields I2 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

and sum over all α. We get



One can be sure that


This follows from the relation


which can be derived from the equation (1.444.7) of the book by Gradshteyn & Ryzhik (2015)


if we differentiate both of its parts with respect to ψ and use the expression for Φα(r) in the form (A.14). From (B.3) we have, multiplying both sides by einw and integrating over w,


Therefore 𝔄 = (2π)2∫d2ϖf0(ϖ)=δM = 0, where δM is the perturbed mass. Using the explicit form fn = 0(ϖ) (50) it is easy to verify analytically that the integral ∫d2ϖf0(ϖ) vanishes. From (B.5) it follows that the sum , included in 𝔐β, also vanishes, since this sum does not contain a term with n = 0. So 𝔐β also vanishes. As a result, (B.1) turns into the equation

that is, equation (25).

All Figures

thumbnail 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 the text
thumbnail 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.

In the text
thumbnail 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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.