Issue 
A&A
Volume 684, April 2024



Article Number  A58  
Number of page(s)  10  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202348556  
Published online  04 April 2024 
Exploring the dynamics of collisionless spherical stellar systems using the matrix method: Insights from the dilation mode
^{1}
Zentrum für Astronomie der Universität Heidelberg, Astronomisches RechenInstitut, Mönchhofstr 1214, 69120 Heidelberg, Germany
email: epolyachenko@uniheidelberg.de
^{2}
Institute of Astronomy Russian Academy of Sciences, 48 Pyatnitskaya st., 119017 Moscow, Russia
^{3}
Institute of SolarTerrestrial Physics, Russian Academy of Sciences, Siberian Branch, PO Box 291 Irkutsk 664033, Russia
email: shukhman@iszf.irk.ru
Received:
10
November 2023
Accepted:
18
December 2023
Context. Analytical solutions to the perturbed equations that govern selfgravitating collisionless stellar systems are crucial for both code testing and theoretical insights. For spheres, a solution has been known for years that corresponds to the entire object’s shift from the origin. We recently introduced a new exact stationary solution, relevant for models with a single length parameter. This solution, referred to as the scaleinvariant or dilation mode, has led to insights regarding the concept of perturbation energy within the linear theory framework.
Aims. Our aim is to use Hénon’s isochrone model as an example to verify the ability of the standard matrix method to successfully predict the existence of a dilation mode, and to explore its potential application as a test disturbance.
Methods. We used the standard matrix method for radial perturbations and applied CluttonBrock potentialdensity pairs to determine the properties of the perturbations.
Results. In this particular case of stationary radial perturbations, the typical relationship between the perturbations of the distribution function and the potential fails. This discrepancy poses a challenge when attempting to use the dilation mode as a test. When using CluttonBrock pairs with the matrix method, a mass conservation equation as an additional equation to the ordinary set of linear equations is required. With this added equation, it’s possible to obtain the needed test: identical vanishing of the determinant of this modified set of equations with an increasing number of included basis functions.
Key words: galaxies: clusters: general / galaxies: kinematics and dynamics / galaxies: star clusters: general
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), 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 wellknown 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 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 selfconsistent 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 secondorder quantity in amplitude, which is not calculable within linear theory. We’ve shown that the accurate expression for perturbation energy, constructed with secondorder perturbations, doesn’t coincide with the wellknown expression for perturbation energy, constructed as a bilinear form of the firstorder perturbations, as is given in Sect. 5.4.2 of Binney & Tremaine (2008). These two forms of secondorder 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 wellestablished 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 halfplane, ω, 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
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 ∫ℱ d^{3}r d^{3}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 selfconsistent 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
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:
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 f_{1}. Using Eq. (10), we can rewrite Eq. (12) as
The perturbations of firstorder Φ, ρ, 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. CluttonBrock 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 potentialdensity 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 symbol^{2}. 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 CluttonBrock (1973) functions (hereafter CB) as biorthonormal potentialdensity 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
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.
Fig. 1. Basis functions of potential for α = 1 and α = 10 from the CluttonBrock 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/r^{3} 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, H_{0} = −ℰ, depends only on the action variables J = (J_{r}, J_{θ}, J_{φ}) in the form H_{0} = H_{0}(J_{r}, L)= − ℰ, where L = J_{θ} + J_{φ} is the angular momentum, and the radial action J_{r} is
The radial action and angular momentum are conserved along the unperturbed orbit. The radial frequency that corresponds to J_{r} is
It should be noted that in the isochrone model an analytic expression for Hamiltonian H_{0}(J_{r}, 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 ℰ, 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 actionangle 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 actionangle 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}:
Here,
Further, ϖ ≡ (ℰ,L) stands for the set of variables ℰ and L. The integration ∫d^{2}ϖ (…) 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 d^{3}r = 4πr^{2}dr, and taking into account that d^{3}r d^{3}v = d^{3}J d^{3}w = (2π)^{2}d^{2}ϖ 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 equations^{3}
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 twodimensional 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 semifinished matrix ), and then subtracting the contribution of the zero harmonic numerically:
where
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, 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 , while the correct expression for f_{n = 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 f_{n = 0}:
The nonvanishing 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 righthand 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, , 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 d^{3}r d^{3}v = d^{3}J d^{3}w = (2π)^{2}d^{2}ϖ dw, we get
Expanding f in Eq. (52) into harmonics (see Eqs. (31) and (32)), we obtain
Using Eq. (14) for the harmonics, f_{n}(ϖ), we have
or, taking into account that ,
where δ_{nm} is the Kronecker delta. Substituting Eq. (55) into the righthand 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 actionangle 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 lefthand side of Eq. (56) is
By shifting it to the righthand side of Eq. (56), we obtain
where
By utilizing the expression (Eq. (49)) for the zero harmonic f_{n = 0} and the expression (57) for , we can formulate the final system of inhomogeneous linear equations (Eq. (59)) for the coefficients, C^{α}, as follows:
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 lefthand 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 nontrivial 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 equation^{4}
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^{α}.
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 CluttonBrock (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 selfconsistent model but with adjusted energy and scale length, its new DF is
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^{−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 nonradial 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 secondorder perturbations Φ_{2}, ρ_{2}, and f_{2}. 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., CluttonBrock 1973; Polyachenko & Shukhman 1981; Lilley et al. 2018).
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 lefthand side and on the righthand 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.
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
 Bertin, G. 2014, Dynamics of Galaxies, 2nd edn. (Cambridge University Press) [Google Scholar]
 Bertin, G., Pegoraro, F., Rubini, F., & Vesperini, E. 1994, ApJ, 434, 94 [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
 CluttonBrock, M. 1973, Ap&SS, 23, 55 [Google Scholar]
 Fridman, A. M., & Polyachenko, V. L. 1984, Physics of Gravitating Systems. I  Equilibrium and stability (New York: Springer) [Google Scholar]
 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]
 Hénon, M. 1960, Ann. Astrophys., 23, 474 [NASA ADS] [Google Scholar]
 Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
 Jaffe, W. 1983, MNRAS, 202, 995 [NASA ADS] [Google Scholar]
 Kalnajs, A. J. 1971, ApJ, 166, 275 [CrossRef] [Google Scholar]
 Kalnajs, A. J. 1976, ApJ, 205, 751 [Google Scholar]
 Kalnajs, A. J. 1977, ApJ, 212, 637 [Google Scholar]
 King, I. R. 1966, AJ, 71, 64 [Google Scholar]
 Lau, J. Y., & Binney, J. 2021, MNRAS, 507, 2241 [Google Scholar]
 Lilley, E. J., Sanders, J. L., Evans, N. W., & Erkal, D. 2018, MNRAS, 476, 2092 [Google Scholar]
 Michie, R. W., & Bodenheimer, P. H. 1963, MNRAS, 126, 269 [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
 Polyachenko, E. V. 2004, MNRAS, 348, 345 [Google Scholar]
 Polyachenko, E. V. 2005, MNRAS, 357, 559 [Google Scholar]
 Polyachenko, V. L., & Shukhman, I. G. 1981, Sov. Ast., 25, 533 [Google Scholar]
 Polyachenko, E. V., & Shukhman, I. G. 2015, MNRAS, 451, 601 [Google Scholar]
 Polyachenko, E. V., & Shukhman, I. G. 2023, Astron. Rep., 67, 1156, [arXiv:2311.05551] (PS23) [Google Scholar]
 Polyachenko, E. V., Shukhman, I. G., & Borodina, O. I. 2021, MNRAS, 503, 660 [Google Scholar]
 Ramond, P., & Perez, J. 2021, J. Math. Phys., 62, 112704D [Google Scholar]
 Saha, P. 1991, MNRAS, 248, 494 [Google Scholar]
 Tremaine, S. 2005, ApJ, 625, 143 [Google Scholar]
 Van Kampen, N. G. 1955, Physica, 21, 949 [Google Scholar]
 Weinberg, M. D. 1991, ApJ, 368, 66 [Google Scholar]
 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, e^{inw}, quite simply, without passing to the actionangle 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 righthand side of (A.2) yields
or, in a compact form,
where
Since the integrands (A.6) and (A.7) depend only on ℰ and r, we proceed to integration over these variables:
where
The integrals I_{1, 2}(r) are taken analytically. Below, we show how to take them. We give the result,
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
Finally, we show how the integrals, I_{1, 2}(r), are calculated.
(i) Calculation of I_{1}(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 I_{1} in the form (A.12).
(ii) Calculation of I_{2}(r). Integrating by parts, we get
Again, using (A.18) and (A.19) yields I_{2} 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
where
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 e^{−inw} and integrating over w,
Therefore 𝔄 = (2π)^{2}∫d^{2}ϖ f_{0}(ϖ)=δM = 0, where δM is the perturbed mass. Using the explicit form f_{n = 0}(ϖ) (50) it is easy to verify analytically that the integral ∫d^{2}ϖ f_{0}(ϖ) 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
Fig. 1. Basis functions of potential for α = 1 and α = 10 from the CluttonBrock set. The number of nodes of Φ^{α}(r) is equal to α − 1. 

In the text 
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 
Fig. 3. Expansion coefficients, C^{α}, of the perturbed potential Eq. (21) over CluttonBrock (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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.