Radial orbit instability in systems of highly eccentric orbits: Antonov problem reviewed

TL;DR

This paper refines the proof of radial orbit instability in spherical stellar systems with nearly radial orbits, showing that instability grows rapidly and invalidates orbital approximation in highly anisotropic systems.

Contribution

It introduces non-singular models to analyze instability, linking previous proofs with new integral equations for systems close to purely radial orbits.

Findings

Maximum growth rates tend to infinity as radial anisotropy increases.

Instability occurs for both even and odd spherical harmonics.

Unstable modes develop rapidly with eigenfrequencies comparable to orbital frequencies.

Abstract

Stationary stellar systems with radially elongated orbits are subject to radial orbit instability -- an important phenomenon that structures galaxies. Antonov (1973) presented a formal proof of the instability for spherical systems in the limit of purely radial orbits. However, such spheres have highly inhomogeneous density distributions with singularity $\sim 1/r^2$, resulting in an inconsistency in the proof. The proof can be refined, if one considers an orbital distribution close to purely radial, but not entirely radial, which allows to avoid the central singularity. For this purpose we employ non-singular analogs of generalised polytropes elaborated recently in our work in order to derive and solve new integral equations adopted for calculation of unstable eigenmodes in systems with nearly radial orbits. In addition, we establish a link between our and Antonov's approaches and uncover the meaning of infinite entities in the purely radial case. Maximum growth rates tend to infinity as the system becomes more and more radially anisotropic. The instability takes place both for even and odd spherical harmonics, with all unstable modes developing rapidly, i.e. having eigenfrequencies comparable to or greater than typical orbital frequencies. This invalidates orbital approximation in the case of systems with all orbits very close to purely radial.

Equations

$$F(E,L)=\frac{H(L_T-L)}{L_T^2}{F}_0(E),$$

$$F_0(E)=\frac{N(L_T,q)}{4{\pi }^3}(-2E{)}^q.$$

$$F(E,L)=\delta (L^2)F_0(E)=\frac{\pi }{r^2}\delta (v_{\theta })\delta (v_{\phi })F_0(E),$$

$$\begin{gathered} f_1(t,r,\theta ,\phi ,v_r,v_{\theta },v_{\phi })=A(t,r,\theta ,\phi ,v_r)\delta (v_{\theta })\delta (v_{\phi }) \\ +B(t,r,\theta ,\phi ,v_r){\delta }^{\prime }(v_{\theta })\delta (v_{\phi })+C(t,r,\theta ,\phi ,v_r)\delta (v_{\theta }){\delta }^{\prime }(v_{\phi }), \end{gathered}$$

$$\begin{gathered} \frac{\partial f_1}{\partial t}+v_r\frac{\partial f_1}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial f_1}{\partial \theta }+\frac{v_{\phi }}{r\sin \theta }\frac{\partial f_1}{\partial \phi } \\ +\left(\frac{v_{\theta }^2+v_{\phi }^2}{r}-\frac{d{\Phi }_0}{dr}\right)\frac{\partial f_1}{\partial v_r}-\left(\frac{v_r{v}_{\theta }}{r}-\cot \theta \frac{v_{\phi }^2}{r}\right)\frac{\partial f_1}{\partial v_{\theta }} \\ -\left(\frac{v_r{v}_{\phi }}{r}+\cot \theta \frac{v_{\theta }{v}_{\phi }}{r}\right)\frac{\partial f_1}{\partial v_{\phi }}=\frac{\partial {\Phi }_1}{\partial r}\frac{\partial F}{\partial v_r} \\ +\frac{1}{r}\frac{\partial {\Phi }_1}{\partial \theta }\frac{\partial F}{\partial f{v}_{\theta }}+\frac{1}{r\sin \theta }\frac{\partial {\Phi }_1}{\partial \phi }\frac{\partial F}{\partial v_{\phi }}, \end{gathered}$$

$$\frac{\partial A}{\partial t}+\hat{D}A+\frac{2v_r}{r}A-\frac{1}{r}\frac{\partial B}{\partial \theta }-\frac{\cot \theta }{r}B$$

$$-\frac{1}{r\sin \theta }\frac{\partial C}{\partial \phi }=\frac{\pi }{r^2}{F}_0^{\prime }(E)v_r\frac{\partial {\Phi }_1}{\partial r},$$

$$\frac{\partial B}{\partial t}+\hat{D}B+\frac{3v_r}{r}B=\frac{\pi }{r^3}{F}_0(E)\frac{\partial {\Phi }_1}{\partial \theta },$$

$$\frac{\partial C}{\partial t}+\hat{D}C+\frac{3v_r}{r}C=\frac{\pi }{r^3\sin \theta }{F}_0(E)\frac{\partial {\Phi }_1}{\partial \phi }$$

$$\hat{D}=v_r\frac{\partial }{\partial r}-\frac{d{\Phi }_0}{dr}\frac{\partial }{\partial v_r},$$

$${\Phi }_1\equiv \chi (t,r)P_l(\cos \theta ),$$

$$A(t,r,\theta ,\phi ,v_r)=\frac{\mathcal{A}(t,r,v_r)}{r^2}{P}_l(\cos \theta ),$$

$$B(t,r,\theta ,\phi ,v_r)=\frac{\mathcal{B}(t,r,v_r)}{r^3}\frac{d{P}_l(\cos \theta )}{d\theta }.$$

$${\rho }_1(t,r,\theta ,\phi )\equiv \Pi (t,r)P_l(\cos \theta )$$

$$\Pi (t,r)=\frac{1}{r^2}\int \mathcal{A}(t,r,v_r)d{v}_r.$$

$$\frac{\partial \mathcal{A}}{\partial t}+\hat{D}\mathcal{A}+\frac{l(l+1)}{r^2}\mathcal{B}=\pi v_r\frac{d\chi }{dr}{F}_0^{\prime }(E),$$

$$\frac{\partial \mathcal{B}}{\partial t}+\hat{D}\mathcal{B}=\pi \chi (r)F_0(E),$$

$$\chi (r)=-\frac{4\pi G}{2l+1}\int d{r}^{\prime }\int d{v}_r^{\prime }\mathcal{A}(r^{\prime },v_r^{\prime }){\mathcal{F}}_l(r,r^{\prime }),$$

$$w=\Omega (E)\underset{0}{\overset{r}{\int }}\frac{d{r}^{\prime }}{v_r(E,r^{\prime })},$$

$$v_r=\pm \sqrt{2[E-{\Phi }_0(r)]}.$$

$$\{\chi(t,w,E),{\cal A}(t,w,E),{\cal B}(t,w,E)\}=\\ =\sum\limits_{n=-\infty}^{\infty}\{\Phi_{n/2}(t,E),A_{n/2}(t,E),B_{n/2}(t,E)\}\,e^{inw/2}\$$

$$\theta \to \pi -\theta ,\phi \to \pi +\phi .$$

$$\frac{\partial \mathcal{A}}{\partial t}+{\Omega }_1\frac{\partial \mathcal{A}}{\partial w}+\frac{l(l+1)}{r^2(E,w)}\mathcal{B}=\pi {\Omega }_1\frac{\partial \chi (t,E,w)}{\partial w}{F}_0^{\prime }(E),$$

$$\frac{\partial \mathcal{B}}{\partial t}+{\Omega }_1\frac{\partial \mathcal{B}}{\partial w}=\pi \chi (t,E,w)F_0(E),$$

$$\chi (t,E,w)=-\frac{4\pi G}{2l+1}\int \frac{d{E}^{\prime }}{{\Omega }_1(E^{\prime })}$$

$$\times \int d{w}^{\prime }\mathcal{A}(t,E^{\prime },w^{\prime }){\mathcal{F}}_l[r(E,w),r^{\prime }(E^{\prime },w^{\prime })].$$

$$r(2\pi -w)=r(-w)=r(w),$$

$$\begin{gathered} \{{\Phi }_{n/2}(t,E),A_{n/2}(t,E),B_{n/2}(t,E)\} \\ =\frac{1}{4\pi }\underset{-2\pi }{\overset{2\pi }{\int }}\{\chi (t,r),\mathcal{A}(t,r,v_r),\mathcal{B}(t,r,v_r)\}{e}^{-inw/2}dw \end{gathered}$$

$$\begin{gathered} \{{\Phi }_k(t,E),A_k(t,E),B_k(t,E)\} \\ =\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}\{\chi (t,r),\mathcal{A}(r,v_r),\mathcal{B}(r,v_r)\}\cos (kw)dw, \end{gathered}$$

$$-i(\omega -k\Omega )A_k+l(l+1)\sum _{k^{\prime }=-\infty }^{\infty }{p}_{k-k^{\prime }}{B}_{k^{\prime }}=i\pi k{\Phi }_k{F}_0^{\prime }(E),$$

$$-i(\omega -k\Omega )B_k=\pi {\Phi }_k{F}_0(E)$$

$${\Phi }_k(E)=-\frac{2G}{2l+1}\int \frac{d{E}^{\prime }}{\Omega (E^{\prime })}\sum _{k^{\prime }=-\infty }^{\infty }{\mathcal{K}}_{k{k}^{\prime }}^{\mathrm{even}}(E,E^{\prime })A_{k^{\prime }}(E^{\prime }),$$

$$p_k(E)=\frac{1}{2\pi }\oint \frac{dw}{r^2(E,w)}{e}^{-ikw}=\frac{1}{2\pi }\oint \frac{\cos (kw)}{r^2(E,w)}dw$$

$${\mathcal{K}}_{k{k}^{\prime }}^{\mathrm{even}}(E,E^{\prime })=4\underset{0}{\overset{\pi }{\int }}dw\underset{0}{\overset{\pi }{\int }}d{w}^{\prime }\cos (kw)\cos (k^{\prime }{w}^{\prime }){\mathcal{F}}_l(r,r^{\prime }).$$

$$p_{0}=\Big{\langle}\frac{1}{r^{2}}\Big{\rangle}\equiv\frac{1}{2\pi}\oint\frac{dw}{r^{2}}=\frac{\Omega}{\pi}\int\limits_{0}^{r_{\rm max}(E)}\frac{dr}{r^{2}\,\sqrt{2E-2\Phi_{0}(r)}}\ ,$$

Full text

1 Introduction

The radial orbit instability (ROI), first mentioned in a preprint by Polyachenko & Shukhman (1972), plays an important role in the evolution of initially spherically symmetric and axisymmetric systems leading to bar-like perturbations. It has been widely studied both analytically (Antonov, 1973; Polyachenko & Shukhman, 1981; Palmer & Papaloizou, 1987; Weinberg, 1991; Saha, 1991; Palmer, 1994; Polyachenko et al., 2011, 2015) and numerically (Polyachenko, 1981; Merritt, 1985, 1987; Barnes et al., 1986; Aguilar & Merritt, 1990; Bertin et al., 1994; Meza & Zamorano, 1997; Trenti & Bertin, 2006). There are two basic candidates for a physical mechanism of ROI: an analog of Jeans instability in the anisotropic media and an orbital approach based on tendency of any pair of orbits to align under their mutual gravity. Discussion on these topics can be found in Polyachenko & Shukhman (2015).

A distinct approach is suggested by Maréchal & Perez (2010), who give an example of dissipation-induced ROI. A comprehensive modern review on ROI can be found in Maréchal & Perez (2012), who also suggest a new symplectic method for exploring stability of equilibrium gravitating systems.

Antonov (1973) presents a first formal proof of ROI for purely radial motion using the Lyapunov method. However, his proof is doubtful: the Lyapunov function is ill-defined due to a divergence of its time derivative at the lower limit of integration. Although the main conclusion of the paper is correct, a rigorous examination of the purely radial case is still needed. The goal of this paper is to reconsider the Antonov problem by applying our technique of an eigenvalue problem in the form of integral equations.

For this purpose, we shall use a general family of models

[TABLE]

where $H(x)$ is the Heaviside function. We retain an arbitrary form for $F_{0}(E)$ whenever possible, otherwise we admit a polytropic law

[TABLE]

Here $E=\frac{1}{2}(v_{r}^{2}+v_{\perp}^{2})+\Phi_{0}(r)\leq 0$ and $L=rv_{\perp}$ are the energy and absolute value of the angular momentum, respectively; $\Phi_{0}(r)$ is the unperturbed gravitational potential. The additive constant in $\Phi_{0}$ is chosen so that the potential vanishes at the outer radius of sphere $R$; the normalization constant $N(L_{T},q)$ is chosen so that the total mass of the system is $M$. In the calculations below we shall assume that $M=R=G=1$. Equilibrium properties of family (1.1) with polytropic dependence from energy (1.2), called softened polytrope models, are specially built to consider the limit of purely radial motion and studied in our paper (Polyachenko et al., 2013). Stability properties of some series (fixed $q$) are studied in Polyachenko & Shukhman (2015).

The polytropic law includes a series of mono-energetic models, in which all stars have zero total energy, at the limit $q\to-1$ (e.g., Gelfand & Shilov, 1959). A limit $L_{T}\to 0$ in this series gives a well-known Agekyan (1962) model which was employed in the Antonov’s work. The model is particularly useful in our case, since it provides the simplest eigenvalue equations, yet preserving all features of interest.

As is already said, our proof is based on analysis and solution of characteristic equations for eigenmodes – spherical harmonics and corresponding complex frequencies $\omega$, such that ones with the positive imaginary parts give unstable solutions. In Section 2 we derive the equation for a model with purely radial orbits, using delta-function expansion technique (Fridman & Polyachenko, 1984). The unperturbed distribution function (DF) of the purely radial system is proportional to the Dirac delta-function of the angular momentum (2.1), while the perturbed DF is a linear combination of the delta-function and its derivatives (2.2). The linearised kinetic equation and Poisson equation provide matrix equations (2.36) and (2.48), for even and odd spherical harmonics, respectively. Both of them contain infinite entities $p_{k}$ defined by (2.29) which are a manifestation of the central singularity.

In Section 3 we use the integral equation technique for the two-parametric family of models (1.1) with nearly radial orbits (Polyachenko et al., 2007; Polyachenko & Shukhman, 2015). Since this family includes the purely radial model of Section 2, we can get a link between different parts of the integral equations obtained in Sections 2 and 3, as the control parameter $L_{T}$ in the DF approaches zero. In particular, we infer the meaning of the infinite entities, eqs. (3.11, 3.12).

Then, this finding helps us in Section 4 to reduce further the legitimate integral equations of Section 3 for nearly radial orbits to fairly compact limiting integral equations ($L_{T}\ll 1$), for even and odd spherical harmonics, (4.12) and (4.13), respectively. They allow to prove existence of the aperiodic even unstable spherical solutions, and absence of the odd unstable spherical solutions. The analytical results are accompanied in Section 5 with numerical eigenmodes’ calculation for series $q=-1$ and $q=-1/2$.

In Section 6, we show how the orbital approach breaks down in spherical systems with orbits very close to radial. Comparison of our numerical results with qualitative results by (1991) shows that orbital approach is satisfactory for the systems with orbits of moderate eccentricity only.

Lastly, Section 7 contains a summary and conclusion. Appendix A is devoted to Antonov’s ‘proof’ of the existence of ROI (in our terms and notations) with the help of Lyapunov function, as a reminder and demonstration of difficulties appearing in the investigation of the systems with pure radial systems. Appendix B clarifies the sense of diverging coefficients which appear in the equations for purely radial models with the help of limiting procedure from models with finite dispersion over the angular momentum ($L_{T}\neq 0$).

2 Pure radial motion: $\delta$-function expansion

Purely radial orbits possess zero angular momentum, $L=0$. Thus, for systems with purely radial motion, we demand

[TABLE]

where $\delta(x)$ is the Dirac delta-function. The analysis for instability prescribes the following ansatz for the perturbed DF:

[TABLE]

where $\delta^{\prime}$ denotes a derivative of the delta-function. From the linearized kinetic equation

[TABLE]

relations between decomposition coefficients $A$, $B$, and $C$ can be obtained:

[TABLE]

Here ${\hat{D}}$ is a differential operator,

[TABLE]

and $\Phi_{1}$ denotes a perturbed potential depending on the polar angle $\theta$ only through Legendre polynomials $P_{l}$,

[TABLE]

since the eigenmode spectrum does not depend on the azimuthal number $m$ (e.g., Fridman & Polyachenko, 1984; Bertin et al., 1994).

Substitution to (2.4–2.6) gives $C=0$, while $A\propto P_{l}(\cos\theta)$, and $B\propto{\rm d}P_{l}(\cos\theta)/{\rm d}\theta$. It is convenient to introduce new functions ${\cal A}$ and ${\cal B}$ independent of the angles:

[TABLE]

The perturbed density

[TABLE]

is an integral from ${\cal A}$ over the radial velocity,

[TABLE]

The eqs. (2.4–2.6) and the Poisson equation for the new functions take the form:

[TABLE]

with ${\cal F}_{l}(r,r^{\prime})={r_{<}^{l}}/{r_{>}^{l+1}}$, $r_{<}={\rm min}(r,r^{\prime})$ and $r_{>}={\rm max}(r,r^{\prime})$.

Below, we shall explore the system of eqs. (2.13–2.15) in terms of action–angle variables. Radial orbits can be treated as highly eccentric ellipses with vanishingly small minor axis, see Fig. 1. The stellar position is fixed by four variables, three of which determine the orbit length and orientation (e.g., energy $E$, angles $\theta$ and $\varphi$), and the last one – radial angle variable $w$ – sets the position along the orbit,

[TABLE]

where $\Omega(E)$ is the frequency of radial oscillations, $v_{r}$ is the radial velocity:

[TABLE]

During full revolution, star’s angular variable changes in the range $-2\pi\leq w\leq 2\pi$. Therefore, the most general functions of $w$ have a period of $4\pi$, and their Fourier expansions should read

[TABLE]

with $n$ running over all integers. Thus, variables $(t,r,v_{r})$ are changed to $(t,E,w)$.

As the star travels from the upper part of the orbit $0<w<2\pi$ to the lower part $-2\pi<w<0$, the polar and azimuthal angles change discontinuously:

[TABLE]

This results in additional factor in case of odd spherical harmonics, so further analysis should be done separately for even and odd cases.

2.1 Equations for even spherical harmonics

In action–angle variables, the evolutionary equations for perturbations in purely radial systems are

[TABLE]

Radius of a star as a function of angle $w$ obeys the following symmetry conditions:

[TABLE]

thus coefficients

[TABLE]

vanish for odd $n$, and equal to

[TABLE]

otherwise ($k=n/2$). In this case, periodicity changes to $2\pi$ due to symmetry of potential and density perturbations with respect to transformation $w\to-w$.

Assuming that perturbations are $\propto\exp(-i\omega t)$, one can obtain from (2.20)–(2.22)

[TABLE]

and

[TABLE]

where $k$ and $k^{\prime}$ are integers,

[TABLE]

and

[TABLE]

Eq. (2.29) emphasises an issue arising in systems with purely radial orbits – the integrals diverge in the centre ($w\to 0,|2\pi|$). Thus, these expressions for $p_{k}$ require an interpretation. Note that a similar difficulty appeared in Antonov (1973), but then no adequate attention has been paid. For example, $\varpi(E)\equiv p_{0}(E)$ is $1/r^{2}$, averaged along the orbit:

[TABLE]

and diverges evidently at $r=0$, since singularity of $\Phi_{0}$ is weaker than $1/r^{2}$. We plan to tackle the issue employing a family of models with nearly radial orbits and study the system of interest by considering more and more radially anisotropic systems (see Section 3).

With (2.26) and (2.27), one can exclude $A_{k}$ and $B_{k}$ from the equation in favour of $\Phi_{k}$:

[TABLE]

For some $F_{0}(E)$ (e.g., polytropes (1.2) with $q<0$) the integral from the term including $dF_{0}/dE^{\prime}$ diverges. In this case one should use the Lagrangian form, which is obtained formally by integration by parts and omission of the surface term:

[TABLE]

(see Polyachenko & Shukhman (2015) for details). For further analysis, it is convenient to use an alternative form of the last equation. Using the identity provided $m\neq n^{\prime}$

[TABLE]

one can have

[TABLE]

In the final equation, we have separated the last term which retains even in the case of radial oscillations $l=0$.

To find out the meaning of integrals (2.29) in expressions for $p_{k}$, a specific model is not important, since $p_{k}(E)$ depends on the orbit, but not on the orbit distribution over the phase space. For simplicity we consider a monoenergetic model corresponding to $q\to-1$ limit in (1.2), which leads to algebraic equations:

[TABLE]

where $\Phi_{k}$ now denotes $\Phi_{k}(E=0)$.

2.2 Equations for odd spherical harmonics

The jump of the polar angle (2.19) gives rise to additional factor

[TABLE]

in case of the odd spherical harmonics, i.e.

[TABLE]

The functions to be expanded

[TABLE]

are antisymmetric, i.e.

[TABLE]

Now the expansion coefficients for even $n$ vanish, while for odd $n$ one has ($n=2k+1$):

[TABLE]

From the eqs. similar to (2.20) and (2.21), it follows that expansion coefficients $B_{n/2}(t,E)$ also vanish for even $n$. For a new set of variables $\{\bar{\Phi}_{k},\bar{A}_{k},\bar{B}_{k}\}\equiv\{\Phi_{k+1/2},A_{k+1/2},B_{k+1/2}\}$, one obtains equations for the odd spherical harmonics $l$:

[TABLE]

and

[TABLE]

where $k$ and $k^{\prime}$ are integers; $p_{k}(E)$ is given by (2.29), i.e. the same as for the even $l$;

[TABLE]

Eliminating $\bar{A}_{k}(E)$ and $\bar{B}_{k}(E)$ in favour of $\bar{\Phi}_{k}(E)$, one obtains the equations similar to (2.36) for even $l$:

[TABLE]

3 Nearly radial orbits: integral equations

The integral equations (2.36) and (2.48) are of no use, since they contain infinite coefficients $p_{k}$. In this section we consider nearly radial series of ‘dispersed’ Agekyan models ($q=-1$) with a control parameter $L_{T}$ which includes the purely radial (Agekyan) model of the previous section as a limiting case $L_{T}\to 0$. We shall see, that the dispersed models allow for a well-defined integral equations, and their solutions indeed give infinitely large growth rates in the limit of the purely radial case.

The integral equations for the nearly radial models in the Lagrangian form are (Polyachenko & Shukhman, 2015):

[TABLE]

where $\Omega_{l_{1}^{\prime}l_{2}^{\prime}}(E^{\prime},L^{\prime})\equiv l_{1}^{\prime}\,\Omega_{1}(E^{\prime},L^{\prime})+l_{2}^{\prime}\,\Omega_{2}(E^{\prime},L^{\prime})$; differentiation operator $\partial/\partial E^{\prime}$ acts both on $\Omega_{l_{1}^{\prime}l_{2}^{\prime}}(E^{\prime},L^{\prime})$ and the last row; coefficients $D_{l}^{k}$ vanish for odd $|l-k|$ and

[TABLE]

otherwise. For the given models, the right-hand side can be written as the sum of two terms:

[TABLE]

Due to orbit symmetry, $r(w)=r(-w)$, the kernel functions $\Pi_{l_{1},\,l_{2};\,l_{1}^{\prime},\,l_{2}^{\prime}}$ and unknown expansion coefficients of the potential $\phi_{l_{1}\,l_{2}}$ can be expressed in integral forms with integration reduced from $[-\pi,\pi]$ to $[0,\pi]$:

[TABLE]

[TABLE]

The angle $\Theta_{l_{1}l_{2}}(E,L,w)$ is

[TABLE]

with

[TABLE]

denoting the azimuthal change of the particle coordinate as it passes from the pericentre to the current radius $r$; $\Psi(r)\equiv-\Phi_{0}(r)$. At the apocentre, $\delta\varphi(E,L;\pi)=(\Omega_{2}/\Omega_{1})\,\pi$.

Further, we shall expand the functions entering eq. (3.1) considering $L_{T}$ as a small parameter. For nearly radial orbits, the precession velocity

[TABLE]

is small compared to frequencies $\Omega_{1,2}$. So, it can be separated out in the linear combination

[TABLE]

The angle $\Theta_{l_{1}\,l_{2}}(E,L;w)$ can be written as a sum

[TABLE]

where the expression in the square brackets retains in the limit $L\to 0$, while

[TABLE]

vanishes.

In Appendix B, we give details of the expansion of eq. (3.1) in $L_{T}$ and $\beta$ for the even spherical harmonics $l$. It should be compared with eq. (2.36) for the systems with purely radial orbits. The equations coincide entirely, if $p_{0}(E)$ is substituted by the limiting ratio $\Omega_{\rm pr}(L_{T})/L_{T}$, i.e.

[TABLE]

and coefficients $p_{k}(E)$ for $k\neq 0$ are understood as the limits

[TABLE]

Given that $(L_{T}/\pi)\int_{0}^{\pi}dw/{r^{2}}=\Omega_{2}$, one can have

[TABLE]

where the right-hand side converges in the usual sense. Note that expansion of (3.1), not given here, and comparison with (2.48) for the odd harmonics lead to the same results (3.11) and (3.12).

The obtained relation between $p_{k}$ and the limiting value of $\varpi\equiv\Omega_{\textrm{pr}}(L_{T})/L_{T}$ implies that $p_{k}$ are infinitely large. Indeed, in purely radial systems the density is necessarily singular, at least not weaker than $1/r^{2}$ (Bouvier & Janin, 1968; , 1984). Thus the potential and gravitational force are also singular and linear law of the precession rate $\Omega_{\textrm{p}r}(L)$ is no longer valid. In particular, for the softened polytropes all purely radial models ($q\leq 1/2$) give $\varpi\to\infty$ in the limit $L_{T}\to 0$ (see Fig.9b in Polyachenko et al., 2013). Besides, for dispersed Agekyan model, we found numerically that $\Omega_{\rm pr}(L_{T})\approx 0.316\,(L_{T})^{0.26}$, i.e., $\varpi\approx 0.316/(L_{T})^{0.74}$ (see below Sect. 5.1).

4 Limiting integral equations

Eqs. (3.11) and (3.13) for $p_{k}(E)$ show that infinitely large coefficients occur in (2.26)–(2.28) and (2.44)–(2.46) as $L_{T}$ goes to zero. This enables us to obtain a simplified counterparts of the stability equations. We shall start from the equations in the form (2.20)–(2.22), and assume everywhere that $L_{T}\ll 1$. The right hand side in (2.20) should be omitted since it does not contain $p_{0}$. Then, one should neglect the difference between $p_{k}$ and $p_{0}=\varpi\equiv\Omega_{\textrm{pr}}(L_{T})/L_{T}\gg 1$ since $p_{k}/p_{0}-1={\cal O}(1/p_{0})\ll 1$. The expansion $1/r^{2}=\sum p_{k}\,\exp(ikw)$ then turns into

[TABLE]

so that (2.20) and (2.21) turn into

[TABLE]

Now changing $\partial/\partial t$ to $-i\omega$ and solving the equations, taking into account symmetry of functions ${\cal R}(w)$ and $\chi(w)$, one obtains for ${\cal A}$:

[TABLE]

or

[TABLE]

where $\nu\equiv\omega/\Omega$ is the dimensionless frequency. For ${\cal B}$ the solution is

[TABLE]

In particular, ${\cal B}(E,0)$ is

[TABLE]

where

[TABLE]

Multiplying the Poisson equation (2.22) by \cos\bigl{[}\nu\,(w-\pi)\bigr{]} and integrating from [math] to $\pi$, one obtains an integral equation for even $l$:

[TABLE]

where $\Omega^{\prime}$ and $\nu^{\prime}$ denote $\Omega(E^{\prime})$ and $\omega/\Omega^{\prime}$, and the kernel is

[TABLE]

The analogous equation of odd $l$ has the form:

[TABLE]

We shall refer further to eqs. (4.9) and (4.11) as limiting integral equations. Note that they lack the advantage of the linear eigenvalue problem, since frequency $\omega$ enters into the kernel function and into the argument of sine and cosine in denominators. However, they retain their forms during the change $\omega\to-\omega$, so both equations should depend on $\omega^{2}$.

We need to emphasise that $\varpi$ depends on $L_{T}$, and that it is assumed that $L_{T}\ll 1$ and $\varpi(E,L_{T})\gg 1$. This is the only variable dependent on $L_{T}$, in all other places $L_{T}\to 0$ limit leads to finite quantities, so there we assume $L_{T}=0$.

Relative simplicity of the limiting integral equations (4.9) and (4.11) allows us to demonstrate analytically existence of aperiodic unstable solutions ($\omega=i\gamma$ with $\gamma>0$) for even $l$ and their absence for odd $l$. Introducing $\sigma(E)\equiv\gamma/\Omega(E)$ one obtains equivalent equations for $\gamma$:

for even $l$,

[TABLE]

for odd $l$,

[TABLE]

In both equations

[TABLE]

Redefinition of the eigenfunction

[TABLE]

allows one to symmetrize the integral equations. Using an integral representation for ${\cal F}_{l}(r,r^{\prime})$ through Bessel functions (e.g., Gradshteyn & Ryzhik, 2015, eq. 6.574),

[TABLE]

it can be proven that the kernel of the symmetrized equation for even $l$, ${\cal Q}_{\gamma}^{\rm even}$ is positive. So the eigenvalue problem (4.12) can be rewritten as

[TABLE]

Here $\Lambda_{n}(\gamma)$ $(n=0,1,2,...)$, are a set of positive eigenvalues of the linear problem depending on $\gamma$ as a parameter. The eigenvalues $\Lambda_{n}$ can be ordered so that $\Lambda_{0}>\Lambda_{1}>\Lambda_{2}>...$, and larger $n$ correspond to eigenfunctions with larger number of nodes ($n=0$ eigenfunction has the largest scale). The needed values of $\gamma$ satisfy

[TABLE]

In the limit $\gamma\to 0$ frequencies $\sigma$ and $\sigma^{\prime}$ are vanishingly small, and the kernel

[TABLE]

is large, so many $\Lambda_{n}$ are greater than 1. On the other hand, for $\gamma\gg 1$

[TABLE]

and the kernel takes a form:

[TABLE]

Due to rapidly decreasing exponents the kernel is small, and thus $\Lambda_{n}$ are small. When $\gamma$ is changing from zero to infinity, many $\Lambda_{n}$ cross the unity value. Since for a given $\gamma$, the eigenfunction with the largest scale has the largest eigenvalue, $\Lambda_{0}$ will be the first to cross unity as $\gamma$ increases, so $\gamma_{0}>\gamma_{1}>\gamma_{2}>...\,$.

An equation analogous to (4.16) for odd $l$ has a negative kernel, ${\cal Q}_{\gamma}^{\rm odd}<0$. It means that all $\Lambda_{n}$ are negative for any $\gamma$, and no aperiodic solution is possible.

5 Numerical results

5.1 Aperiodic modes in the dispersed Agekyan model ($q=-1$)

For the dispersed Agekyan model

[TABLE]

which corresponds to $q=-1$ in (1.2) the integral equation (4.12) is reduced to an algebraic one,

[TABLE]

where $\sigma={\gamma}/{\Omega}$,

[TABLE]

It is known (Polyachenko & Shukhman, 2015) that this equation has only one even aperiodic solution $\gamma$, which is large when $\varpi$ is large. Thus, keeping in the hyperbolic functions the leading exponents only, one obtains the characteristic equation for even $l$ aperiodic modes,

[TABLE]

For this model, the function $\varpi(L_{T})$ can be approximated by the power law

[TABLE]

obtained numerically in the range $-3<\lg L_{T}<-1$. With this approximation formula, one can find solutions at arbitrary small $L_{T}$. Fig. 2 shows the dependence of the growth rate $\gamma$ of the unstable aperiodic solution for $l=2$. As expected, $\gamma$ is large as $L_{T}\to 0$ and it scales approximately as $\gamma\sim\varpi$ for very small $L_{T}$, and $\gamma\sim\varpi^{1/2}$ for $L_{T}\sim 0.1$. Recall that the dynamic frequency $\Omega=2.16$ is of the order unity, thus the obtained growth rates obey the inequality $|\gamma|>\Omega$ for $L_{T}<10^{-2}$.

5.2 Oscillatory modes in the Agekyan model

This approximation formula (5.5) allows us to calculate oscillatory unstable solutions in the form of even and odd spherical harmonics using

[TABLE]

for even $l$, and

[TABLE]

for odd $l$, where $\nu={\omega}/{\Omega}$ and

[TABLE]

The results for the first three harmonics in a wide range of frequencies and $\lg L_{T}=-2...-6$ are presented in panels of Fig. 3. As expected, the aperiodic solutions are absent for odd modes. The growth rates of aperiodic solutions (for $l=2$) rapidly increase as $L_{T}\to 0$, so for most values of $L_{T}$ the apriodic solutions are outside the (middle) panel.

Growth rates of the oscillatory solutions show weak dependence on $\textrm{Re}\,\omega$, especially for the smallest $L_{T}$ when $\gamma=\textrm{Im}\,\omega\approx 1$. Real parts of frequencies obey approximately $\textrm{Re}\,\omega/\Omega=n-1/4$ in the limit $L_{T}\to 0$, but these limiting values approach from different sides (in case of even $l$ – from the right, and in case of odd $l$ – from the left). In all cases the oscillatory solutions obey $|\omega|\gtrsim\Omega$.

5.3 Series $q=-1/2$

In this section we study a series of models with nontrivial dependence of the DF on the energy. For $q=-1/2$ the potential can be obtained in an analytical form for arbitrary $L_{T}$ (Polyachenko et al., 2013) and this explains our choice of $q$.

In the purely radial limit $\Phi=\ln r$, $\Omega=\sqrt{2\pi}e^{-E}$. If $L_{T}$ is small but finite, there is a small radius $r_{1}={\cal O}(L_{T})$ which separates two intervals. From $r_{1}$ to approximately 1, the potential is close to $\Phi=\ln r$, but in the interval $[0,r_{1}]$ it behaves like $-\sin{kr}/r$, with $k\propto L_{T}^{-1}$. However, for energies in the range $[E_{c},0]$, where $E_{c}=E_{c}(L_{T})$ is the energy of the particle on the circular orbit with angular momentum $L_{T}$ (see Fig. 4), the pericentre distance is not less than $r_{1}$, i.e. one can use the potential for the purely radial case to calculate the precession rate. The azimuth change $g(\alpha)$ during the pericentre passage is

[TABLE]

where $\alpha\equiv L/L_{\rm circ}(E)$, $L_{\rm circ}(E)=e^{E-1/2}$, $\mu=[\ln\left(1/\alpha\right)]^{-1}$ (Touma & Tremaine, 1997; Polyachenko et al., 2013) and thus the precession rate of the particle with energy $E$ and angular momenta $L$ is

[TABLE]

where

[TABLE]

However, this formula is valid for nearly radial orbits $E>E_{c}$ only and fails for the circular orbit $E=E_{c}$ where $\Omega_{\textrm{pr}}=0.293/L_{T}$. So, instead of (5.10), we shall calculate the precession rate numerically, using $\Phi=\ln r$ for the potential, and

[TABLE]

where

[TABLE]

and

[TABLE]

With new variables $z=e^{E}$, $z_{c}=\sqrt{e}\,L_{T}$ and $\nu=\omega/\Omega_{1}(z,L_{T})$, where

[TABLE]

the limiting integral equations (4.9) and (4.11) become

[TABLE]

for even $l$, and

[TABLE]

for odd $l$, where

[TABLE]

and the kernel is given by eq. (4.10).

The instability growth rates of the aperiodic modes $\omega=i\gamma$ ($l=2$) are presented in Fig. 5. Contrary to the Agekyan model (Fig. 2), in this case we have many aperiodic solutions.

Results of our calculations of oscillatory eigenmodes for $l=1...3$ spherical harmonics are presented in Fig. 6. Panels of the figure show both real and imaginary parts of $\omega$ of the first two modes versus the control parameter $L_{T}$. The difference between two successive real parts is $\approx 2.3$, and the overall behaviour resembles one of the oscillatory modes in the dispersed Agekyan model (see Fig. 3).

6 The orbital approach in systems with nearly radial stellar orbits

In this section we shall analyse validity of the orbital approach in studying systems with purely radial and nearly radial orbits. Recall that the orbital approach turns from consideration of a particle trajectory to precessing motion of the closed orbital wires. An angle between two successive apocentres of the particle on the radial orbit in the scale free potentials $\Phi\propto r^{s}$ is

[TABLE]

(Touma & Tremaine, 1997). The potentials of the softened polytropes in the limit of purely radial motion $L_{T}\to 0$ diverge in the centre as (Polyachenko et al., 2013)

[TABLE]

i.e. weaker than any negative power $s$, so $\delta\varphi=\pi$, and according to (3.6), $\Omega_{1}(E,L=0)=2\Omega_{2}(E,L=0)$. The precession rate of the nearly radial orbits $L\ll 1$,

[TABLE]

is slow, $\Omega_{\textrm{p}r}\ll\Omega_{1,2}$.

Now consider the motion of a particle on a nearly radial orbit in presence of a weak non-rotating slowly growing bar potential

[TABLE]

Switching to the action–angle variables in the orbital plane, ${\mathbf{x}}={\mathbf{x}}({\mathbf{I}},{\mathbf{w}})$, the small perturbation due to the bar can be written as the Fourier series over radial angle $w_{1}$ of the unperturbed orbit:

[TABLE]

The phases $lw_{1}+mw_{2}=(l\Omega_{1}+m\Omega_{2})t$ vary quickly for all $l$ except $l=-1$. Thus, omitting quickly oscillating terms one obtains an ‘averaged’ hamiltonian

[TABLE]

which possesses an adiabatic invariant $J\equiv I_{2}+I_{1}/2$ and ‘slow’ angle variable $\overline{w}_{2}=w_{2}-\frac{1}{2}\,w_{1}$ (Lynden-Bell, 1979; Polyachenko, 2004). The equations of motion are

[TABLE]

In particular, eq. (6.8) gives the change of the angular momentum perpendicular to the orbital plane, $I_{2}$, and eq. (6.10) describes the apsidal precession. The requirement of adiabaticity implies

[TABLE]

(1991) obtained an expression for the growth rate in the monoenergetic model ($q=-1$) with pure radial orbits in the framework of spoke approximation, when the orbital wires turn into spokes. In our notations the growth rate is

[TABLE]

where $\Omega$ is the radial frequency,

[TABLE]

$\rho_{\rm lin}(r)=1/|v_{r}|=1/\sqrt{2\Psi(r)}$ is a linear density of the spoke; \varpi\equiv\bigl{[}d\Omega_{\rm pr}/dL\bigr{]}_{L=0}. However, as we argued in Sections 3 and 5, the last parameter grows infinitely, as we turn to more and more radially anisotropic systems.

It is interesting to note that if we formally assume the scaled growth rate of the mode $\sigma=\gamma/\Omega$ to be small in Eq. (5.2) for monoenergetic model ($q=-1$), and use identity (4.15) for ${\cal F}_{l}(r,r^{\prime})$, we obtain exactly the same expression for the growth rate (6.12) found by (1991) in the spoke approximation. This fact justifies the spoke approximation for systems with sufficiently small $\varpi$ (moderately elongated orbits), but not for the very eccentric orbits! Note also that the growth rate $\gamma$ for $q=-1$ series scales as $\varpi^{1/2}$ for not too small $L_{T}$ in agreement with (6.12), but for very small $L_{T}$ grows even faster than $\varpi$.

Hence stability study of the spherical systems with nearly radial or purely radial orbits cannot be made in the framework of the orbital approach (and the spoke approximation in particular), since $\gamma$ grows with $\varpi$ and condition (6.11) fails.

7 Summary and Conclusions

Using a new technique based on integral eigenvalue equations, we reconsider here a well-known work on radial orbit instability by Antonov (1973) in which spherical models with purely radial motion are studied. The Antonov problem cannot be correctly solved in the purely radial models due to singularity in the centre. Thus series of models with parameter $L_{T}$ controlling orbit eccentricity including purely radial model (corresponding to $L_{T}=0$) should be used.

The derived integral equations involve an only large quantity $\varpi\equiv\Omega_{\textrm{pr}}(L_{T})/L_{T}$ in case of small $L_{T}$ tending to infinity as $L_{T}$ goes to zero. This quantity coincides with the Lynden-Bell derivative $[\partial\Omega_{\textrm{pr}}/\partial L]_{L=0}$ playing a crucial role in theory of radial orbit instability (Lynden-Bell, 1979).

We investigated stability of the spherically symmetric models with respect to perturbations $\propto\chi(r)\,P_{l}(\cos\theta)$ and obtained numerical solutions for two series of softened polytropic models $F(E,q)\propto H(L_{T}-L)(-2\,E)^{q}$ ($H(x)$ is the Heaviside function) allowing the purely radial limit (Polyachenko et al., 2013; Polyachenko & Shukhman, 2015).

The first one, $q=-1$, is a dispersed Agekyan model. The instability exists both for even and odd spherical harmonics $l$, for which multiple oscillatory modes with $\textrm{Re}\,\omega\approx(n-1/4)\,\Omega$ are found, where $n=1,2,...$; $\Omega\approx 2.16$ is the radial frequency of particles (in units $G=M=R=1$). The modes growth rates $\gamma\equiv\textrm{Im}\,\omega\approx 1$. Besides, we found aperiodic modes $\textrm{Re}\,\omega=0$ for even spherical harmonics with growth rates tending to infinity as $L_{T}\to 0$.

The second series $q=-1/2$ provides an analytic potential for any value of parameter $L_{T}$, and relatively simple formulae for the radial frequency and the precession rate for nearly radial models. As with the previous series, we found multiple oscillatory modes for even and odd spherical harmonics. A characteristic feature of this model is multiple aperiodic modes with growth rates increasing as $L_{T}\to 0$. We conclude that in all cases (both series, aperiodic and oscillatory modes, even and odd $l$) $|\omega|$ values are of the order of or larger than $\Omega$.

There are several interpretations for the physical mechanism of radial orbit instability (ROI). One relates ROI to the well-known Jeans instability in anisotropic medium for which insufficient velocity dispersion perpendicular to the radial direction cannot resist gravitational clusterization (Polyachenko & Shukhman, 1972). Another one is connected to precession dynamics of eccentric orbits that attract to each other, provided $[\partial\Omega_{\textrm{pr}}/\partial L]_{L=0}>0$ (Lynden-Bell, 1979; Merritt, 1987). This point of view can be justified only for the so-called ‘slow modes’ which satisfy ‘slow’ integral equation in which only one resonance term $\propto[2\Omega_{\textrm{pr}}-\omega]^{-1}$ is retained (Palmer, 1994). In turn, this implies (i) even $l$ only, and ‘slowness’ of the mode, i.e. $|\omega|$ should be much less than the dynamical frequencies, e.g. $\Omega$ (Polyachenko & Shukhman, 2015). As we saw, none of the solutions obtained in this work satisfy any of these requirements, and we must conclude that the orbital interpretation is limited.

Using the energy approach, Maréchal & Perez (2010) argue that instability in sufficiently anisotropic systems can be induced by dissipation inevitably present in the real stellar systems. The energy approach claims that if the second order variation of energy due to the perturbation, $H^{(2)}$, is negative, then the system may be unstable. If, in addition, a small dissipation takes place, the system is guaranteed to be unstable, with the growth rate proportional to the dissipation. Note, however, that the energy approach makes no conclusions for systems without dissipation in the case of negative sign of $H^{(2)}$. In other words, it is of little help for highly anisotropic spherical systems subject to very strong collisionless (i.e., non-dissipative) radial orbit instability, which is apparently more important than the instability potentially induced by dissipation.

Similar to Antonov (1973), we consider here non-radial perturbations independent of parity $l$, and an instability mechanism independent of the suggestion of slowness. However, strong singularity of central density inherent to the system with purely radial orbits (Bouvier & Janin, 1968; , 1984) leads to singularity of the potential, consequently infinite $\varpi$ and the growth rates $\gamma$ for even aperiodic modes. We suppose that this instability is manifestation of Jeans instability modified due to periodic radial motion of stars along their orbit.

It is worth recalling, in this context, the argument against our interpretation of ROI, raised for the first time by Merritt (1987). According to the virial theorem, the growth rate of Jeans instability is of the same order as the inverse crossing time, $\sim(G\rho)^{1/2}$. Since radially anisotropic systems are also strongly radially inhomogeneous, he claims that “unstable mode would scarcely begin to grow before the particles contributed to it had moved away from their initial positions, to regions of very different density and velocity dispersion”. The growth rates obtained in our calculations, however, are large compared to the inverse crossing time, which prevent particles from being escaped before the instability takes over the system. Thus, the virial estimate and the entire argument are not valid for the systems with orbits very close to purely radial.

Acknowledgments

We thank Dr. J. Perez for his comments when reviewing the paper, and Dr. R. Moetazedian for his help in improving the English language. This work was supported by the Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject A6) of the German Research Foundation (DFG), and by the Volkswagen Foundation under the Trilateral Partnerships grant No. 90411. The authors acknowledge financial support by the Russian Basic Research Foundation, grants 15-52-12387, 16-02-00649, and by Department of Physical Sciences of RAS, subprogram ‘Interstellar and intergalactic media: active and elongated objects’.

Appendix A The ‘proof’ of the instability using Lyapunov function for $l\gg 1$

In this Appendix we reproduce the original proof by Antonov, made with the aid of Lyapunov function for the case of large $l$, but in the notations and terms adopted in the present work.

When the radial derivatives of the perturbed potential can be neglected compared to the angular derivatives,

[TABLE]

or

[TABLE]

eqs. (2.13) and (2.14) can be simplified as

[TABLE]

where ${\hat{D}}=\nu(E)\,\partial/\partial w\ \ \textrm{and}\ \nu(E)=\Omega_{1}(E,L=0)$. The Poisson equation (2.22) then can be reduced to an algebraic equation:

[TABLE]

Now we introduce new variables $A$ and $B$:

[TABLE]

and the systems (A.1) and (A.2) can be rewritten as

[TABLE]

[TABLE]

Following Antonov (1973), we construct the Lyapunov function

[TABLE]

Differentiating over time, one obtains:

[TABLE]

Now rewriting (A.3)

[TABLE]

we have

[TABLE]

Since $B\,({\hat{D}}\,A)+A\,({\hat{D}}\,B)$ is a full derivative over $w$,

[TABLE]

the corresponding integral

[TABLE]

vanishes and we obtain

[TABLE]

or finally

[TABLE]

The last equation is the full analog of the Antonov’s expression for $dF/dt$ (but expressed in our variables).111 Note that in the cited paper by Antonov (1973) this expression (following eq. (7)) contains a misprint: the second term in the r.h.s. of the expression for $dF/dt$ should read: $2\pi G\int r^{2}\,dr\left[\int dE_{0}\rho_{E_{0}}\,(\xi_{+}+\xi_{-})\right]^{2}$.

The proof is based on the evident positiveness of both terms in (A.9). However, as we already noted in the main text, the first term diverges at $r=0$, so rigorously speaking such a proof of the radial orbit instability is invalid.

Appendix B Integral equations in the limit of small $L_{T}$ (even $l$).

Clarifying a sense of diverging coefficients $p_{k}$

In this appendix we restrict ourselves to the relatively compact derivation for the case of dispersed Agekyan model ($q=-1$), although generalisation to arbitrary $F(E)$ is possible. Besides, we shall consider even spherical harmonics $l$ only, but the desired relations used for interpretation of diverging integrals in the delta function technique are universal and valid for odd $l$ as well. So we assume

[TABLE]

where

[TABLE]

and starting from the integral equation in the Lagrange form,

[TABLE]

Here we denote

[TABLE]

where $\Theta_{l_{1}l_{2}}(E,L,w)$ is an angle,

[TABLE]

and

[TABLE]

is the azimuthal angle as the particle travels from pericentre to current radius $r$, $\Psi$ is the relative potential, $\Psi(r)\equiv-\Phi_{0}(r)$. In particular, in the apocentre ($w=\pi$) this angle is $\delta\varphi(E,L;\pi)=(\Omega_{2}/\Omega_{1})\,\pi$. The kernel functions are

[TABLE]

Note that the symmetry of the radial function $r(2\pi-w)=r(w)$ allows one to reduce integration in eqs (B.4) and (B.7) over full range of the angle variable to the interval $[0,\pi]$.

The r.h.s. of (B.3) can be divided into two parts

[TABLE]

where

[TABLE]

and

[TABLE]

Now we shall expand the integral equation entities on the small parameter $L_{T}$. The linear combination of frequencies can be rewritten through the precession rate $\Omega_{\rm pr}$,

[TABLE]

For $\Theta_{l_{1}\,l_{2}}(E,L;w)$ one can write

[TABLE]

where

[TABLE]

or

[TABLE]

Here terms proportional to $\Omega_{\rm pr}$ and $L$ are considered to be small and vanishing as $L$ approaches zero. To be clear, we assume the lower limit in the integral in (B.13) $(L/\Omega_{1})\int_{w}^{\pi}{dw^{\prime}}/{r^{2}(w^{\prime})}$ is not too close to zero, otherwise this integral becomes of the order unity, since for $w=0$ it equals to $\pi\,(\Omega_{2}/\Omega_{1})\approx\frac{1}{2}\,\pi$. However, the range of $w$ where the integral becomes $\sim 1$ is very small for $L\to 0$, and we shall see below that this bring no difficulties in further integrations. In (B.12) we take into account that the angle $\delta\varphi$ changes from zero to $\approx\pi/2$ in the centre, and then remains almost constant in the remaining part of the orbit. Angle $\beta$ is the remaining part of angle $\delta\varphi$ gained from $r$ to $r_{\rm max}$. Thus $\beta$ is small as long as $r\gg r_{\rm min}$ in (B.14), and contribution to $\delta\varphi$ gained near the centre is taken into account by the term $-\frac{1}{2}\,l_{2}\,\pi$ in the square brackets in (B.11).

For the even $l$, values of $l_{2}$ in the integral equation are even, so the sum $l_{1}+{\textstyle\frac{1}{2}}\,l_{2}$ is an integer. Introducing new indices

[TABLE]

one can switch in expressions for $Q_{1}$ and $Q_{2}$ from double summation over $l_{1}^{\prime}$ and $l_{2}^{\prime}$ to summation over $n^{\prime}$ and $l_{2}^{\prime}$,

[TABLE]

[TABLE]

For $\phi_{l_{1}l_{2}}$ one obtains, providing $\Omega_{\rm pr}$ and $\alpha$ are small,

[TABLE]

where

[TABLE]

and

[TABLE]

Similarly, for the kernel functions

[TABLE]

[TABLE]

and

[TABLE]

Since $\beta(w^{\prime})$ in (B.20) and (B.23) is multiplied by $\sin(n^{\prime}w^{\prime})$, which vanishes at $w=0$, the uncertainty in $\beta$ at $w\approx 0$ does not lead to any difficulties.

Now it is easy to relate eqs. (B.3) and (2.35). In the leading order over $L$, $\phi_{l_{1}l_{2}}$ coincides with $(-1)^{l_{2}/2}\Phi_{n}$ and the kernel functions $\Pi_{l_{1}l_{2};l_{1}^{\prime}l_{2}^{\prime}}$ coincide with ${\cal K}^{\rm even}_{nn^{\prime}}\cdot(-1)^{l_{2}/2+l_{2}^{\prime}/2}$ of eq. (2.35). Using the identities

[TABLE]

one can show that $Q_{E}$ turns into the last term containing the energy derivative. The remaining term, $Q_{L}$, vanishes in the leading order ${\cal O}(1/L_{T})$,

[TABLE]

because of the second identity in (B.24). To proceed further, we have to expand $Q_{L}$ to the next order ${\cal O}(L_{T}^{0})$ and compare it with the first square bracket in (2.35).

Small additional terms $(\delta\Phi)_{n^{\prime}}$ and $(\delta{\cal K})_{nn^{\prime}}$ can be expanded over functions of the leading order. According to (B.19)

[TABLE]

and from (B.20) one has

[TABLE]

where

[TABLE]

Similarly, for $(\delta{\cal K})_{nn^{\prime}}$ one obtains

[TABLE]

Summarising, for $L=0$, $L^{\prime}\ll 1$ one obtains:

[TABLE]

[TABLE]

[TABLE]

Using these expressions in (B.10), we obtain in the order ${\cal O}(L_{T}^{0})$:

[TABLE]

Summing up over $l_{2}^{\prime}$ with the help of (B.24) allows us to reduce (B.33) to

[TABLE]

The expression ${\bar{Q}}_{L}$ should be compared with the first term in r.h.s. of Eq. (2.36), which can be rewritten as

[TABLE]

In particular, comparison of the underbraced terms in (B.34) and (B.35) gives that $p_{0}$ should be associated with $\lim\limits_{L_{T}\to 0}\Omega_{\rm pr}(L_{T})/L_{T}$, i.e.

[TABLE]

Next, from (B.13) and (B.28)

[TABLE]

for $m\neq\pm n$ and

[TABLE]

for $m=\pm n$. Then, introducing

[TABLE]

one can have

[TABLE]

Now it is not difficult to show that ${\bar{Q}}_{L}$ completely coincides with ${Q}_{L}^{\rm pure\ radial}$ if one associate $p_{k}$ as limiting values of $P_{k}(L_{T})$:

[TABLE]

In particular, for $k=0$ one obtains (B.36)

[TABLE]

From definition of $p_{k}$,

[TABLE]

with the integral converging in the usual sense, thus $p_{k\neq 0}-p_{0}={\cal O}(1)$.