1 INTRODUCTION

One of the traditional methods for studying the perturbation dynamics of equilibrium models of spherical stellar systems is to study the evolution of small perturbations. As a rule, the main question of interest to researchers is whether the equilibrium state, described by the distribution function (DF) \(F({\mathbf{r}},{\mathbf{v}})\) of stars and the gravitational potential \({{\Phi }_{0}}({\mathbf{r}})\), is stable or unstable. Along with methods consisting in finding general stability criteria by deriving the corresponding theorems (see, e.g., monograph [1], hereinafter BT, and references therein), there is a method for solving a linearized eigenvalue problem. For this, assuming that the perturbations of the gravitational potential \(\Phi ({\mathbf{r}},t)\) and DF \(f({\mathbf{r}},{\mathbf{v}},t)\) are small and proportional to \(\exp ( - i\omega t)\), the eigenvalues ω of the linearized system of equations, consisting of the collisionless kinetic equation and the Poisson equation, are found. The presence of eigenvalues with \({\text{Im}}(\omega ) > 0\) means that the system is unstable.

Finding the eigenvalues \(\omega \) is a rather computationally expensive task. Except for a few models, where the equilibrium potential is harmonic (see, e.g., [2‒5]), it is solved using the so-called matrix methods. In these methods, the problem reduces to numerically finding the roots \(\omega \) of a certain determinant \(\mathcal{D}(\omega ) \equiv \det \left\| {{{D}^{{\alpha \beta }}}(\omega )} \right\|\) = 0, \(\alpha ,\beta = 1,2,3...\). For disk models, the matrix method was first proposed by Kalnajs [6], and for spherical systems, by Polyachenko and Shukhman [7]. It consists in expanding the amplitudes \(\hat {\Phi }({\mathbf{r}})\) and \(\hat {\rho }({\mathbf{r}})\) of the perturbed potential and density \(\Phi ({\mathbf{r}},t) = \hat {\Phi }({\mathbf{r}}){\kern 1pt} {{e}^{{ - i\omega t}}}\) and \(\rho ({\mathbf{r}},t) = \hat {\rho }({\mathbf{r}}){\kern 1pt} {{e}^{{ - i\omega t}}}\) in the so-called biorthonormal set of potential–density basis pairs \({{\Phi }^{\alpha }}(r)\) and \({{\rho }^{\alpha }}(r)\) and obtaining a system of linear equations for the expansion coefficients \({{C}^{\alpha }}\). The equality of the determinant of this system to zero leads to the desired dispersion relation. This method works for systems with an integrable Hamiltonian \({{H}_{0}}\), i.e., for equilibrium stellar systems, the potential \({{\Phi }_{0}}({\mathbf{r}})\) of which allows a conversion from the coordinate–velocity variables \(({\mathbf{r}},{\mathbf{v}})\) to the action–angle variables \(({\mathbf{J}},{\mathbf{w}})\). In the alternative matrix method proposed by E. Polyachenko (see [8, 9]), the original system of linearized equations reduces to a standard linear eigenvalue problem of the form \(\omega {{f}_{n}}({\mathbf{J}}) = \sum\nolimits_{n{\kern 1pt} '}^{} {\int \,d{\mathbf{J}}{\kern 1pt} '{{K}_{{nn{\kern 1pt} '}}}({\mathbf{J}},{\mathbf{J}}{\kern 1pt} '){\kern 1pt} {{f}_{{n{\kern 1pt} '}}}({\mathbf{J}}{\kern 1pt} ')} \), where \({{f}_{n}}({\mathbf{J}})\) are the harmonics of the Fourier expansion of the perturbed DF in angular variables \({\mathbf{w}}\) and \({{K}_{{nn{\kern 1pt} '}}}({\mathbf{J}},{\mathbf{J}}{\kern 1pt} ')\) is the kernel.

Recently, a number of works have appeared [10–13], studying the dynamics of perturbations against the background of equilibrium models not with the aim of studying their stability, as happened over the previous decades and was reflected in numerous works and monographs (see, e.g., [1, 14, 15] and references therein), but with the aim of studying fluctuations of density and potential around equilibrium and their effect on the processes of slow relaxation, as well as their role in the process of N-body numerical simulation of processes in stellar systems. The equilibrium models under their consideration are knowingly stable systems, and perturbations in them can be caused either by the noise associated with a finite number \(N\) of particles [12] or by an external action. In this case, weakly damped oscillations are of interest, which can last for many characteristic crossing times, being almost identical to real neutral eigenmodes [16, 17].

In stable equilibrium spherical systems, there are no discrete modes with Im\((\omega ) > 0\). Due to the reversibility of the collisionless kinetic equation in time, there are no damped discrete modes, with Im\((\omega ) < 0\). The presence of discrete neutral modes, with Im\((\omega ) = 0\), is possible only in rare situations. This is due to the presence of resonances of disturbance waves with the orbital motion of stars of the type \(\omega - {\mathbf{n}}\Omega = 0\), where \(\Omega ({\mathbf{J}}) = ({{\Omega }_{1}},{{\Omega }_{2}},{{\Omega }_{3}})\) are the frequencies of the orbital motion and \({\mathbf{n}} = ({{n}_{1}},{{n}_{2}},{{n}_{3}})\) are integer numbers; therefore, neutral discrete modes are possible only in the presence of “gaps” in the phase space that are free from resonance.Footnote 1 It turns out that, for such equilibrium models, the complete system of eigenmodes is represented exclusively by the continuous spectrum of van Kampen modes [19] with a real frequency \(\omega \). Note that a perturbation that decays exponentially according to the so-called Landau damping [20], in which the frequency \({{\omega }_{{\text{L}}}}\) has a negative imaginary part, \({{\omega }_{{\text{L}}}} = {\text{Re}}({{\omega }_{{\text{L}}}}) + i{\kern 1pt} {\text{Im}}({{\omega }_{{\text{L}}}})\), \({\text{Im}}({{\omega }_{{\text{L}}}}) < 0\), is not a true damped eigenmode, but represents a continuum superposition of singular van Kampen modes. To distinguish a Landau-damped perturbation from a true eigenmode, we will call it a quasi-mode. We traced in more detail the dynamics of initial perturbations in the form of a superposition of van Kampen modes and its relation with Landau quasi-modes for infinite homogeneous gravitating systems in [21] and, for the case of shear fluid flows, in [22].

The presence of weakly damped Landau quasi-modes plays a significant role for stable systems. In terms of van Kampen modes, their presence means that their amplitude is especially large when the van Kampen mode frequencies \(\omega \) are close to the real part of the Landau quasi-mode frequency: \(\omega \approx {\text{Re}}({{\omega }_{{\text{L}}}})\). On the other hand, their presence allows oscillations excited, say, by the close passage of a disturbing external body, to last for a rather long time practically undamped [16].

Therefore, the search for Landau quasi-modes for stable systems is of interest. However, from a practical point of view, finding Landau quasi-modes poses significant difficulties. The fact is that the dispersion relation obtained by any of the matrix methods described above [7–9] holds only in the upper half-plane of the complex variable \(\omega \), while the frequencies of Landau quasi-modes lie in the lower half-plane of \(\omega \). This is related to the principle of causality and was repeatedly described in the literature, starting with the pioneering work of Landau [20] (see also BT [1]). In order to use the dispersion relation \(\mathcal{D}(\omega ) = 0\) to find the frequencies of Landau quasi-modes, it is necessary to perform an analytical continuation of the function \(D(\omega )\) into the lower half-plane of the complex variable \(\omega \). Landau [20] was the first to perform this procedure for a homogeneous electron plasma. To do this, he deformed the integration contour over the (only one in his problem) variable velocity \(v\), shifting it into the complex plane of \(v\) so that it passed below all possible resonance points \({{v}_{c}} \equiv {{\omega }_{{\text{L}}}}{\text{/}}k\). This procedure is called the Landau–Lin bypass rule, since Lin [23] derived the same bypass rule for shear flows of an inviscid fluid, but based not on the principle of causality (meaning that the perturbation must disappear in the distant past), like Landau, but on the principle of dissipativity (i.e., by adding an infinitesimal positive viscosity to the inviscid Euler equation).

The problem of finding the analytical continuation \(\mathcal{D}(\omega )\) for equilibrium spherical stellar systems is much more difficult than in a homogeneous plasma [20], in an infinite homogeneous gravitating medium [21], or in shear fluid flows [22]. Firstly, even in the simplest case of self-consistent models, we are dealing with at least a two-dimensional phase space in the variables of action \({\mathbf{J}}\) rather than with a one-dimensional one, where we have to work with integrals containing only a single velocity component, parallel to a fixed direction of the wave vector \({\mathbf{k}}\). The second problem making the analytical continuation more difficult is associated with the presence in the integrals over the phase volume of not a single resonant denominator of the form \(1{\text{/}}(\omega - kv)\), but an infinite number of them of the form \(1{\text{/}}[\omega - {\mathbf{n}} \cdot {\mathbf{\Omega }}{\kern 1pt} ({\mathbf{J}})]\). The first of these problems was overcome by Barré et al. [24] by considering a one-dimensionally nonuniform system with an artificial, rather simple one-dimensional interaction potential between particles (not gravitational), reducing the problem to a one-dimensional one, albeit with a large number of resonant denominators of the type \(1{\text{/}}[\omega - n{\kern 1pt} \Omega (J)]\).

For spherical systems with a real gravitational potential (more precisely, for King’s models [25]), an attempt to construct an analytic continuation of the determinant \(\mathcal{D}(\omega )\) into the lower half-plane was made by Weinberg [16]. To do this, he approximated the function \(\mathcal{D}(\omega )\) in the upper half-plane by a sum of fractional-rational functions that admit a simple analytic continuation to the lower half-plane and, having obtained an approximate expression for the analytic continuation of \(\mathcal{D}(\omega )\), found (for certain model parameters) the frequencies of weakly damped Landau quasi-modes. Although the results of this work are widely cited in the literature, from our point of view, they are not convincing enough.

Another way to find Landau exponential decay is to directly solve the evolution equation (more precisely, a system of equations) for the Fourier harmonics of the perturbed DF \({{f}_{n}}(J,t)\). To do this, it is necessary to specify the initial DF \(f({\mathbf{J}},0)\) and the corresponding perturbed potential \(\Phi ({\mathbf{r}},0)\). If the equilibrium state under consideration contains a Landau quasi-mode, it must manifest itself for any choice of the initial DF, since the determinant \(\mathcal{D}(\omega )\) depends only on the properties of the unperturbed system, and the presence of zeros in its lower half-plane means that, asymptotically, perturbations of the density and potential must decay exponentially. This follows from the fact that the zeros of \(\mathcal{D}(\omega )\) are poles of the Laplace image of the perturbation. The corresponding procedure for solving the evolution equation for an infinite homogeneous medium [21] and for shear fluid flows [22] was performed explicitly and demonstrated full agreement of the asymptotic behavior of the amplitude \({{\hat {\rho }}_{k}}(t)\) of stellar density perturbations or the amplitude of total vorticity across the channel, \({{N}_{k}}(t) = \int \,dy{\kern 1pt} {{\hat {\zeta }}_{k}}(y,t)\), with the Landau damping with a frequency \({{\omega }_{{\text{L}}}}\) found from the condition \(\mathcal{D}(\omega ) = 0\).

The aforesaid means that the problem of numerical study of the perturbation in stable systems depends rather strongly on the selected codes and calculation parameters: the grid on the phase plane, the number of retained Fourier harmonics in the variables \({\mathbf{w}}\), and the number of retained basis functions. Therefore, the presence of a test perturbation for code verification is highly desirable. One such test perturbation has long been known. It consists in a shift of the spherical system as a whole. If this shift occurs, say, along the \(z\) axis by a small distance \(\xi \), then the perturbation in the density and the potential that arise in this case are \(\rho (r,\theta ) = - \xi {\kern 1pt} \rho _{0}^{'}(r){\kern 1pt} \cos \theta \) and \(\Phi (r,\theta ) = - \xi {\kern 1pt} \Phi _{0}^{'}(r){\kern 1pt} \cos \theta \), respectively. This is a dipole shear perturbation corresponding to a spherical harmonic \({{P}_{{l = 1}}}(\cos \theta ) = \cos \theta \). Obviously, the natural frequency \(\omega \) corresponding to this perturbation is zero. This test has been repeatedly used previously when verifying codes in the analysis of stability (see, e.g., [26–28]).

In this paper, we propose another simple test perturbation, which admits an exact solution. It works for a very specific class of spherical models, namely, for models described by DFs that depend only on energy \(E\) and contain a single length parameter \(b\).

This exact solution allows us to simultaneously clarify another issue concerning the correct definition of the concept of perturbation energy. The fact is that the perturbation energy, being a quantity quadratic in amplitude, at first glance, in principle, cannot be calculated in linear theory. However, it can be shown that the system of linearized kinetic equation and Poisson’s equation admits a quadratic integral of motion, whose form is very similar to the total perturbation energy, which, strictly speaking, is not the same as the true energy, since its calculation requires knowledge of the perturbations of the DF, the potential, and the density of the second order. The test perturbation proposed allows one to calculate the perturbation energy in any order and perform a comparison of these two second-order “energies”.

In Section 2, we present the idea of a test perturbation and give several examples of models of self-consistent equilibrium DFs, whose spectrum contains the considered mode. In Section 3, we will consider in more detail the concept of perturbation energy, which can be constructed within linear theory, and, using a scale-invariant perturbation as an example, compare the accurately calculated energy (taking into account the second-order perturbations) with the generally accepted expression for the perturbation energy in li-near theory. The results obtained are discussed in S-ection 4.

2 THE IDEA OF A TEST PERTURBATION AND SEVERAL EXAMPLES OF RELEVANT DF MODELS

Let the spherical model be described by an equilibrium DF containing a single characteristic scale in the radial variable \(r\). It will be called the scale factor and denoted by \(b\). For such models, the unperturbed potential and density have the form

$$\begin{gathered} {{\Phi }_{0}}(r,b) = \frac{{MG}}{b}{\kern 1pt} \phi \left( {\frac{r}{b}} \right), \\ {{\rho }_{0}}(r,b) = \frac{M}{{{{b}^{3}}}}{\kern 1pt} \varrho \left( {\frac{r}{b}} \right), \\ \end{gathered} $$
(1)

where \(\phi (x)\) and \(\varrho (x)\) are related by Poisson’s equation:

$$\frac{1}{{{{x}^{2}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{d}{{dx}}\left[ {{{x}^{2}}{\kern 1pt} \frac{{d\phi (x)}}{{dx}}} \right] = 4\pi {\kern 1pt} \varrho (x).$$
(2)

The distribution function is

$$\begin{gathered} {{F}_{0}}(\mathcal{E},b) = \frac{1}{{{{{(MG{\kern 1pt} b)}}^{{3/2}}}}}{\kern 1pt} \mathcal{F}(\mathcal{E}), \\ 0 \leqslant \mathcal{E} \leqslant \Psi (0) \equiv - \phi (0), \\ \end{gathered} $$
(3)

where the dimensionless energy \( - \mathcal{E}\) should also be regarded as a function of \(v\), \(r\), and the scale factor \(b\):

$$\begin{gathered} \mathcal{E} = \mathcal{E}(r,v;b) = - \frac{b}{{MG}}\left[ {\frac{1}{2}{{v}^{2}} + {{\Phi }_{0}}(r,b)} \right] \\ = - \frac{b}{{MG}}\left[ {\frac{1}{2}{{v}^{2}} + \frac{{MG}}{b}\phi \left( {\frac{r}{b}} \right)} \right]. \\ \end{gathered} $$
(4)

Here, the dimensionless function \(\mathcal{F}(\mathcal{E})\) is normalized so that \(\int \,\mathcal{F}{\kern 1pt} {{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{{{\kern 1pt} 3}}}{\mathbf{v}} = 1\).

It is absolutely clear that, if we fix the total mass \(M\), but change \(b\), we will get the same equilibrium model, but only with a different scale factor, \(b + \delta b\). But this means that the eigenfrequency \(\omega \) of the mode corresponding to such expansion/contraction is equal to zero. This fact can serve as a test of various codes when studying the perturbation dynamics in spherical systems.

Let us give several examples of models that have this form.

2.1 Isochrone Hénon model [29]

For it, the function \(\phi (x)\) entering into the potential is

$$\phi (x) = - \frac{1}{{1 + a}}{\kern 1pt} ,\quad a = \sqrt {1 + {{x}^{2}}} {\kern 1pt} ,$$
(5)

the density is

$$\varrho (x) = \frac{1}{{4\pi }}{\kern 1pt} {\kern 1pt} \frac{{1 + 2{\kern 1pt} a}}{{{{{(1 + a)}}^{2}}{\kern 1pt} {{a}^{3}}}}{\kern 1pt} ,$$
(6)

and the distribution function is

$$\begin{gathered} {{\mathcal{F}}_{{{\text{He}}{\kern 1pt} '{\text{non}}}}}(\mathcal{E}) = \frac{1}{{\sqrt{2} {\kern 1pt} {{{(2\pi )}}^{3}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\sqrt{\mathcal{E}} }}{{{{{[2{\kern 1pt} (1 - \mathcal{E})]}}^{4}}}}\left[ {64{{\mathcal{E}}^{4}} - 240{{\mathcal{E}}^{3}}} \right. \\ \, + 320{{\mathcal{E}}^{2}} - 66\mathcal{E} + 27 + 3{\kern 1pt} (16{{\mathcal{E}}^{2}} \\ \, + 28\mathcal{E} - 9)\left. {\frac{{\arcsin \sqrt{\mathcal{E}} }}{{\sqrt{\mathcal{E}{\kern 1pt} (1 - \mathcal{E})} }}} \right],\quad 0 \leqslant \mathcal{E} \leqslant \frac{1}{2}. \\ \end{gathered} $$
(7)

2.2 Hernquist model [30]

The potential for it is

$$\phi (x) = - \frac{1}{{1 + x}}{\kern 1pt} ,$$
(8)

the density is

$$\varrho (x) = \frac{1}{{2\pi }}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{x{\kern 1pt} {{{(1 + x)}}^{3}}}}{\kern 1pt} ,$$
(9)

and the distribution function is

$$\begin{gathered} {{\mathcal{F}}_{{{\text{Hernquist}}}}}(\mathcal{E}) = \frac{1}{{\sqrt{2} {\kern 1pt} {{{(2\pi )}}^{3}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\sqrt{\mathcal{E}}}}{{{{{(1 - \mathcal{E})}}^{2}}}} \\ \times \;\left[ {(1 - 2\mathcal{E}){\kern 1pt} (8{{\mathcal{E}}^{2}} - 8\mathcal{E} - 3) + \frac{{3\text{arcsin}\sqrt{\mathcal{E}}}}{{\sqrt{\mathcal{E}{\kern 1pt} (1 - \mathcal{E})} }}} \right], \\ 0 \leqslant \mathcal{E} \leqslant 1. \\ \end{gathered} $$
(10)

2.3 Jaffe model [31]

The potential for it is

$$\phi (x) = - \ln \left( {1 + \frac{1}{x}} \right),$$
(11)

the density is

$$\varrho (x) = \frac{1}{{4\pi }}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{{{x}^{2}}{\kern 1pt} {{{(1 + x)}}^{2}}}}{\kern 1pt} ,$$
(12)

and the distribution function is

$$\begin{gathered} {{\mathcal{F}}_{{{\text{Jaffe}}}}}(\mathcal{E}) = \frac{1}{{2{{\pi }^{3}}}}[{{F}_{ - }}(\sqrt{2\mathcal{E}} ) - \sqrt{2} {\kern 1pt} {{F}_{ - }}(\sqrt{\mathcal{E}}) \\ - \;\sqrt{2} {\kern 1pt} {{F}_{ + }}(\sqrt{\mathcal{E}}) + {{F}_{ + }}(\sqrt{2\mathcal{E}} )],\quad 0 \leqslant \mathcal{E} < \infty , \\ \end{gathered} $$
(13)

where

$${{F}_{ \pm }}(x) = {{e}^{{ \mp {{x}^{2}}}}}\int\limits_0^x \,dy{\kern 1pt} {{e}^{{ \pm {{y}^{2}}}}}.$$
(14)

2.4 Plummer model [32]

The potential for it is

$$\phi (x) = - \frac{1}{{\sqrt {1 + {{x}^{2}}} }}{\kern 1pt} ,$$
(15)

the density is

$$\varrho (x) = \frac{3}{{4\pi }}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{{{{(1 + {{x}^{2}})}}^{{5/2}}}}}{\kern 1pt} ,$$
(16)

and the distribution function is

$$\begin{gathered} {{\mathcal{F}}_{{{\text{Plummer}}}}}(\mathcal{E}) = A{\kern 1pt} {{\mathcal{E}}^{{7/2}}}, \\ A = \frac{3}{7}{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{{2}^{7}}}}{{\sqrt 2 {\kern 1pt} {{{(2\pi )}}^{3}}}}{\kern 1pt} , \\ 0 \leqslant \mathcal{E} \leqslant \Psi (0) = 1. \\ \end{gathered} $$
(17)

2.5 Polytropes

Note that the Plummer model is a special case of a series of polytropic models with a distribution function

$${{\mathcal{F}}_{{{\text{polytropes}}}}}(\mathcal{E}) = {{A}_{n}}{\kern 1pt} {{\mathcal{E}}^{{n - 3/2}}}$$
(18)

and density

$$\begin{gathered} \varrho (x) = {{\Lambda }_{n}}{\kern 1pt} {{A}_{n}}{\kern 1pt} {{\Psi }^{n}}(x), \\ {{\Lambda }_{n}} = \frac{1}{{n!}}{\kern 1pt} {{(2\pi )}^{{3/2}}}{\kern 1pt} \Gamma \left( {n - \frac{1}{2}} \right), \\ \end{gathered} $$
(19)

corresponding to \(n = 5\). In the case of arbitrary \(n\), these models have no explicit analytical expression for the potential \(\phi (x)\), but there is a corresponding second-order nonlinear equation for the potential, following from Poisson’s equation (Lane–Emden equation). Analysis of this equation shows (see BT) that polytropic models with \(n > 5\) have infinite mass and are irrelevant. However, models with \(1{\text{/}}2 < n < 5\), which, although they have (unlike the models given above) a finite radius \(b\), must nevertheless also contain a scale-invariant mode in the spectrum, since this radius is the only length scale in the model. In particular, for \(n = 1\), when the Lane–Emden equation becomes linear, there is an analytical solution with finite radius and mass:

$$\phi (x) = - \frac{{\sin (\pi x)}}{{\pi {\kern 1pt} x}},\quad \varrho (x) = \frac{{\sin (\pi x)}}{{4{\kern 1pt} x}},\quad x \leqslant 1,$$
(20)
$$\mathcal{F}(\mathcal{E}) = \frac{{\sqrt 2 }}{{16\pi }}{{\mathcal{E}}^{{ - 1/2}}},\quad 0 \leqslant \mathcal{E} \leqslant 1.$$
(21)

Note that models with \(n < 3{\text{/}}2\) have a positive sign of the derivative with respect to the energy \(E = - \frac{{MG}}{b}\mathcal{E}\), i.e., \(\mathcal{F}{\kern 1pt} '(\mathcal{E}) < 0\) and, in principle, may be unstable. We will not discuss this issue in more detail here.

Consider the change in model parameters, associated with variations in the scale factor \(b\), \(b = {{b}_{0}} + \delta b\):

$$\begin{gathered} \Phi (r,b) = {{\Phi }_{0}}(r,{{b}_{0}}) + \epsilon {{\Phi }_{1}}(r,{{b}_{0}}) \\ + \;{{\epsilon }^{2}}{\kern 1pt} {{\Phi }_{2}}(r,{{b}_{0}}) + \mathcal{O}({{\epsilon }^{3}}), \\ \end{gathered} $$
(22)
$$\begin{gathered} \rho (r,b) = {{\rho }_{0}}(r,{{b}_{0}}) + \epsilon {\kern 1pt} {{\rho }_{1}}(r,{{b}_{0}}) \\ + \;{{\epsilon }^{2}}{\kern 1pt} {{\rho }_{2}}(r,{{b}_{0}}) + \mathcal{O}({{\epsilon }^{3}}), \\ \end{gathered} $$
(23)
$$\begin{gathered} F(\mathcal{E},b) = {{F}_{0}}(\mathcal{E},{{b}_{0}}) + \epsilon {{f}_{1}}(\mathcal{E},{{b}_{0}}) \\ + \;{{\epsilon }^{2}}{\kern 1pt} {{f}_{2}}(\mathcal{E},{{b}_{0}}) + \mathcal{O}({{\epsilon }^{3}}). \\ \end{gathered} $$
(24)

Here, \(\epsilon = \delta b{\text{/}}{{b}_{0}} \ll 1\) is the expansion parameter. Note that we have prepared an expansion of all quantities up to the second order. Although, linear theory does not require the knowledge of second-order quantities, we do this in order to obtain a correct expression for potential and kinetic energies, which are quantities of the second order in the perturbation amplitude and cannot be calculated simply as a bilinear form from first-order quantities. Setting \(G = M = 1\), we have the following expressions for the potential,

$$\begin{gathered} {{\Phi }_{1}} = - \frac{1}{{{{b}_{0}}}}(x{\kern 1pt} \phi ){\kern 1pt} ', \\ {{\Phi }_{2}} = \frac{1}{{2{\kern 1pt} {{b}_{0}}}}[{\kern 1pt} 2{\kern 1pt} (x{\kern 1pt} \phi ){\kern 1pt} '\; + x{\kern 1pt} (x{\kern 1pt} \phi ){\kern 1pt} ''], \\ \end{gathered} $$
(25)

the density,

$$\begin{gathered} {{\rho }_{1}} = - \frac{1}{{b_{0}^{3}}}({\kern 1pt} 3{\kern 1pt} \varrho + x{\kern 1pt} \varrho {\kern 1pt} '), \\ {{\rho }_{2}} = \frac{1}{{b_{0}^{3}}}\left( {6{\kern 1pt} \varrho + 4{\kern 1pt} x{\kern 1pt} \varrho {\kern 1pt} '\; + \frac{1}{2}{{x}^{2}}{\kern 1pt} \varrho {\kern 1pt} ''} \right), \\ \end{gathered} $$
(26)

and the distribution function,

$${{f}_{1}} = \frac{1}{{b_{0}^{{3/2}}}}\left\{ { - \frac{3}{2}{\kern 1pt} \mathcal{F} + \mathcal{F}{\kern 1pt} '(\mathcal{E}){\kern 1pt} [\mathcal{E} + (x{\kern 1pt} \phi ){\kern 1pt} ']} \right\},$$
(27)
$$\begin{gathered} {{f}_{2}} = \frac{1}{{2b_{0}^{{3/2}}}}\left\{ {\frac{{15}}{4}{\kern 1pt} \mathcal{F}} \right. - \mathcal{F}{\kern 1pt} '{\kern 1pt} \text{[}3{\kern 1pt} \mathcal{E} + 3{\kern 1pt} (x\phi ){\kern 1pt} ' + x{\kern 1pt} (x\phi ){\kern 1pt} ''] \\ \left. {\mathop + \limits_{_{{_{{}}}}} \;\mathcal{F}{\kern 1pt} '\,{{{[\mathcal{E} + (x\phi ){\kern 1pt} ']}}^{2}}} \right\}. \\ \end{gathered} $$
(28)

The prime in the designations of functions means the derivative with respect to the corresponding argument. Henceforward, without loss of generality, we can accept \({{b}_{0}} = 1\) and write, using (25),

$${{f}_{1}} = \left\{ { - \frac{3}{2}\mathcal{F}(\mathcal{E}) + \mathcal{F}{\kern 1pt} '(\mathcal{E}){\kern 1pt} [\mathcal{E} - {{\Phi }_{1}}(x)]} \right\}.$$
(29)

It is easy to verify that the perturbation of the system’s mass due to such perturbations of the DF and density is indeed equal to zero in both the first and second orders: \(\int {{\rho }_{{1,2}}}{{d}^{3}}{\mathbf{r}} = \int {{{d}^{3}}{\mathbf{r}}} \int {{{f}_{{1,2}}}} {{d}^{3}}{\mathbf{v}} = 0\).

First-order perturbations of \({{\Phi }_{1}}\), \({{\rho }_{1}}\), and \({{F}_{1}}\) comprise a test perturbation corresponding to the eigenfrequency \(\omega = 0\). If we analyze the perturbation dynamics by solving a system of evolution equations for the amplitudes of the Fourier harmonics \({{f}_{n}}\) of the perturbed DF, specifying the initial DF in the form (27) and the potential in the form (25), we must obtain \(\partial {{f}_{n}}{\text{/}}\partial t = 0\).

Indeed, the linearized kinetic equation for radially perturbed \(f \equiv {{F}_{1}}\) and \(\Phi \equiv {{\Phi }_{1}}\) in the action–angle variables has the form

$$\frac{{\partial f}}{{\partial t}} = - \Omega \frac{\partial }{{\partial w}}(f + \mathcal{F}{\kern 1pt} '\Phi ),$$
(30)

where \(\Omega \equiv {{\Omega }_{R}}(\mathcal{E},L)\) is the frequency corresponding to the radial action \({{J}_{R}}\), \({{\Omega }_{R}} = \partial {{H}_{0}}{\text{/}}\partial {{J}_{R}}\); \(L = {{J}_{\theta }}\; + \;{\text{|}}{{J}_{\phi }}{\kern 1pt} {\text{|}}\) is the angular momentum; and \(w \equiv {{w}_{R}}\) is the angular variable conjugate to the radial action, \(dw{\text{/}}dt = \Omega \). In harmonics, we have

$$\frac{{\partial {\kern 1pt} {{f}_{n}}}}{{\partial t}} = - i{\kern 1pt} n\Omega {\kern 1pt} ({{f}_{n}} + \mathcal{F}'{{\Phi }_{n}}),$$
(31)

where

$$\begin{gathered} {{f}_{n}}(\mathcal{E},L;t) = \oint \,dw{\kern 1pt} f(\mathcal{E},L,w,t){\kern 1pt} {{e}^{{ - inw}}}, \\ {{\Phi }_{n}}(\mathcal{E},L;t) = \oint \,dw{\kern 1pt} \Phi (\mathcal{E},L,w,t){\kern 1pt} {{e}^{{ - inw}}}, \\ \Phi (\mathcal{E},L,w,t) \equiv \Phi (r(\mathcal{E},L,w),t). \\ \end{gathered} $$

From (29), we have

$${{f}_{n}}(\mathcal{E},L) = - \mathcal{F}{\kern 1pt} '(\mathcal{E}){\kern 1pt} {{\Phi }_{n}}(\mathcal{E},L),\quad n \ne 0,$$
(32)

for nonzero harmonics. With the help of (31), we verify that \(\partial {{f}_{n}}{\text{/}}\partial t = 0\) for all \(n\), as it should be.

Thus, we have shown that the perturbed DF for the test perturbation, indeed, remains constant when solving the evolution equation, or is an eigenfunction with the eigenvalue \(\omega = 0\).

Equation (31) was actually tested on the isochrone model (7), for which it suffices to simply obtain analytical expressions relating the radial coordinate \(r\) to the action–angle variables, or, which is equivalent here, to the variables \(\mathcal{E}\) and \(L\) and the radial angular variable w.Footnote 2 Knowing the parametric relation of \(r\) with \(w\),

$$\begin{gathered} r(\mathcal{E},L,\xi ) = \sqrt {{{{\left( {\frac{{1 - p\cos \xi }}{{2\mathcal{E}}}} \right)}}^{2}} - 1} , \\ w = \xi - p\sin \xi ,\quad p = \sqrt {{{{(1 - 2\mathcal{E})}}^{2}} - 2\mathcal{E}{\kern 1pt} {{L}^{2}}} , \\ - \pi \leqslant \xi \leqslant \pi \quad {\text{and}}\quad - {\kern 1pt} \pi \leqslant w \leqslant \pi , \\ \end{gathered} $$
(33)

it is possible to numerically perform an expansion in Fourier harmonics of the radial angular variable \(w\). Setting functions (27) and (25), respectively, as the initial perturbation of the DF \({{f}_{1}}\) and potential \({{\Phi }_{1}}\) and expanding them in harmonics,

$${{f}_{1}}(\mathcal{E},L,w) = (2\pi {{)}^{{ - 1}}}\sum \,{{f}_{n}}(\mathcal{E},L;0){\kern 1pt} {{e}^{{inw}}},$$
$${{\Phi }_{1}}(\mathcal{E},L;w) = (2\pi {{)}^{{ - 1}}}\sum \,{{\Phi }_{n}}(\mathcal{E},L;0){\kern 1pt} {{e}^{{inw}}},$$

indeed, we find that \({{f}_{n}}(\mathcal{E},L;t) = {{f}_{n}}(\mathcal{E},L;0)\).

In addition, for this model, we can control the conservation of the total mass, i.e., the vanishing of the integral of the zero harmonic of the perturbed DF over the admissible region of the phase plane \({\kern 1pt} \varpi \equiv (\mathcal{E},L)\) in the model, \({{M}_{1}} = (2\pi {{)}^{2}}\int \,{{d}^{2}}\varpi {\kern 1pt} {{f}_{{n = 0}}}(\mathcal{E},L) = 0\):

$$\begin{gathered} {{M}_{1}} = (2\pi {{)}^{2}}{\kern 1pt} \int\limits_0^{1/2} \frac{{d\mathcal{E}}}{{\Omega (\mathcal{E})}}\int\limits_0^{L_{{{\text{circ}}}}^{2}(\mathcal{E})} {\kern 1pt} d({{L}^{2}})\left\{ {2\pi \left[ { - \frac{3}{2}{\kern 1pt} \mathcal{F}(\mathcal{E})} \right.} \right. \\ \left. {\left. {\mathop + \limits_{_{{}}} \;\mathcal{E}\mathcal{F}{\kern 1pt} '(\mathcal{E})} \right] - \mathcal{F}{\kern 1pt} '(\mathcal{E}){\kern 1pt} {{\Phi }_{{n = 0}}}(\mathcal{E},L)} \right\} = 0, \\ \end{gathered} $$

where \({{L}_{{{\text{circ}}}}}(\mathcal{E}) = \frac{{1 - 2\mathcal{E}}}{{\sqrt {2\mathcal{E}} }}\) is the line of circular orbits and \({{\Phi }_{{n = 0}}}(\mathcal{E},L) = 2\pi {\kern 1pt} \frac{{{{{(2\mathcal{E})}}^{{3/2}}}}}{{\sqrt {4 + {{L}^{2}}} }}{\kern 1pt} \).

3 PERTURBATION ENERGY

In this section, we want to verify on a test perturbation the correctness of the expression for energy, which, being a second-order quantity, is nevertheless constructed from only first-order perturbations. We begin with gravitational potential energy. For the models under consideration, with a single scale parameter \(b\), we have

$$\begin{gathered} W = \frac{{G{{M}^{2}}}}{b}V, \\ V = - {\kern 1pt} \frac{1}{2}\int {\kern 1pt} {{x}^{2}}{\kern 1pt} {{[\phi {\kern 1pt} '(x)]}^{2}}{\kern 1pt} dx. \\ \end{gathered} $$
(34)

Expanding \(W\) in a Taylor series in \(\epsilon = \delta b{\text{/}}{{b}_{0}}\) and setting \(G = M = {{b}_{0}} = 1\), we have \(\delta W = \epsilon {\kern 1pt} {{W}_{1}} + {{\epsilon }^{2}}{\kern 1pt} {{W}_{2}} + \mathcal{O}({{\epsilon }^{3}})\), where

$${{W}_{1}} = - V,\quad {{W}_{2}} = V.$$
(35)

We emphasize that expressions (35) provide correct expressions for the first- and the second-order perturbations of the potential energy. In particular, it is with the expression for \({{\epsilon }^{2}}{\kern 1pt} {{W}_{2}}\) that we have to compare the bilinear form \(\tilde {W} = \frac{1}{2}{\kern 1pt} G{\kern 1pt} {{\epsilon }^{2}}{\kern 1pt} \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{\rho }_{1}}({\mathbf{r}}){\kern 1pt} {{\Phi }_{1}}({\mathbf{r}})\) = \( - {{(8\pi )}^{{ - 1}}}{\kern 1pt} {{\epsilon }^{2}}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{[\nabla {{\Phi }_{1}}({\mathbf{r}})]}^{2}}\), which is usually called in the literature the disturbance potential energy. On the other hand, knowing explicit expressions (25) for \({{\Phi }_{{1,2}}}({\mathbf{r}})\), we can check the correctness of expression (35) by direct integration. In the first order, \({{W}_{1}} = - \int \,dx{\kern 1pt} {{x}^{2}}{\kern 1pt} \Phi _{0}^{'}(x){\kern 1pt} \Phi _{1}^{'}(x)\), or, after the substitution \(\Phi _{0}^{'} = \phi {\kern 1pt} '\), \(\Phi _{1}^{'} = - (x\phi ){\kern 1pt} ''\),

$${{W}_{1}} = \frac{1}{2}\int {\kern 1pt} {{x}^{2}}\phi {\kern 1pt} {{'}^{2}}{\kern 1pt} dx = - V,$$
(36)

as it should be according to (35). In the second order,

$${{W}_{2}} = - \frac{1}{2}\int {\kern 1pt} {{x}^{2}}[\Phi {{_{1}^{'}}^{2}}(x) + 2\Phi _{0}^{'}(x){\kern 1pt} \Phi _{2}^{'}(x)]{\kern 1pt} dx,$$
(37)

where, given that (see (25)) \(\Phi _{2}^{'} = \frac{1}{{2{\kern 1pt} {{x}^{2}}}}[{{x}^{3}}(x\phi {\kern 1pt} '')]{\kern 1pt} '\), we obtain, after a chain of integrations by parts,

$${{W}_{2}} = - \frac{1}{2}{\kern 1pt} \int {{x}^{2}}\phi {\kern 1pt} {{'}^{2}}{\kern 1pt} dx = V < 0.$$
(38)

Again it turned out that \({{W}_{2}} = V\), as it should be according to (35). Thus, we have seen for a test perturbation that the correct expression for the second-order potential energy \({{W}_{2}}\) is obtained only by taking into account the contribution to (37) of the second-order potential \({{\Phi }_{2}}\). This means that the bilinear form

$${{\epsilon }^{2}}{\kern 1pt} {{\tilde {W}}_{2}} = - \frac{1}{2}{{\epsilon }^{2}}\int {{x}^{2}}{\kern 1pt} \Phi {{_{1}^{'}}^{2}}(x){\kern 1pt} dx,$$
(39)

comprised only of the first-order perturbations, is not the second-order potential energy.

On the other hand, it can be shown, without going beyond the linear approximation, that, for perturbations in systems with an ergodic \(F(E)\) (\(F{\kern 1pt} '{\kern 1pt} (E) < 0\)) in the absence of external forces, there is a quadratic integral of motion:Footnote 3

$${{\tilde {E}}_{{{\text{pert}}}}} = {{\epsilon }^{2}}{\kern 1pt} ({{\tilde {T}}_{2}} + {{\tilde {W}}_{2}}),$$
(40)

where

$$\begin{gathered} {{{\tilde {T}}}_{2}} = \frac{1}{2}\int \frac{{f_{1}^{2}{\kern 1pt} {{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}}}{{ - d{{F}_{0}}{\text{/}}dE}}{\kern 1pt} , \\ {{{\tilde {W}}}_{2}} = \frac{1}{2}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{\Phi }_{1}}({\mathbf{r}}){\kern 1pt} {{\rho }_{1}}({\mathbf{r}}) = - {\kern 1pt} \frac{1}{{8\pi }}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{[\nabla {{\Phi }_{1}}({\mathbf{r}})]}^{2}}. \\ \end{gathered} $$
(41)

This expression for \({{\tilde {E}}_{{{\text{pert}}}}}\) can be obtained by considering the work performed on the system by an external force \( - \epsilon {\kern 1pt} \nabla {{\Phi }_{{{\text{ext}}}}}\), which is considered a first-order quantity (see [33], as well as [1], Section 5.4.2):

$$\frac{{d{{{\tilde {E}}}_{{{\text{pert}}}}}}}{{dt}} = - {{\epsilon }^{2}}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{1}}({\mathbf{r}},{\mathbf{v}},t){\kern 1pt} {\mathbf{v}} \cdot \nabla {{\Phi }_{{{\text{ext}}}}}{\kern 1pt} .$$
(42)

That is why the quantity \({{\tilde {E}}_{{{\text{pert}}}}}\) is usually associated with the total perturbation energy, with \({{\epsilon }^{2}}{\kern 1pt} {{\tilde {T}}_{2}}\) and \({{\epsilon }^{2}}{\kern 1pt} {{\tilde {W}}_{2}}\) being the kinetic and potential parts, respectively.Footnote 4

However, our example with the test perturbation shows that \({{\epsilon }^{2}}{\kern 1pt} {{\tilde {W}}_{2}}\) is not true potential energy, implying \({{\epsilon }^{2}}{\kern 1pt} {{\tilde {T}}_{2}}\) might not be true kinetic energy either.

Indeed, the correct expression for the second-order kinetic energy—let us call it \({{T}_{2}}\)—has the form:

$${{T}_{2}} = {{\epsilon }^{2}}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{2}}({\mathbf{r}},{\mathbf{v}})\frac{{{{v}^{2}}}}{2}{\kern 1pt} ,$$
(43)

and the correct expression for the total second-order perturbation energy, which we denote by \({{E}_{{{\text{pert}}}}}\) (without \(\widetilde {\,\,\,}\)), has the form

$${{E}_{{{\text{pert}}}}} = {{T}_{2}} + {{W}_{2}},$$
(44)

or

$$\begin{gathered} {{E}_{{{\text{pert}}}}} = {{\epsilon }^{2}}\left\{ {\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{2}}({\mathbf{r}},{\mathbf{v}}){\kern 1pt} \frac{{{{v}^{2}}}}{2}} \right.{\kern 1pt} \; \\ \left. {\mathop + \limits_{_{{_{{_{{_{{_{{}}}}}}}}}}}^{} \;\left[ {{{{\tilde {W}}}_{2}} + \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{{f}_{2}}({\mathbf{r}},{\mathbf{v}}){\kern 1pt} {{\Phi }_{0}}({\mathbf{r}})} \right]} \right\}. \\ \end{gathered} $$
(45)

Note that the expression in square brackets in formula (45) is a correct expression for potential energy \({{W}_{2}}\), taking into account the second-order contribution. This equation can be written as

$${{E}_{{{\text{pert}}}}} = {{\epsilon }^{2}}\left[ {\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{2}}({\mathbf{r}},{\mathbf{v}}){\kern 1pt} E + {{{\tilde {W}}}_{2}}} \right]{\kern 1pt} ,$$
(46)

where \(E = \frac{1}{2}{{v}^{2}} + {{\Phi }_{0}}({\mathbf{r}})\) is the energy of the star, which is an integral of the motion of the unperturbed system. We see that the correct expression (46) for the perturbation energy, which includes all contributions, differs from the formula (40), accepted in the literature (see, e.g., [1]), in different expressions for both potential and kinetic energy: \({{\tilde {T}}_{2}} \ne {{T}_{2}}\), \({{\tilde {W}}_{2}} \ne {{W}_{2}}\). Nevertheless, it turns out that the time derivatives of their sums are the same. In other words, although \({{W}_{2}}(t) \ne {{\tilde {W}}_{2}}(t)\) and \({{T}_{2}}(t) \ne {{\tilde {T}}_{2}}(t)\), the sum \({{W}_{2}} + {{T}_{2}}\) is equal to the sum \({{\tilde {W}}_{2}} + {{\tilde {T}}_{2}}\) up to an additive constant. Let us show it.

In the second order of the kinetic equation, we have

$$\begin{gathered} \left( {\frac{\partial }{{\partial t}} + {\mathbf{v}}\frac{\partial }{{\partial {\mathbf{r}}}} - \frac{{\partial {{\Phi }_{0}}}}{{\partial {\mathbf{r}}}}\frac{\partial }{{\partial {\mathbf{v}}}}} \right){\kern 1pt} {{f}_{2}} \\ = \frac{{\partial {{\Phi }_{2}}}}{{\partial {\mathbf{r}}}}\frac{{\partial {{F}_{0}}}}{{\partial {\mathbf{v}}}} + \frac{{\partial {{\Phi }_{1}}}}{{\partial {\mathbf{r}}}}\frac{{\partial {{f}_{1}}}}{{\partial {\mathbf{v}}}} + \frac{{\partial {{\Phi }_{{{\text{ext}}}}}}}{{\partial {\mathbf{r}}}}\frac{{\partial {{f}_{1}}}}{{\partial {\mathbf{v}}}}{\kern 1pt} . \\ \end{gathered} $$
(47)

Here, we have also added to the right-hand side an external potential \({{\Phi }_{{{\text{ext}}}}}({\mathbf{r}})\), which should serve as a source of change in the total energy, since the gravitational force \( - \nabla {{\Phi }_{{{\text{ext}}}}}\) does work on the stars of the system. Multiplying both sides of (47) by \(E\), integrating over the phase volume, and taking into account that \(E\) is the integral of unperturbed motion, we get

$$\begin{gathered} \frac{d}{{dt}}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{2}}({\mathbf{r}},{\mathbf{v}}){\kern 1pt} E \\ = \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} E{\kern 1pt} \left( {\frac{{\partial {{\Phi }_{2}}}}{{\partial {\mathbf{r}}}}\frac{{\partial {{F}_{0}}}}{{\partial {\mathbf{v}}}} + \frac{{\partial {{\Phi }_{1}}}}{{\partial {\mathbf{r}}}}\frac{{\partial {{f}_{1}}}}{{\partial {\mathbf{v}}}} + \frac{{\partial {{\Phi }_{{{\text{ext}}}}}}}{{\partial {\mathbf{r}}}}\frac{{\partial {{f}_{1}}}}{{\partial {\mathbf{v}}}}} \right). \\ \end{gathered} $$
(48)

The first term on the right-hand side of (48) vanishes due to the antisymmetry of the integrand with respect to \({\mathbf{v}}\), since \(E{\kern 1pt} \partial {{F}_{0}}{\text{/}}\partial {\mathbf{v}} = E{\kern 1pt} {\mathbf{v}}{\kern 1pt} F_{0}^{'}(E)\). The second term, after a chain of transformations, turns into \( - \frac{{d{{{\tilde {W}}}_{2}}}}{{dt}}\). Indeed, for it we have from (48):

$$\begin{gathered} \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} \frac{{\partial {{\Phi }_{1}}}}{{\partial {\mathbf{r}}}}\int \,{{d}^{3}}{\mathbf{v}}{\kern 1pt} E{\kern 1pt} \frac{{\partial {{f}_{1}}}}{{\partial {\mathbf{v}}}} = - \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} \frac{{\partial {{\Phi }_{1}}}}{{\partial {\mathbf{r}}}}\int \,{{d}^{3}}{\mathbf{v}}{\kern 1pt} ({\mathbf{v}}{\kern 1pt} {{f}_{1}}) \\ = \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{\Phi }_{1}}({\mathbf{r}})\frac{\partial }{{\partial {\mathbf{r}}}}\int \,{{d}^{3}}{\mathbf{v}}{\kern 1pt} ({\mathbf{v}}{\kern 1pt} {{f}_{1}}) = - \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{\Phi }_{1}}\frac{{\partial {{\rho }_{1}}}}{{\partial t}} \\ = - \int \,{{d}^{3}}{\mathbf{r}}\frac{{\partial {{\Phi }_{1}}}}{{\partial t}}{{\rho }_{1}}({\mathbf{r}}) = - \frac{1}{2}\int \,{{d}^{3}}{\mathbf{r}}\left[ {{{\Phi }_{1}}\frac{{\partial {{\rho }_{1}}}}{{\partial t}} + \frac{{\partial {{\Phi }_{1}}}}{{\partial t}}{{\rho }_{1}}} \right] \\ = - \frac{d}{{dt}}\left[ {\frac{1}{2}\int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} ({{\rho }_{1}}{\kern 1pt} {{\Phi }_{1}})} \right] = - \frac{{d{{{\tilde {W}}}_{2}}}}{{dt}}{\kern 1pt} . \\ \end{gathered} $$
(49)

The third term turns into \( - \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{1}}{\kern 1pt} {\mathbf{v}}{\kern 1pt} \nabla {{\Phi }_{{{\text{ext}}}}}\):

$$\begin{gathered} \int \,{{d}^{3}}{\mathbf{r}}\frac{{\partial {{\Phi }_{{{\text{ext}}}}}}}{{\partial {\mathbf{r}}}}\int \,{{d}^{3}}{\mathbf{v}}{\kern 1pt} E\frac{{\partial {{f}_{1}}}}{{\partial {\mathbf{v}}}} \\ = - \int \,{{d}^{3}}{\mathbf{r}}\frac{{\partial {{\Phi }_{{{\text{ext}}}}}}}{{\partial {\mathbf{r}}}}\int \,{{d}^{3}}{\mathbf{v}}{\kern 1pt} {\mathbf{v}}{\kern 1pt} {{f}_{1}} = - \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{1}}{\kern 1pt} {\mathbf{v}}{\kern 1pt} \nabla {{\Phi }_{{{\text{ext}}}}}. \\ \end{gathered} $$
(50)

As a result, combining (48)–(50), we find that the rate of change in the system’s energy, caused by the work of the external force, is

$$\begin{gathered} \frac{{d{{E}_{{{\text{pert}}}}}}}{{dt}} \equiv {{\epsilon }^{2}}\frac{d}{{dt}}\left[ {{\kern 1pt} \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{2}}\frac{{{{v}^{2}}}}{2} + {{W}_{2}}} \right] \\ = - {{\epsilon }^{2}}{\kern 1pt} \int \,{{d}^{3}}{\mathbf{r}}{\kern 1pt} {{d}^{3}}{\mathbf{v}}{\kern 1pt} {{f}_{1}}{\kern 1pt} {\mathbf{v}}{\kern 1pt} \nabla {{\Phi }_{{{\text{ext}}}}}. \\ \end{gathered} $$
(51)

Comparing the right-hand sides of (42) and (51), we find \(d{{E}_{{{\text{pert}}}}}{\text{/}}dt = d{{\tilde {E}}_{{{\text{pert}}}}}{\text{/}}dt\). Thus, the true total perturbation energy \({{E}_{{{\text{pert}}}}}(t)\) differs from the bilinear construction \({{\tilde {E}}_{{{\text{pert}}}}}(t)\), usually called the perturbation energy, by a constant. This means that both of these second-order quantities, in the absence of external forces, are conserved during the evolution of the system. Therefore, the construction \({{\tilde {E}}_{{{\text{pert}}}}}(t)\) expressed by relations (40) and (41), by analogy with the linear theory of shear fluid flows, can be called pseudoenergy. Recall that, in the theory of shear flows, there is also the concept of a pseudoenergy integral, which is constructed as a bilinear form of first-order perturbations. Pseudoenergy differs from true energy, which must be calculated taking into account second-order perturbations (see [22, 34]).

Note that pseudoenergy and true second-order energy may differ in sign. In particular, it can be shown that the pseudoenergy of the eigenmodes of systems with a decreasing ergodic DF, \(F_{0}^{'}(E) < 0\) (i.e., the van Kampen modes [12]), is positive, although the true energy can be of any sign. This circumstance may be important from the viewpoints of attempts to construct the thermodynamics of star clusters based on the involvement of van Kampen waves [12]. For the success of such attempts, the positive sign of the energy is critical. However, for example, for our test perturbation, the true energy is negative. This follows from the virial relation \(2T + W = 0\), which must be satisfied in all orders of perturbation theory due to the stationarity of the perturbation. In particular, in the second order, we have \(2{{T}_{2}} + {{W}_{2}} = 0\). Hence \({{E}_{{{\text{pert}}}}} = {{\epsilon }^{2}}({{T}_{2}} + {{W}_{2}}) = \frac{1}{2}{{\epsilon }^{2}}{{W}_{2}}\). But, as follows from (38), \({{W}_{2}} < 0\).

4 CONCLUSIONS

We have shown that, for testing codes in a numerical study of the dynamics of small perturbations in spherical stellar systems, there is a control radial stationary perturbation for which the distribution function, density, and gravitational potential can be expressed in explicit form. This perturbation is relevant for ergodic systems, in the models of which there is a single scale factor of the dimension of length. Examples of several models of this type known in the literature have been given. When solving the eigenvalue problem (by any of the known matrix methods), this perturbation must give the eigenfrequency \(\omega = 0\) and the corresponding known eigenfunctions, and, when solving the initial value problem for the perturbation of the DF \(f({\mathbf{r}},{\mathbf{v}};t)\), must confirm the conservation of the DF at each point of the phase space, \(f({\mathbf{r}},{\mathbf{v}};t) = f({\mathbf{r}},{\mathbf{v}};0)\), if the test value is taken as the initial perturbation.

In addition, in this work, we have analyzed the concept of perturbation energy, which appears in the linear theory of perturbations of collisionless stellar systems. It is known that, although the true perturbation energy, being a quantity of the second order in the perturbation amplitude, in principle, cannot be calculated within linear theory, it is possible to construct a bilinear form from first-order quantities, which is an integral of the linearized equations. This quantity is very similar to energy and comprises the sum of two contributions, which are usually called the kinetic and potential perturbation energies. For a test perturbation as an example, for which we can obtain expressions in any order of perturbation theory, we have found that the expressions for the kinetic and potential energies obtained within linear theory differ from the correct expressions for the kinetic and potential energies obtained taking into account second-order perturbations. In this work, it has been shown that, although the integral of motion represented by the correct expression for the perturbation energy and the integral corresponding to the energy constructed within linear theory (pseudoenergy) do not coincide, they differ only by a time-independent value (constant). However, these quantities may differ in sign, which may be important for problems related to the application of van Kampen modes to stellar systems.