Skip to Main Content

Article Navigation

Journal Article

Stellar dynamics around a massive black hole – III. Resonant relaxation of razor-thin axisymmetric discs

Abstract

We study the resonant relaxation (RR) of an axisymmetric, low-mass (or Keplerian) stellar disc orbiting a more massive black hole (MBH). Our recent work on the general kinetic theory of RR is simplified in the standard manner by the neglect of ‘gravitational polarization’ and applied to a razor-thin axisymmetric disc. The wake of a stellar orbit is expressed in terms of the angular momenta exchanged with other orbits, and used to derive a kinetic equation for RR under the combined actions of self-gravity, 1 PN and 1.5 PN general relativistic effects of the MBH and an arbitrary external axisymmetric potential. This is a Fokker–Planck equation for the stellar distribution function (DF), wherein the diffusion coefficients are given self-consistently in terms of contributions from apsidal resonances between pairs of stellar orbits. The physical kinetics is studied for the two main cases of interest. (1) ‘Lossless’ discs in which the MBH is not a sink of stars, and disc mass, angular momentum and energy are conserved: we prove that general H-functions can increase or decrease during RR, but the Boltzmann entropy is (essentially) unique in being a non-decreasing function of time. Therefore, secular thermal equilibria are maximum entropy states, with DFs of the Boltzmann form; the two-ring correlation function at equilibrium is computed. (2) Discs that lose stars to the MBH through an ‘empty loss cone’: we derive expressions for the MBH feeding rates of mass, angular momentum and energy in terms of the diffusive fluxes at the loss-cone boundaries.

Galaxy: centre, galaxies: kinematics and dynamics, galaxies: nuclei

1 INTRODUCTION

Dense star clusters orbit massive black holes (MBHs) in galactic nuclei. Within the sphere of influence of the MBH, the inverse-square Newtonian gravity of the MBH is the dominant force on a star when its speed is non-relativistic. Were cluster self-gravity to be neglected completely, every star would orbit the MBH in a closed Keplerian ellipse. The orbits are not closed when the small perturbation due to self-gravity is taken into account; in addition to oscillations over the (fast) Kepler orbital time, the orbital elements of the ellipses also evolve over longer time-scales. This long-term behaviour is called secular evolution, for which there is a well-developed theory in planetary dynamics (Murray & Dermott 1999). In contrast to a planetary system, a low-mass star cluster orbiting an MBH is a many-body system, for which both the dynamics and the statistical mechanics are important. In Sridhar & Touma (2016a,b), hereafter referred to as Papers I and II, respectively, we presented a general theory of the secular evolution of these systems.

Let us consider a stellar system of size R with N ≫ 1 stars, each of mass m_⋆, orbiting an MBH of mass M_• . The typical Kepler orbital period, T_kep = 2π(R³/GM_•)^1/2, is the shortest time-scale in the problem. In common with other large N stellar systems, a Keplerian system may, in the first approximation, be treated in the ‘collisionless’ limit, N → ∞ and m_⋆ → 0, while M = Nm_⋆ is held constant. In this limit each star moves like a massless test particle that is acted upon by the MBH and the smooth mean gravitational field of the stellar system. For a low-mass stellar system, the total mass in stars, M = Nm_⋆, is much smaller than M_• ; the mass ratio, ε = M/M_• ≪ 1, is an important small parameter. The Keplerian orbital elements of every stellar orbit evolve over secular times, T_sec = ε⁻¹T_kep. In Paper I, we exploited the fact that there are two different time scales, T_kep and T_sec ≫ T_kep; thus, we began with the collisionless Boltzmann equation (CBE; Binney & Tremaine 2008), valid for a general stellar system, and averaged it over the fast Kepler orbital phase using the method of multiple time-scales. This resulted in a reduced description where each star is replaced by a Gaussian ring, which is a Keplerian ellipse whose eccentricity, inclination, apsidal angle and nodal longitude evolve over times T_sec, whereas the semi-major axis remains constant. The distribution function (DF) of Gaussian rings obeys a secular CBE, and hence many conclusions followed from general considerations: a secular Jeans theorem can be used to construct a rich variety of secular collisionless equilibria with greater ease than in the full, non-averaged dynamics; the linear instability of these equilibria is also simpler to study; and an arbitrary initial condition may be expected to undergo ‘violent’ (collisionless) relaxation to some coarse-grained secular equilibrium state, over times |${\sim }\, {\rm \text{a}\, few\,} T_{\rm sec}$|⁠.

When N is large but finite, as is the case in practice, a collisionless equilibrium state is realized only in an approximate sense. The small graininess inherent in the sampling of such a state will make the stellar system evolve in time due to gravitational ‘collisions’ between the stars. Two-body encounters between stars enable exchanges of both energy and angular momentum with a coherence time of the order of the orbital time T_kep. The incoherent sum of many such encounters results in kinetic evolution over the two-body relaxation time T_relax = ε⁻²NT_kep : some stars gain energy and increase their orbital radii while those that lose energy go into more tightly bound orbits, leading to a gravothermal catastrophe (Binney & Tremaine 2008). However T_relax is usually greater than the Hubble time for galactic nuclei of interest. Rauch & Tremaine (1996) proposed a more efficient process, called resonant relaxation (RR), which can be thought of as driven by gravitational interactions between pairs of Gaussian rings. Since the Keplerian energies (equivalently the semi-major axes) of Gaussian rings are conserved, they can only exchange angular momentum. This has a coherence time ∼ T_sec which is longer than the coherence time (T_kep) of two-body relaxation. The incoherent sum of many such encounters results in kinetic evolution over the RR time T_res = NT_sec = εT_relax , which can be shorter than the Hubble time for galactic nuclei of interest.¹ RR leads to a genuine equilibrium within a Hubble time (which can be expected to evolve very slowly due to two-body relaxation). In Paper II we presented a first-principles theory of RR, which is based on the rigorous kinetic theory of Gilbert (1968).

Gilbert's theory applies to general stellar systems and is based on a systematic 1/N expansion of the BBGKY hierarchy of equations:² the O(1) secular theory describes collisionless mean-field dynamics, and the collisional theory emerges at O(1/N). In order to derive a theory for RR in Paper II, we orbit-averaged Gilbert's kinetic equations by using the method of multiple scales that was introduced in Paper I. A resonantly relaxing stellar system can be thought of as passing through a sequence of quasi-stationary, collisionless equilibria. The RR theory of Paper II places no restriction on either the spatial geometry or the orbital structure of the stellar system, so both of these could evolve with time. Scalar and vector RR are naturally treated on par as contributions arising from the apsidal and nodal resonances between stellar orbits. But this very generality of our theory of RR required a formulation in terms of the abstract, yet compact, Poisson-Bracket notation. The goal of the present work is to make concrete this abstraction and derive explicit formulae that describe the physical kinetics of an astrophysically interesting system.

The model stellar system must be as simple as possible, so we are able to see in explicit form how the kinetic theory works; the system must also be realistic enough to describe interesting physical phenomena such as the feeding of mass, energy and angular momentum to the MBH. Hence, we begin by considering stellar systems whose DFs possess high spatial symmetry. The simplest three-dimensional DFs are those corresponding to spherical, non-rotating systems. These could be reasonable approximations to the stellar nuclei of many galaxies, and it is important to work out their RR kinetics and the rates at which mass and stellar orbital energy are fed to the MBH. However, for the simplest case of a non-spinning (Schwarzschild) MBH, the inflow of stars cannot carry, by symmetry, a net angular momentum. This limitation can be overcome only by relaxing the requirement of complete spherical symmetry: one could introduce some net rotation in the DF or/and consider a spinning (Kerr) MBH. But an exactly spherical and rotating system is somewhat unnatural, and the exclusion of Schwarzschild MBHs seems rather extreme. Hence, it is useful to study RR in non-spherical and rotating DFs. This can be done – and it is important to do so – but it comes at the expense of working with a more complicated orbital structure. We need a simpler system in order to display the physical kinetics of RR as transparently as possible. Our choice is an axisymmetric Keplerian disc of stars confined to a plane, with the spin of the MBH (when it is non zero) directed perpendicular to the disc plane. We note that this limits our study to scalar RR.

Discs are naturally rotating and we may expect RR to feed angular momentum, in addition to mass and energy, to the MBH through a loss cone. Discs around compact objects are ubiquitous in the Universe; a two-dimensional system is, of course, an idealization of a realistic disc but it is often useful in the study of disc dynamics (Binney & Tremaine 2008). In particular, the exploration of the (secular) statistical mechanics of two-dimensional Keplerian discs by Touma & Tremaine (2014) brings important insights into the nature of secular thermal equilibrium, whereas a corresponding detailed study does not exist for three-dimensional systems. They consider both axisymmetric discs and non-axisymmetric discs. In this paper, we make a beginning by providing the kinetic theory of the approach to secular thermal equilibrium of razor-thin axisymmetric discs.

In Section 2, we recall the salient results of Papers I and II for a general Keplerian system orbiting an MBH. This generality is necessary because we introduce in Section 3 the ‘Passive Response Approximation’ (PRA) in the context of RR (the PRA is a standard simplification of the kinetic equations for plasmas and gravitating systems, wherein gravitational polarization is ignored). From Section 4 onward we specialize to the case of a razor-thin axisymmetric disc. The wake produced by a stellar orbit is expressed in terms of the angular momenta exchanged between two stellar orbits, now viewed as two gravitationally interacting Gaussian rings. In Section 5 we obtain an explicit form of the collision integral and use this to write the RR kinetic equation in the ‘conservation’ form. The kinetic equation is then cast in the Fokker–Planck form and the diffusion coefficients identified. The collision integral, the RR probability current density and the diffusion coefficients get contributions only from orbits that are in apsidal resonance. Formulae for the total mass, energy and angular momentum of the stellar disc are derived. In Section 6, we consider a ‘lossless’ stellar disc – a disc in which there is no loss of stars to the MBH. We demonstrate that the kinetic equation conserves total mass, angular momentum and energy, and prove an H-theorem on the non-decreasing nature of the time evolution of the Boltzmann entropy. Therefore, RR drives the system towards a stationary state of maximum entropy, which is a secular thermal equilibrium described by a DF of the Boltzmann form. The RR probability current density vanishes at equilibrium, but particle correlations do not. We calculate the two-ring correlation function at secular thermal equilibrium. In Section 7, we consider the disc with a loss cone and derive explicit formulae for the rates at which mass, angular momentum and energy are fed to the MBH. Conclusions follow in Section 8.

2 SECULAR THEORY OF KEPLERIAN STELLAR SYSTEMS

Here, we provide a brief review of some of the results of Papers I and II, which are needed later for elaborating the RR of axisymmetric discs. We consider a model system composed of N ≫ 1 stars, each of mass m_⋆, orbiting an MBH of mass M_• . The stellar system is Keplerian when the dominant gravitational force on the stars is the inverse-square Newtonian force of the MBH. An isolated and only mildly relativistic system is Keplerian when its total mass in stars, M = Nm_⋆, is much smaller than M_• ; then the mass ratio ε = M/M_• ≪ 1 is a natural small parameter. The limiting case of negligible stellar self-gravity, ε → 0 , reduces to the Kepler problem of each star orbiting the MBH independently on a fixed Keplerian ellipse. An orbit is then conveniently described in six-dimensional phase space, in terms of the evolution of its Delaunay action-angle variables {I, L, L_z; w, g, h}. The three actions are |$I\,=\,\sqrt{GM_\bullet a\,}\,$|⁠, |$L\,=\, I\sqrt{1-e^2\,}$| the magnitude of the angular momentum and L_z = L cos i the z-component of the angular momentum. The three angles conjugate to them are, respectively, w the orbital phase, g the angle to the periapse from the ascending node and h the longitude of the ascending node. The Kepler orbital energy E_k(I) = −1/2(GM/I)² depends only on the action I. Hence, when ε = 0, all the Delaunay variables except w are constant in time; w itself advances at the (fast) rate 2π/T_kep = (GM_•/a³)^1/2.

When ε ≪ 1 and non-zero, the self-gravity is small but its effects build up over long times. Secular theory seeks to describe the average behaviour of quantities over (secular) times T_sec = ε⁻¹T_kep , by averaging over the fast phase w. It should be noted that the stellar system must be stable on the fast time T_kep in order for the secular description to make sense (see footnote 4 of Paper I). Since w is no longer present in the secular theory, its conjugate action I is a conserved quantity; hence the semi-major axis of every stellar orbit is constant in time. Then, the stellar system can be described as a self-gravitating system of N Gaussian rings. Each ring is a point in five-dimensional ring space with coordinates |${\cal R}\equiv \lbrace I, L, L_z; g, h\rbrace$|⁠, representing a Keplerian ellipse of given semi-major axis, eccentricity, inclination, periapse angle and nodal longitude. A ring can deform over secular times T_sec , under the combined action of mutual self-gravity between the rings, secular relativistic corrections and external gravitational fields (which can vary over times T_sec). This is represented by an orbit |${\cal R}(\tau )$|⁠, with I = constant, parametrized by the slow time variable τ = ε(time).

2.1 Collisionless dynamics

Over times approximately several T_sec ≪ NT_sec , the stellar system can be thought of as ‘collisionless’. As noted in the Introduction, the collisionless limit corresponds to assuming that the system is composed of an infinite number of stars, each of infinitesimal mass, and the whole having a mass M equal to the total stellar mass (N → ∞ , m_⋆ → 0 with M = Nm_⋆ held constant). The Keplerian stellar system is described by |$\,F({\cal R}, \tau )$|⁠, the ring DF. When the MBH is not considered to be a sink of stars, the DF is normalized as

\begin{equation} \int F({\cal R}, \tau )\,{\rm d}{\cal R}\,=\, 1, \end{equation}

(1)

so we may regard F as a probability DF in five-dimensional ring space. Ring orbits are governed by the secular Hamiltonian:

\begin{equation} H({\cal R}, \tau ) \,=\, \Phi ({\cal R}, \tau ) \,+\, H^{\rm rel}(I, L, L_z) \,+\, \Phi ^{\rm tid}({\cal R}, \tau ), \end{equation}

(2)

which is the sum of three terms accounting for self-gravity, general relativistic corrections to the MBH's field, and external gravitational sources. Since time is measured by the slow time variable ε(time), each of the three pieces is scaled by ε⁻¹. Thus, the mean-field potential of all the rings is |$\,\varepsilon \Phi ({\cal R}, \tau )$|⁠, where

\begin{eqnarray} \Phi ({\cal R}, \tau ) \,=\, \int F({\cal R}^{\prime }, \tau )\,\Psi ({\cal R}, {\cal R}^{\prime }) \,{\rm d}{\cal R}^{\prime },\quad \mbox{ring mean-field potential} \end{eqnarray}

(3a)

\begin{eqnarray} \Psi ({\cal R}, {\cal R}^{\prime }) \,=\, -GM_\bullet \oint \oint \frac{{\rm d}w}{2\pi }\, \frac{{\rm d}w^{\prime }}{2\pi }\,\frac{1}{\left|{r} - {r}^{\prime }\right|}, \quad \mbox{`bare' inter-ring potential} \end{eqnarray}

(3b)

is determined self-consistently by the DF itself.³ The second term, H^rel, determines the general relativistic 1 PN (post-Newtonian) Schwarzschild precession of the apses and 1.5 PN Lense–Thirring precession of the apses and nodes. The third term, Φ^tid, is the scaled tidal field due to an external source – see equations (30)–(34) of Paper I for more details on all three terms. The motion of each ring, |${\cal R}(\tau )$|⁠, is governed by the following Hamiltonian equations:

\begin{eqnarray} &&{I = \sqrt{GM_\bullet a} = \mbox{constant},} \nonumber \\ &&{\frac{{\rm d}L}{{\rm d}\tau } = -\,\frac{\mathrm{\partial} H}{\mathrm{\partial} g},\quad \frac{{\rm d}g}{{\rm d}\tau } \,=\, \frac{\mathrm{\partial} H}{\mathrm{\partial} L}\, \quad \frac{{\rm d}L_z}{{\rm d}\tau } \,=\, -\,\frac{\mathrm{\partial} H}{\mathrm{\partial} h},\quad \frac{{\rm d}h}{{\rm d}\tau } \,=\, \frac{\mathrm{\partial} H}{\mathrm{\partial} L_z}. } \end{eqnarray}

(4)

During collisionless evolution, the DF is constant along the orbit of every ring. Hence, |$F({\cal R}, \tau )$| obeys the secular CBE:

\begin{equation} \frac{{\rm d}F}{{\rm d}\tau } \,\equiv \, \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,+\, \left[\,F,\,H\,\right] \,=\, 0, \end{equation}

(5)

where [ , ] is the four-dimensional Poisson Bracket,

\begin{equation} \left[\,\chi _1,\,\chi _2\,\right] \,\stackrel{{\rm def}}{=}\, \left(\frac{\mathrm{\partial} \chi _1}{\mathrm{\partial} g}\frac{\mathrm{\partial} \chi _2}{\mathrm{\partial} L} - \frac{\mathrm{\partial} \chi _1}{\mathrm{\partial} L}\frac{\mathrm{\partial} \chi _2}{\mathrm{\partial} g}\right) \,+\, \left(\frac{\mathrm{\partial} \chi _1}{\mathrm{\partial} h}\frac{\mathrm{\partial} \chi _2}{\mathrm{\partial} L_z} - \frac{\mathrm{\partial} \chi _1}{\mathrm{\partial} L_z}\frac{\mathrm{\partial} \chi _2}{\mathrm{\partial} h}\right), \end{equation}

(6)

whose action is restricted to the four-dimensional I = constant surfaces in the five-dimensional |${\cal R}$|-space. An important general property, valid for quite general |$F({\cal R}, \tau )$| and |$H({\cal R}, \tau )$|⁠, is as follows: the probability for a ring to be in (I, I + dI) is a conserved quantity. In other words, the PDF in one-dimensional I-space, defined by

\begin{equation} P(I) \,=\, \int \, \text{d}L\,\text{d}L_z\,\text{d}g\,\text{d}h\, F(I,\, L,\, L_z,\, g,\, h,\,\tau ), \end{equation}

(7)

is independent of τ, as can be verified directly using the ring CBE equation (5).

Secular collisionless equilibria |$F= F_0({\cal R})$| are stationary (i.e. τ-independent) solutions of the ring CBE equation (5) and satisfy |$[F_0({\cal R}),\,H_0({\cal R})] = 0\,$|⁠. These can be constructed by a secular Jeans’ theorem which states that F₀ is a function of |${\cal R}$| only through the time-independent integrals of motion of |$H_0({\cal R})$|⁠, and any (positive and normalized) function of the time-independent integrals of |$H_0({\cal R})$| is a stationary solution of equation (5). Linear response and stability properties of a secular equilibrium |$F_0({\cal R})$| can be studied by looking at the behaviour of small perturbations. For a small perturbation |$F_1({\cal R}, \tau )$|⁠, the change in the self-gravitational potential is εΦ₁, which is related to F₁ through the Poisson integral equation (3a):

\begin{equation} \Phi _1({\cal R}, \tau ) \,=\, \int {\rm d}{\cal R}^{\prime }\,F_1({\cal R}^{\prime }, \tau )\,\Psi ({\cal R}, {\cal R}^{\prime }). \end{equation}

(8)

As discussed in Paper I, one can also allow for a non-zero change in the external tidal field |$\Phi _1^{\rm tid}$|⁠. When |$\Phi _1^{\rm tid} = 0,$| the perturbed Hamiltonian H₁ = Φ₁ is entirely due to self-gravity. Substituting F = F₀ + F₁ and H = H₀ + Φ₁ into equation (5), linearizing and using equation (8) for Φ₁, we obtain the following linearized CBE (or LCBE) determining the secular stability of an equilibrium |$F_0({\cal R})$|⁠:

\begin{equation} \frac{\mathrm{\partial} F_1}{\mathrm{\partial} \tau } \,+\, \left[\,F_1,\,H_0\,\right] \,+\, \int {\rm d}{\cal R}^{\prime }\left[\,F_0,\,\Psi ({\cal R}, {\cal R}^{\prime })\,\right]\,F_1({\cal R}^{\prime }, \tau ) \,\,=\,0. \end{equation}

(9)

The equilibrium is said to be secularly stable if there are no growing solutions to equations (9).⁴ A Fourier transform in time gives a linear integral eigenvalue problem for secular stability that is much simpler than the corresponding problem for full stability. In Paper I, a stability result was proved for razor-thin axisymmetric discs; this is discussed further in Section 3 and later.

2.2 Resonant relaxation

Over times ∼ T_res = NT_sec, the small O(1/N) effects, due to ring–ring gravitational encounters (‘collisions’), build up and cannot be neglected. The evolution of |$F({\cal R}, \tau )$| is now governed by the ring kinetic equation in the BBGKY form:

\begin{equation} \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,+\, \left[\,F,\, H - \frac{\Phi ({\cal R}, \tau )}{N}\,\right] \,=\, {\cal C}[F], \end{equation}

(10)

where

\begin{equation} {\cal C}[F] \,=\, \frac{1}{N}\int \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right]\,{\rm d}{\cal R}^{\prime } \end{equation}

(11)

is the ‘collision integral’ and |$(1/N)F^{(2)}_{\rm irr}$| is the irreducible part of the two-ring correlation function. When compared with the secular CBE of equation (5), two new O(1/N) terms appear in the kinetic equation (10). The first term, Φ/N, accounts for the fact that the mean-field acting on any given ring is due only to the other (N − 1) rings. The second term, |${\cal C}[F]$|⁠, drives irreversible evolution due to RR. We briefly note how these terms arise from the BBGKY hierarchy: the equation for the one-particle DF depends on the two-particle DF. Gilbert (1968) expressed the two-particle DF as the sum of a ‘reducible’ part and an ‘irreducible’ part. The reducible part is equal to the product of two one-particle DFs which, evidently, has no information on particle correlations. This corrects the mean-field potential by accounting for the fact that only (N − 1) particles act on any given particle. The irreducible part contains information about correlations between particles that develop due to gravitational ‘collisions’, which eventually drives kinetic evolution. Upon orbit-averaging Gilbert's equations, these give (see Paper II) the terms Φ/N and |$(1/N)F^{(2)}_{\rm irr}$|⁠.

The kinetic equation applies to the collisional evolution of F whenever it satisfies

\begin{equation} \left[\,F({\cal R}, \tau ),\, H - \frac{\Phi ({\cal R}, \tau )}{N}\,\,\right] \,\simeq \, 0,\quad \mbox{being of the order of $1/N$ or smaller.} \end{equation}

(12)

In other words, |$F({\cal R}, \tau )$| may be thought of as passing through a sequence of secular collisionless equilibria in a quasi-static manner. The dependence on τ is very slow, occurring over the RR time-scale T_res = NT_sec . Each of these equilibria must be secularly stable, i.e. the LCBE equation (9) for linear perturbations about |$F({\cal R}, \tau )$| must not have growing modes.⁵ For equation (10) to be a closed kinetic equation for the DF |$F({\cal R}, \tau )$|⁠, we still need to specify |$F^{(2)}_{\rm irr}$| in terms of F. The BBGKY equation determining |$F^{(2)}_{\rm irr}$| depends on the three-ring DF. In general, the three-ring DF satisfies an equation that depends on the four-ring DF. More precisely, it is the irreducible part of the three-ring DF that depends on the four-ring DF; the reducible part depends only on the one-particle and two-particle DFs and does not need information on the four-ring DF. It turns out that the reducible part is O(1/N), whereas the irreducible part is O(1/N²), which is of higher order than the theory aims to be. Hence, the BBGKY hierarchy closes at the level of the two-ring DF; in other words, the equations for F and |$F^{(2)}_{\rm irr}$| form a closed system.

However, the equation for |$F^{(2)}_{\rm irr}$| depends on the 11 independent variables |${\cal R}$|⁠, |${\cal R}^{\prime }$| and τ in a complicated manner. The resolution of this difficulty lies in the fact that the irreducible part of the two-particle correlator can be expressed in terms of the linear response of the purely collisionless problem, which involves fewer independent variables. We describe this briefly below in the context of secular dynamics. Consider the following gedanken experiment due to Rostoker (1964) and Gilbert (1968): from the N rings distributed according to the DF F, select one ring and place it an orbit |${\cal R}^{\prime }(\tau ^{\prime })$| such that it arrives at the location |${\cal R}^{\prime }$| at time τ΄ = τ. This small perturbation to the system will induce an additional response throughout the stellar system. We write the net perturbation to F, at any instant τ΄ < τ, as

\begin{equation} F_1({\cal R}, \tau ^{\prime }) \,=\, -\frac{F({\cal R}, \tau ^{\prime })}{N} \,+\, \frac{\delta ({\cal R}- {\cal R}^{\prime }(\tau ^{\prime }))}{N} \,+\, \frac{W({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime })}{N}, \end{equation}

(13)

where |${\cal R}^{\prime }(\tau ^{\prime })$| is the location at time τ΄ of the ring which arrives at |${\cal R}^{\prime }$| at time τ. The first term on the right side can be thought of as the decrement suffered by F when one ring is removed from the system, and the second term accounts for the placing of this ring in a particular orbit. The last term is the collisionless response of the system: the quantity, |$W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )$|⁠, will be referred to as the wake at location |${\cal R}$|⁠, at time τ, due to a ring located at |${\cal R}^{\prime }$|⁠. Since the experiment does not alter the total mass of the system, there must be no net mass either in F₁ or in W. Also the average value of W at any location |${\cal R}$| must vanish. Hence, the wake satisfies the two identities:

\begin{eqnarray} \int W({\cal R}, {\cal R}^{\prime }, \tau )\,{\rm d}{\cal R}\,&=& 0,\quad \mbox{zero mass in the wake of ring {\cal R}^{\prime }\,}\nonumber \\ \int W({\cal R}, {\cal R}^{\prime }, \tau )F({\cal R}^{\prime }, \tau )\,{\rm d}{\cal R}^{\prime } \,&=& 0, \quad \mbox{zero net wake at {\cal R} due to all the rings.} \end{eqnarray}

(14)

The ring selection in the Rostoker–Gilbert experiment is imagined to be made in the distant past so that the system experiences an ‘adiabatic’ stimulus. Hence, the wake also builds up gradually from zero in the distant past.

Since N ≫ 1 the perturbation is effectively infinitesimal, so F₁ must satisfy the LCBE equation (9). Substituting equation (13) into (9), we obtain the following ring wake equation:

\begin{eqnarray} \frac{\mathrm{\partial} W}{\mathrm{\partial} \tau ^{\prime }} + \left[\,W({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }),\,H({\cal R}, \tau ^{\prime })\,\right] \,+\, \left[\,F({\cal R}, \tau ^{\prime }),\,\Phi ^{\rm w}({\cal R}, {\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime })\,\right] = \left[\,\Phi ^{\rm p}({\cal R}, {\cal R}^{\prime }(\tau ^{\prime }),\tau ^{\prime }),\,F({\cal R}, \tau ^{\prime })\,\right] ,\quad \mbox{for $\tau ^{\prime } \le \tau $.} \end{eqnarray}

(15)

Here, Φ^w is the gravitational potential due to the wake:

\begin{equation} \Phi ^{\rm w}({\cal R}, {\cal R}^{\prime }, \tau ^{\prime }) \,=\, \int \,W({\cal R}^{\prime \prime }\,\vert \,{\cal R}^{\prime }, \tau ^{\prime })\,\Psi ({\cal R}, {\cal R}^{\prime \prime })\,{\rm d}{\cal R}^{\prime \prime },\quad \mbox{ring wake potential}. \end{equation}

(16)

Φ^p is the difference between the ‘bare’ inter-ring interaction potential and the mean-field potential of equations (3a) and (3b):

\begin{equation} \Phi ^{\rm p}({\cal R}, {\cal R}^{\prime }, \tau ^{\prime }) \,=\, \Psi ({\cal R}, {\cal R}^{\prime }) \,-\, \Phi ({\cal R}, \tau ^{\prime }),\quad \mbox{ring perturbing potential}. \end{equation}

(17)

Equation (15) should be solved with the ‘adiabatic turn-on’ initial condition, |$W({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }) \rightarrow 0\,$| as τ΄ → 0 . The right side of equation (15) is the ‘source term’ for the wake because if it were absent W = 0 would be a solution for all time that is compatible with the initial condition. The third term on the left side depends on the ring wake potential, Φ^w, which is determined in equation (16) as an integral of W over |${\cal R}$|-space. Hence, equation (15) for W is an inhomogeneous linear integro-differential equation in the six independent variables |${\cal R}$| and τ. This is, of course, difficult to solve in full generality; in Section 3, we discuss a method of producing iterative solutions.

Extending the work of Rostoker (1964) in homogeneous plasmas to inhomogeneous stellar systems, Gilbert (1968) proved, using non-trivial operator identities, that the wake function, W, completely determines the irreducible part of the two-ring correlation function, |$F^{(2)}_{\rm irr}$|⁠. In the secular case this can be written as

\begin{eqnarray} F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )\,F({\cal R}^{\prime }, \tau ) \,+\, W({\cal R}^{\prime }\,\vert \,{\cal R}, \tau )\,F({\cal R}, \tau )\,+\, \int W({\cal R}\,\vert \,{\cal R}^{\prime \prime }, \tau )\,W({\cal R}^{\prime }\,\vert \,{\cal R}^{\prime \prime }, \tau )\,F({\cal R}^{\prime \prime }, \tau )\,{\rm d}{\cal R}^{\prime \prime }. \end{eqnarray}

(18)

This decomposition of |$F^{(2)}_{\rm irr}$| allows for a nice physical interpretation: the irreducible part of the two-ring correlator at |$({\cal R}, {\cal R}^{\prime })$| is the sum of three different contributions from the wake function: (a) the wake of |${\cal R}^{\prime }$| at |${\cal R}$|⁠; (b) the wake of |${\cal R}$| at |${\cal R}^{\prime }$|⁠; (c) the product of the wake values at the points |${\cal R}$| and |${\cal R}^{\prime }$| of a third ring at |${\cal R}^{\prime \prime }$|⁠, summed over all locations |${\cal R}^{\prime \prime }$|⁠. All three contributions come with suitable F-weighting. The third term accounts for the contribution of collective effects (‘gravitational polarization’) to the microscopic processes driving RR. This is the nature of the full theory at O(1/N). Substituting equation (18) for |$F^{(2)}_{\rm irr}$| into the kinetic equations (10) and (11), and manipulating the resulting expressions, we can also write the kinetic equation in a form where the dissipative and fluctuating contributions to the ring collision term are displayed explicitly:⁶

\begin{eqnarray} \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,+\, \left[\,F,\, H - \frac{\Phi ({\cal R}, \tau )}{N}\,\right] \,=\, {\cal C}^{\rm dis}[F] \,+\, {\cal C}^{\rm fluc}[F], \end{eqnarray}

(19a)

\begin{eqnarray} {\cal C}^{\rm dis}[F] \,=\, \frac{1}{N}\int \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,F({\cal R}, \tau )W({\cal R}^{\prime }\,\vert \,{\cal R}, \tau )\,\right]\,{\rm d}{\cal R}^{\prime }, \end{eqnarray}

(19b)

\begin{eqnarray} {\cal C}^{\rm fluc}[F] \,=\, \frac{1}{N}\int \,F({\cal R}^{\prime }, \tau )\left[\,\Psi ({\cal R}, {\cal R}^{\prime }) \,+\, \Phi ^{\rm w}({\cal R}, {\cal R}^{\prime }, \tau ),\,W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )\,\right]{\rm d}{\cal R}^{\prime }. \end{eqnarray}

(19c)

It is necessary to specify boundary conditions for the DF and the wake. We consider two cases, one in which the stellar system does not lose stars to the MBH, and the other in which stars are lost to the MBH.

No-loss boundary conditions: when the MBH is not considered to be a sink of stars, the set of N rings is a closed system and the ring DF is normalized as |$\int F({\cal R}, \tau )\,{\rm d}{\cal R}=\ 1\,$|⁠. Subject to this normalization |$F({\cal R}, \tau )$| can take any positive value at any location in |${\cal R}$|-space; in other words, the domain of F is all of |${\cal R}$|-space. Similarly, the domain of the wake function |$W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )$| is all of |${\cal R}$| and |${\cal R}^{\prime }$| spaces.
Empty loss-cone boundary conditions: a star that comes close enough to the MBH will be either torn apart by its tidal field or swallowed whole by it. We can only require that at some initial time τ₀,
\begin{equation} \int F({\cal R}, \tau _0)\,{\rm d}{\cal R}\,=\, 1,\quad \mbox{initial normalization.} \end{equation}
(20)
At later times τ > τ₀, we will have |$\int F({\cal R}, \tau )\,{\rm d}{\cal R}\,\le \, 1\,$|⁠. The loss of stars at small distances from the MBH can be simply modelled by an absorbing barrier: we assume that a star is lost to the MBH when its pericentre distance is smaller than some fixed value r_lc, the ‘loss-cone radius’. When the loss cone is empty, stars belonging to the cluster must necessarily have pericentre radii a(1 − e) larger than r_lc . Hence the domain of |$F({\cal R}, \tau )$| is restricted to regions of |${\cal R}$|-space in which I and L are large enough:
\begin{eqnarray} &&{I_{\rm lc} < I \quad \mbox{where}\quad I_{\rm lc} \,=\, \sqrt{GM_\bullet r_{\rm lc}}}\nonumber \\ &&{L_{\rm lc}(I) < L \,\le \, I, \quad \mbox{where}\quad L_{\rm lc}(I) \,=\, I_{\rm lc}\left[\,2 \,-\, \left(\frac{I_{\rm lc}}{I}\right)^2\,\right]^{1/2}.} \end{eqnarray}
(21)
The boundary conditions on the ring DF and wake functions are
\begin{eqnarray} F({\cal R}, \tau ) &=& 0\quad \mbox{for $I \le I_{\rm lc}\,$ and $\,L \le L_{\rm lc}(I)$,}\nonumber \\ W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,&=& 0\quad \mbox{for $I, I^{\prime } \le I_{\rm lc}\,$ and $\,L \le L_{\rm lc}(I)$, $\,L^{\prime } \le L_{\rm lc}(I^{\prime })$.} \end{eqnarray}
(22)

Equation (15) for W and the kinetic equation for F – either equations (10)–(11) and (18) or (19a)–(19c) – together with suitable boundary conditions are the fundamental equations governing RR.

3 PASSIVE RESPONSE APPROXIMATION

The coupled equations for the DF and the wake, together with either the No-loss or Empty loss-cone boundary conditions, offer a closed and consistent description of RR. One way to make progress is to try to solve equation (15) and obtain W as a function of F. We introduce the ‘convective derivative’, |$\,{\rm d} W/d\tau ^{\prime } \equiv {\mathrm{\partial} W}/{\mathrm{\partial} \tau ^{\prime }} \,+\, \left[\,W({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }),\,H({\cal R}, \tau ^{\prime })\,\right]\,$|⁠, which is the time derivative while following the ring orbits |${\cal R}(\tau ^{\prime })$| and |${\cal R}^{\prime }(\tau ^{\prime })$|⁠. Then equation (15) and the initial condition can be written as

\begin{eqnarray} \frac{{\rm d}W}{{\rm d}\tau ^{\prime }} \,=\, \left[\,\Phi ^{\rm p}({\cal R}, {\cal R}^{\prime }(\tau ^{\prime }),\tau ^{\prime }),\,F({\cal R}, \tau ^{\prime })\,\right] \,+\, \int {\rm d}{\cal R}^{\prime \prime }\left[\,\Psi ({\cal R}, {\cal R}^{\prime \prime }),\,F({\cal R}, \tau ^{\prime })\,\right]\,W({\cal R}^{\prime \prime }\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }), \nonumber \\ [1ex] \mbox{with adiabatic turn-on initial condition}\quad \lim _{\tau ^{\prime }\rightarrow -\infty } W({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }) \,=\, 0. \end{eqnarray}

(23)

The initial condition on W implies that at early times W is small, so the time rate of change of W is dominated by the first term (i.e. the source term) on the right side. Then we can expand the wake perturbatively as

\begin{equation} W \,=\, W^{(1)} \,+\, W^{(2)} \,+\, W^{(3)} \,+\, \cdots , \end{equation}

(24)

and substitute this into equation (23) to get

\begin{eqnarray} \frac{{\rm d}W^{(1)}}{{\rm d}\tau ^{\prime }} \,=\, \left[\,\Phi ^{\rm p}({\cal R}, {\cal R}^{\prime }(\tau ^{\prime }),\tau ^{\prime }),\,F({\cal R}, \tau ^{\prime })\,\right], \quad \mbox{with}\quad \lim _{\tau ^{\prime }\rightarrow -\infty } W^{(1)}({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }) \,=\, 0\, \end{eqnarray}

(25a)

\begin{eqnarray} &&{\frac{{\rm d}W^{(n+1)}}{{\rm d}\tau ^{\prime }} \,=\, \int {\rm d}{\cal R}^{\prime \prime }\left[\,\Psi ({\cal R}, {\cal R}^{\prime \prime }),\,F({\cal R}, \tau ^{\prime })\,\right]\,W^{(n)}({\cal R}^{\prime \prime }\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }),} \nonumber \\ &&{\quad \mbox{with}\quad \lim _{\tau ^{\prime }\rightarrow -\infty } W^{(n+1)}({\cal R}\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }) \,=\, 0, \quad \mbox{for $n=1,2,\ldots \,$.} } \end{eqnarray}

(25b)

The physical meaning of these equations is as follows: the lowest order wake, W⁽¹⁾, is driven directly by the source term, [Φ^p , F], and the role of the self-gravity of the wake on its own evolution is neglected. The stellar system can be thought of as responding passively to the forcing by the source term, and we will refer to this as the ‘PRA’, and to W⁽¹⁾ as the PRA wake function. As W⁽¹⁾ builds up from zero, its gravity acts as the source for W⁽²⁾, and so on; in general W⁽ⁿ⁾ acts as the source for W^{(n + 1)}. These higher order wakes, W⁽²⁾, W⁽³⁾ etc., correct for the neglect of gravitational polarization in the PRA. Direct integration of equations (25b) and (25b) gives the following explicit hierarchy of solutions:

\begin{eqnarray} W^{(1)}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, \int _{-\infty }^{\tau } \left[\,\Phi ^{\rm p}({\cal R}(\tau ^{\prime }), {\cal R}^{\prime }(\tau ^{\prime }),\tau ^{\prime }),\,F({\cal R}(\tau ^{\prime }), \tau ^{\prime })\,\right]\,{\rm d}\tau ^{\prime }, \end{eqnarray}

(26a)

\begin{eqnarray} W^{(n+1)}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, \int _{-\infty }^{\tau }{\rm d}\tau ^{\prime }\,\int {\rm d}{\cal R}^{\prime \prime }\left[\,\Psi ({\cal R}(\tau ^{\prime }), {\cal R}^{\prime \prime }),\,F({\cal R}(\tau ^{\prime }), \tau ^{\prime })\,\right]\,W^{(n)}({\cal R}^{\prime \prime }\,\vert \,{\cal R}^{\prime }(\tau ^{\prime }), \tau ^{\prime }),\quad \mbox{for $n = 1,2,\ldots \,$.} \end{eqnarray}

(26b)

We recall from equation (12), and the discussion following it, that the full RR kinetic theory applies to systems that evolve in a quasi-stationary manner through a sequence of secularly stable (collisionless) equilibria. This property is inherited by the above PRA integral expressions for wakes of various orders. Here F should be regarded as a slowly varying function of τ΄, with time-scale ∼ T_res ≫ T_sec ; this fact is used later in Section 4.2, to simplify in a standard manner the time integral for W⁽¹⁾.

The expansion for the wake, equation (24), can be substituted into equation (18), and |$F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )$| can be expanded to different orders. The lowest order,

\begin{equation} F^{(2,1)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, W^{(1)}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )\,F({\cal R}^{\prime }, \tau ) \,+\, W^{(1)}({\cal R}^{\prime }\,\vert \,{\cal R}, \tau )\,F({\cal R}, \tau ), \end{equation}

(27)

is the PRA form of the irreducible part of the two-ring correlation. When compared with equation (18) we see that the polarization term has dropped out because it is quadratic in W⁽¹⁾ . The two contributions to |$F^{(2,1)}_{\rm irr}$| come from the wake of |${\cal R}^{\prime }$| at |${\cal R}$|⁠, and the wake of |${\cal R}$| at |${\cal R}^{\prime }$|⁠. The PRA kinetic equation in the BBGKY form is

\begin{equation} \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,+\, \left[\,F,\, H - \frac{\Phi ({\cal R}, \tau )}{N}\,\right] \,=\, {\cal C}[F] \,=\, \frac{1}{N}\int \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,1)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right]\,{\rm d}{\cal R}^{\prime }. \end{equation}

(28)

When equation (27) is substituted into (28), we can express the PRA kinetic equation in the fluctuation–dissipation form:

\begin{eqnarray} \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,+\, \left[\,F,\, H - \frac{\Phi ({\cal R}, \tau )}{N}\,\right] \,=\, {\cal C}[F] \,=\, {\cal C}^{\rm dis}[F] \,+\, {\cal C}^{\rm fluc}[F], \end{eqnarray}

(29a)

\begin{eqnarray} {\cal C}^{\rm dis}[F] \,=\, \frac{1}{N}\int \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,F({\cal R}, \tau )W^{(1)}({\cal R}^{\prime }\,\vert \,{\cal R}, \tau )\,\right]\,{\rm d}{\cal R}^{\prime }, \end{eqnarray}

(29b)

\begin{eqnarray} {\cal C}^{\rm fluc}[F] \,=\, \frac{1}{N}\int \,F({\cal R}^{\prime }, \tau )\left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,W^{(1)}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )\,\right]{\rm d}{\cal R}^{\prime }. \end{eqnarray}

(29c)

Comparing with equations (19a)–(19c), we see that, as expected, the wake potential Φ^w has dropped out of the right side of equation (19c) for |${\cal C}^{\rm fluc}[F]$|⁠. Similar to the general case, equations (14), the PRA wake function satisfies the two identities,

\begin{eqnarray} \int W^{(1)}({\cal R}, {\cal R}^{\prime }, \tau )\,{\rm d}{\cal R}\,&=& 0,\quad \mbox{zero mass in the wake of ring ${\cal R}^{\prime }}\nonumber \\ \int W^{(1)}({\cal R}, {\cal R}^{\prime }, \tau )F({\cal R}^{\prime }, \tau )\,{\rm d}{\cal R}^{\prime } \,&=& 0, \quad \mbox{zero net wake at {\cal R} due to all the rings.} \end{eqnarray}

(30)

Formula (26a) for W⁽¹⁾ and the kinetic equation for F – either equations (27) and (28) or (29a)–(29c) – are the basic equations governing RR in the PRA. They should be solved with either the no-loss or the empty loss-cone boundary conditions on F and W⁽¹⁾ discussed in the previous section.

The PRA has a history in both plasma physics and stellar dynamics. The Landau (1936) equations for electrostatic plasmas is a PRA of the equations of Balescu (1960) and Lenard (1960). The general collisional theory for stellar systems, which was initiated by Gilbert (1968), was simplified by Polyachenko & Shukhman (1982) through a PRA applied to systems with integrable mean-field Hamiltonians. They computed the integral for the irreducible part of the two-particle correlation function and showed that only the resonant part contributes to the kinetic evolution. They also proved an H-theorem and derived a Fokker–Planck equation for spherical stellar systems. There has been some work on the applications of the kinetic equation (in either the Gilbert form or the Polyachenko–Shukhman form) but we do not review this literature. Even in the more tractable PRA form, the Polyachenko–Shukhman kinetic equation applies to inhomogeneous systems, and computing interactions between realistic stellar orbits can be complicated. Further simplification can be achieved through a local approximation (Chandrasekhar) by assuming that stars encounter each other on constant velocity orbits. Then one descends from the Polyachenko–Shukhman equation to the familiar description of classical two-body relaxation given in Binney & Tremaine (2008).

The PRA wake is given in equation (26a) as a formal time integral over the histories of the two ring orbits |${\cal R}(\tau ^{\prime })$| and |${\cal R}^{\prime }(\tau ^{\prime })$| in the quasi-statically deforming secular equilibrium, |$F({\cal R}, \tau )$|⁠. The formalism is general and does not require that the secular Hamiltonian, |$H({\cal R}, \tau )$|⁠, be integrable (at fixed τ). However, matters simplify when the secular Hamiltonian is integrable (at fixed τ), because we can proceed with the calculation of the time integral using techniques similar to Polyachenko & Shukhman (1982). In the rest of this paper we apply the PRA equations to the RR of razor-thin axisymmetric discs.

4 WAKES IN A RAZOR-THIN AXISYMMETRIC DISC

A razor-thin (flat) disc corresponds to the idealized case when all N rings are confined to the xy plane. The angular momentum of every ring points along |$\pm \hat{z}$|⁠, so a ring is specified by its semi-major axis, angular momentum and apsidal longitude. Ring space is three-dimensional: we write |${\cal R}= \lbrace I, L, g\rbrace$|⁠, where L stands for the angular momentum which can now take both positive and negative values ( − I ≤ L ≤ I ), and g is the longitude of the periapse.⁷ Earlier formulae apply with |${\rm d}{\cal R}\equiv {\rm d}I\,{\rm d}L\,{\rm d}g$|⁠. Ring space is topologically equivalent to |${R}^3$|⁠, with I the radial coordinate, |$\arccos {(L/I)}$| the colatitude, and g the azimuthal angle.

The secular theory applies to discs that are stable against fast axisymmetric instabilities.⁸ The general non-axisymmetric, time-dependent ring DF is F(I, L, g, τ) which is normalized as

\begin{equation} \int {\rm d}I\,{\rm d}L\,{\rm d}g\,\,F(I, L, g, \tau ) \,=\, 1\, \end{equation}

(31)

for the case of no loss of stars to the MBH. When the MBH is a sink of stars, the above normalization holds only at some initial reference time when the star cluster's mass is M. The secular Hamiltonian is

\begin{equation} H(I, L, g, \tau ) \,=\, \Phi (I, L, g, \tau ) \,+\, H^{\rm rel}(I, L) \,+\, \Phi ^{\rm tid}(I, L, g, \tau ), \end{equation}

(32)

where

\begin{equation} \Phi (I, L, g, \tau ) \,=\, \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,{\rm d}g^{\prime }\,\,\Psi (I, L, g, I^{\prime }, L^{\prime }, g^{\prime })\,F(I^{\prime }, L^{\prime }, g^{\prime }, \tau )\, \end{equation}

(33)

is the self-gravitational potential,

\begin{equation} H^{\rm rel}(I, L) \,=\, -B_1\frac{1}{I^3\vert L\vert } \,+\, B_{1.5} \frac{{\rm Sgn}(L)}{I^3L^2}\, \end{equation}

(34)

is the relativistic correction, and Φ^tid is an external tidal potential.

4.1 Collisionless limit

Ring orbits in the mean-field (or collisionless) limit are determined by the equations of motion:

\begin{equation} I \,=\, \sqrt{GM_\bullet a} \,=\, \mbox{constant},\quad \frac{{\rm d}L}{{\rm d}\tau } \,=\, -\,\frac{\mathrm{\partial} H}{\mathrm{\partial} g},\quad \frac{{\rm d}g}{{\rm d}\tau } \,=\, \frac{\mathrm{\partial} H}{\mathrm{\partial} L}. \end{equation}

(35)

The DF satisfies the ring CBE,

\begin{equation} \frac{{\rm d}F}{{\rm d}\tau } \,=\, \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,-\, \frac{\mathrm{\partial} H}{\mathrm{\partial} g}\frac{\mathrm{\partial} F}{\mathrm{\partial} L} \,+\, \frac{\mathrm{\partial} H}{\mathrm{\partial} L}\frac{\mathrm{\partial} F}{\mathrm{\partial} g} \,=\, \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,+\, \left[\,F,\,H\,\right] \,=\, 0, \end{equation}

(36)

where the Poisson Bracket is now two-dimensional, because it involves only the canonically conjugate pair (L, g). The probability for a ring to be in (I, I + dI) is a conserved quantity, i.e. the PDF in one-dimensional I-space, defined by

\begin{equation} P(I) \,=\, \int \, \text{d}L\,dg\,\, F(I,\, L,\, g, \,\tau ), \end{equation}

(37)

is independent of τ, as can be verified directly using the ring CBE equation (36).

Secular equilibria, both axisymmetric and non-axisymmetric, can be constructed by using the secular Jeans theorem. For an axisymmetric system any F(I, L) is a secular equilibrium because I and L are two isolating integrals of motion of the (integrable) secular Hamiltonian equation (32):

\begin{equation} H(I, L) \,=\, \Phi (I, L) \,+\, H^{\rm rel}(I, L) \,+\, \Phi ^{\rm ext}(I, L). \end{equation}

(38)

Here, we have allowed for an external axisymmetric potential Φ^ext(I, L); since such a field does not cause acceleration of the MBH, we have Φ^tid(I, L) = Φ^ext(I, L) as written above. Using equation (35) the orbit of a ring is given by

\begin{equation} I \,=\, \mbox{constant},\quad L \,=\, \mbox{constant},\quad \frac{{\rm d}g}{{\rm d}\tau } \,\equiv \, \Omega (I, L) \,=\, \frac{\mathrm{\partial} H}{\mathrm{\partial} L}. \end{equation}

(39)

Therefore, every ring has constant semi-major axis and eccentricity, with its apsidal longitude precessing at the constant rate Ω(I, L).

From equation (9) the LCBE determining the secular stability of an axisymmetric, secular equilibrium F(I, L) is

\begin{equation} \frac{\mathrm{\partial} F_1}{\mathrm{\partial} \tau } \,+\, \Omega \frac{\mathrm{\partial} F_1}{\mathrm{\partial} g} \,-\, \frac{\mathrm{\partial} F}{\mathrm{\partial} L}\frac{\mathrm{\partial} }{\mathrm{\partial} g}\int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,{\rm d}g^{\prime }\,\Psi (I, L, g, I^{\prime }, L^{\prime }, g^{\prime })\,F_1(I^{\prime }, L^{\prime }, g^{\prime }, \tau ) \,=\, 0. \end{equation}

(40)

The ‘bare’ ring–ring interaction potential, Ψ(I, L, g, I΄, L΄, g΄), is a complicated function of its arguments. But, even without the knowledge of its specific form, symmetry properties of Ψ (see below) were used in Paper I to obtain some general results on secular stability:

All axisymmetric perturbations of F(I, L) give rise to nearby axisymmetric equilibria.
When (∂F/∂L) does not change sign anywhere in |${\cal R}$|-space, the DF is neutrally stable to modes of all m (where m is the azimuthal wavenumber).

Other stability results for F(I, L) are: non-rotating discs described by DFs that are even functions of L with empty loss cones (F/L non-singular when L → 0 ) are likely to be unstable to m = 1 modes (Tremaine 2005); Mono-energetic counter-rotating discs dominated by nearly radial orbits are prone to non-axisymmetric loss-cone instabilities of all m (Polyachenko, Polyachenko & Shukhman 2007).

4.2 Wake function

The rest of this paper is devoted to the RR of axisymmetric discs described by DFs of the form F(I, L, τ), where the τ dependence is much slower than apse precessional periods. The self-gravitational potential Φ(I, L, τ), and hence the secular Hamiltonian H(I, L, τ) of equation (38), inherits the slow τ dependence of the DF. In the PRA kinetic equation (28) the Poisson Bracket [F , H − Φ/N] = 0 because all the quantities in it are independent of g. Then the PRA kinetic equation for discs in the BBGKY form is

\begin{equation} \frac{\mathrm{\partial} F}{\mathrm{\partial} \tau } \,=\, {\cal C}[F] \,=\, \frac{1}{N}\int \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,1)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right]\,{\rm d}{\cal R}^{\prime }. \end{equation}

(41)

This makes it clear that the RR time-scale – as measured by the slow time variable τ – is N times the apse precession period, 2π/Ω, where Ω(I, L, τ) = ∂H/∂L as in equation (39). The PRA wake function, W⁽¹⁾, of equation (26a) is a time integral over the orbits of the two rings, |${\cal R}(\tau ^{\prime })$| and |${\cal R}^{\prime }(\tau ^{\prime })$|⁠:

\begin{eqnarray} {\cal R}(\tau ^{\prime })\,:\quad I \,=\, \mbox{constant},\quad L \,=\, \mbox{constant},\quad g(\tau ^{\prime }) \,=\, g \,+\, \Omega \left\lbrace \tau ^{\prime } - \tau \right\rbrace ; \end{eqnarray}

(42a)

\begin{eqnarray} {\cal R}^{\prime }(\tau ^{\prime })\,:\quad I^{\prime } \,=\, \mbox{constant},\,\,\, L^{\prime } \,=\, \mbox{constant},\,\,\, g^{\prime }(\tau ^{\prime }) \,=\, g^{\prime } \,+\, \Omega ^{\prime }\left\lbrace \tau ^{\prime } - \tau \right\rbrace. \end{eqnarray}

(42b)

Here we use the shorthand, Ω = Ω(I, L, τ) and Ω΄ = Ω(I΄, L΄, τ). It should be noted that g(τ΄) and g΄(τ΄) refer to two different functions of τ΄, whereas g and g΄ are their values at τ΄ = τ. The integrand involves a Poisson Bracket [Φ^p, F] = [Ψ, F] = (∂F/∂L)(∂Ψ/∂g). Using equations (42a) and (42b) for the ring orbits in equation (26a), the PRA wake function is

\begin{equation*} W^{(1)}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, \int _{-\infty }^{\tau } \frac{\mathrm{\partial} F}{\mathrm{\partial} L}\frac{\mathrm{\partial} }{\mathrm{\partial} g}\Psi \left(I, L, g + \Omega (\tau ^{\prime } - \tau ), I^{\prime }, L^{\prime }, g^{\prime } + \Omega ^{\prime }(\tau ^{\prime } - \tau )\right)\,{\rm d}\tau ^{\prime }. \nonumber \end{equation*}

In the integrand the Ψ-term varies over ring precession times T_sec. We recall from the discussion following equations (26a) and (26b) that F varies over the much longer RR time T_res ≫ T_sec ; hence (∂F/∂L) can be pulled out of the integral. Dropping the superscript ‘1’, the PRA wake function is

\begin{eqnarray} W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) &=& \frac{\mathrm{\partial} F}{\mathrm{\partial} L}\int _{-\infty }^{\tau } \frac{\mathrm{\partial} }{\mathrm{\partial} g}\Psi \left(I, L, g + \Omega (\tau ^{\prime } - \tau ), I^{\prime }, L^{\prime }, g^{\prime } + \Omega ^{\prime }(\tau ^{\prime } - \tau )\right)\,{\rm d}\tau ^{\prime } \nonumber \\ &=& \frac{\mathrm{\partial} F}{\mathrm{\partial} L}\int _{-\infty }^{0} \frac{\mathrm{\partial} }{\mathrm{\partial} g}\Psi \left(I, L, g + \Omega \tau ^{\prime }, I^{\prime }, L^{\prime }, g^{\prime } + \Omega ^{\prime }\tau ^{\prime }\right)\,{\rm d}\tau ^{\prime }. \end{eqnarray}

(43)

The last equation makes it explicit that the τ dependence of the wake comes only from F, Ω and Ω΄, so W changes only over the RR time-scale T_res. We can also rewrite

\begin{eqnarray} W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, -\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\,\Delta L, \end{eqnarray}

(44a)

\begin{eqnarray} \mbox{where}\quad \Delta L \,=\, -\frac{\mathrm{\partial} }{\mathrm{\partial} g}\int _{-\infty }^{0}\Psi \left(I, L, g + \Omega \tau ^{\prime }, I^{\prime }, L^{\prime }, g^{\prime } + \Omega ^{\prime }\tau ^{\prime }\right)\,{\rm d}\tau ^{\prime }. \end{eqnarray}

(44b)

The physical meaning of these expressions is the following: ΔL/N is the net change in the specific angular momentum of ring |${\cal R}$|⁠, due to the torques exerted on it by ring |${\cal R}^{\prime }$|⁠. We recall from equation (13) that W/N is one of three components of the perturbation to the DF, resulting from the Rostoker–Gilbert gedanken experiment with rings. Then equation (44a) is just the standard PRA expression for the net change in the DF at any location |${\cal R}$|⁠. The wake is proportional to the angular momentum exchanged between the rings.

The form of equations (44a) and (44b) for the wake, together with some of the symmetry properties of Ψ, allows us to draw some general conclusions. We recall from Paper I the ‘bare’ ring–ring interaction potential, Ψ, has the following symmetries:

P1: Ψ is a real function which is even in both L and L΄ .
P2: The apsidal longitudes g and g΄ occur in Ψ only in the combination |g − g΄|.
P3: Ψ is independent of both g and g΄ when one of the two rings is circular (i.e. when L = ±I or L΄ = ±I΄ or both).
P4: Ψ is invariant under the interchange of the two rings. This can be achieved by any of the transformations: {I, L}↔{I΄, L΄} , or g↔g΄ or both. In explicit form we have: Ψ(I, L, g, I΄, L΄, g΄) = Ψ(I΄, L΄, g, I, L, g΄), Ψ(I, L, g, I΄, L΄, g΄) = Ψ(I, L, g΄, I΄, L΄, g), which implies that Ψ(I, L, g, I΄, L΄, g΄) = Ψ(I΄, L΄, g΄, I, L, g).

Let ΔL΄/N be the total change in the specific angular momentum of ring |${\cal R}^{\prime }$| due to ring |${\cal R}$|⁠, accrued over the past orbital histories of both rings. Since we consider equal-mass rings we expect ΔL΄ = −ΔL, by the conservation of angular momentum. Properties P2 and P4 can be used to verify this:

\begin{eqnarray} \Delta L^{\prime } \,=\, -\frac{\mathrm{\partial} }{\mathrm{\partial} g^{\prime }}\int _{-\infty }^{0}\Psi \left(I^{\prime }, L^{\prime }, g^{\prime } + \Omega ^{\prime }\tau ^{\prime }, I, L, g + \Omega \tau ^{\prime }\right)\,{\rm d}\tau ^{\prime }\,=\, \frac{\mathrm{\partial} }{\mathrm{\partial} g}\int _{-\infty }^{0}\Psi \left(I, L, g + \Omega \tau ^{\prime }, I^{\prime }, L^{\prime }, g^{\prime } + \Omega ^{\prime }\tau ^{\prime }\right)\,{\rm d}\tau ^{\prime } \,=\, -\Delta L. \end{eqnarray}

(45)

This also implies that we can write

\begin{equation} W({\cal R}^{\prime }\,\vert \,{\cal R}, \tau ) \,=\, -\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }}\,\Delta L^{\prime } \,=\, \frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }}\,\Delta L, \end{equation}

(46)

where we have used the shorthand F΄ = F(I΄, L΄, τ). When equations (44a) and (46) for the wakes are substituted into equation (27), the PRA two-ring correlation function is

\begin{equation} F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, \left\lbrace F\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }} \,-\, F^{\prime }\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\right\rbrace \Delta L \,=\, FF^{\prime }\left\lbrace \frac{\mathrm{\partial} \ln F^{\prime }}{\mathrm{\partial} L^{\prime }} \,-\, \frac{\mathrm{\partial} \ln F}{\mathrm{\partial} L}\right\rbrace \Delta L, \end{equation}

(47)

where we have suppressed the superscript ‘0’. The term in the second parentheses { } vanishes for a DF of the form

\begin{equation} F(I, L) \,=\, \frac{A(I)}{2\pi }\,\exp {\lbrace \gamma L\rbrace }, \end{equation}

(48)

where A(I) is an arbitrary function and γ is a real constant. This implies that |$F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }) = 0$| for all |${\cal R}$| and |${\cal R}^{\prime }$|⁠. Hence the DF of equation (48) is a stationary solution of the kinetic equation (41), and hence does not evolve through RR. We discuss this in more detail in Section 6.2, where this DF will emerge as an infinitely hot secular thermodynamic equilibrium of RR.

We verify here that the wake satisfies some essential properties.

The PRA wake function given by equations (44a) and (44b) must satisfy the two consistency conditions of equations (30). The first one states that there is no net mass in the wake of ring |${\cal R}^{\prime }$|⁠. Since ΔL = ∂{ }/∂g, we have
\begin{equation*} \int W({\cal R}, {\cal R}^{\prime }, \tau )\,{\rm d}{\cal R}\,=\, -\int {\rm d}I\,{\rm d}L \frac{\mathrm{\partial} F}{\mathrm{\partial} L}\int _0^{2\pi }{\rm d}g\,\Delta L \,=\, 0. \nonumber\end{equation*}
The second identity requires that the net wake at any location |${\cal R}$|⁠, due to all the rings, must vanish. Since ΔL = ∂{ }/∂g = −∂{ }/∂g΄ by property P2, we have
\begin{equation*} \int W({\cal R}, {\cal R}^{\prime }, \tau )F({\cal R}^{\prime }, \tau )\,{\rm d}{\cal R}^{\prime } \,=\, -\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,F(I^{\prime }, L^{\prime }, \tau ) \int _0^{2\pi }{\rm d}g^{\prime }\,\Delta L \,=\, 0. \end{equation*}
When one of the two rings is circular, P3 implies that Ψ is independent of g and g΄. Hence ΔL = −ΔL΄ = 0 and |$W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) = W({\cal R}^{\prime }\,\vert \,{\cal R}, \tau ) = 0$|⁠. In other words, a circular ring neither produces a wake nor feels the wake of another ring.

5 KINETIC EQUATION FOR A RAZOR-THIN AXISYMMETRIC DISC

5.1 Collision integral

In order to compute the collision integral it is necessary to do the time integral in equations (44a) and (44b) for the wake. It is convenient to begin with the Fourier-expansion of the ring–ring interaction potential function |$\Psi ({\cal R}, {\cal R}^{\prime })$| in g and g΄. By property P2 we know that Ψ depends only on the difference between g and g΄. We write

\begin{equation} \Psi (I, L, g, I^{\prime }, L^{\prime }, g^{\prime }) \,=\, \exp \lbrace \lambda \tau ^{\prime }\rbrace \sum _{m=-\infty }^{\infty }\,C_m(I, L, I^{\prime }, L^{\prime }) \exp {\left[{\rm i}m (g - g^{\prime })\right]}, \end{equation}

(49)

where λ > 0 is a small positive constant that enforces adiabatic switching by making the ring–ring interaction arbitrarily small in the distant past. The Fourier coefficients, C_m, are known functions of (I, L, I΄, L΄, m), but we do not need explicit expressions in this paper. As in Paper I we note that the symmetry properties of Ψ, given in P1– P4 of the previous section, translate to the following properties of the Fourier coefficients:

F1: the C_m are real functions that are even in L and L΄.
F2: the C_m are even in m: i.e. C_m = C_−m .
F3: when either L = ±I or L΄ = ±I΄ or both, then C_m = 0 for all m ≠ 0.
F4: C_m(I, L, I΄, L΄) = C_m(I΄, L΄, I, L).

Substituting equation (49) into equations (44a) and (44b) the wake is

\begin{eqnarray*} W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, \frac{\mathrm{\partial} F}{\mathrm{\partial} L} \sum _{m \ne 0}\,{\rm i}mC_m \exp {\left[{\rm i}m (g - g^{\prime })\right]} \int _{-\infty }^0 {\rm d}\tau ^{\prime }\exp {\lbrace \lambda + {\rm i}m(\Omega - \Omega ^{\prime })\rbrace \tau ^{\prime }} \,=\, \frac{\mathrm{\partial} F}{\mathrm{\partial} L} \sum _{m \ne 0}\,\frac{{\rm i}mC_m}{\lambda + {\rm i}m(\Omega - \Omega ^{\prime })}\exp {\left[{\rm i}m (g - g^{\prime })\right]}. \nonumber \end{eqnarray*}

Since m runs over all integer values except zero, in the limit λ → 0₊ , we have

\begin{equation*} \frac{{\rm i}m}{\lambda + {\rm i}m(\Omega - \Omega ^{\prime })} \,\longrightarrow \, \frac{1}{\Omega - \Omega ^{\prime }} \,+\, {\rm i}\pi \,{\rm Sgn}(m)\,\delta (\Omega - \Omega ^{\prime }), \end{equation*}

where δ(Ω − Ω΄) is a Dirac delta-function and Sgn(m) = ±1 is the sign of m. Then the wake can be written as the sum of two terms:

\begin{eqnarray} W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, W^{\rm r}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,+\, W^{\rm n}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )\,; \end{eqnarray}

(50a)

\begin{eqnarray} W^{\rm r}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, -2\pi \,\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\,\delta (\Omega - \Omega ^{\prime }) \sum _{m =1}^{\infty }\,C_m\sin \lbrace m(g - g^{\prime })\rbrace , \end{eqnarray}

(50b)

\begin{eqnarray} W^{\rm n}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,=\, \frac{\mathrm{\partial} F}{\mathrm{\partial} L}\left\lbrace \frac{\Psi ({\cal R}, {\cal R}^{\prime }) - C_0}{\Omega - \Omega ^{\prime }}\right\rbrace. \end{eqnarray}

(50c)

Here, W^r is the resonant part of the wake in which the δ(Ω − Ω΄) restricts contributions to pairs of rings that are in apsidal resonance. Wⁿ is the non-resonant part of the wake with no such restriction; it is proportional to (Ψ − C₀) = ∑_{m ≠ 0}C_mexp [im(g − g΄)] which is the part of Ψ that is fluctuating in (g − g΄).

Substituting equations (50a)–(50c) into equation (27) we can decompose |$F^{(2)}_{\rm irr}$| into resonant and non-resonant parts:

\begin{eqnarray} F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, F^{(2, {\rm r})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,+\, F^{(2, {\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,; \end{eqnarray}

(51a)

\begin{eqnarray} F^{(2, {\rm r})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, F^{\prime }\,W^{\rm r}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,+\, F\,W^{\rm r}({\cal R}^{\prime }\,\vert \,{\cal R}, \tau )\,=\, 2\pi \left\lbrace F\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }} \,-\, F^{\prime }\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\right\rbrace \delta (\Omega - \Omega ^{\prime }) \sum _{m =1}^{\infty }\,C_m\sin \lbrace m(g - g^{\prime })\rbrace , \end{eqnarray}

(51b)

\begin{eqnarray} F^{(2, {\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, F^{\prime }\,W^{\rm n}({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) \,+\, F\,W^{\rm n}({\cal R}^{\prime }\,\vert \,{\cal R}, \tau ) \,=\, -\left\lbrace F\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }} \,-\, F^{\prime }\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\right\rbrace \frac{\Psi ({\cal R}, {\cal R}^{\prime }) - C_0}{\Omega - \Omega ^{\prime }}. \end{eqnarray}

(51c)

The collision integral |${\cal C}[F]$| of equation (41) is

\begin{eqnarray} {\cal C}[F] \,&=& \frac{1}{N}\int {\rm d}I^{\prime }{\rm d}L^{\prime }\int _0^{2\pi }{\rm d}g^{\prime } \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right] = \frac{1}{N}\int {\rm d}I^{\prime }{\rm d}L^{\prime }\int _0^{2\pi }{\rm d}g^{\prime } \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,{\rm r})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right] \nonumber \\ &&+\, \frac{1}{N}\int {\rm d}I^{\prime }{\rm d}L^{\prime }\int _0^{2\pi }{\rm d}g^{\prime } \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,{\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right]. \end{eqnarray}

(52)

On the right side, the Poisson Bracket is to be taken over the |${\cal R}$| variables. The next step is to substitute equations (51b) and (51c) into equation (52) and evaluate the two integrals over g΄. The computations are straightforward and are given in the Appendix. Here we state the results: Since |$\left[\Psi , F^{(2,{\rm n})}_{\rm irr}\right] = \mathrm{\partial} \lbrace \,\rbrace /\mathrm{\partial} g^{\prime }\,$|⁠, the integral over g΄ vanishes and there is no contribution to |${\cal C}[F]$| from the non-resonant part of the wake. The resonant part gives

\begin{equation*} \int _0^{2\pi }{\rm d}g^{\prime } \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,{\rm r})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right] \,=\, -4\pi ^2\frac{\mathrm{\partial} }{\mathrm{\partial} L}\left\lbrace \delta (\Omega - \Omega ^{\prime }) \left\lbrace F\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }} - F^{\prime }\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\right\rbrace \sum _{m =1}^{\infty }\,mC_m^2\right\rbrace. \end{equation*}

Then

\begin{equation} {\cal C}[F] \,=\, -\frac{\mathrm{\partial} }{\mathrm{\partial} L} \left\lbrace \frac{4\pi ^2}{N} \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime }) \left\lbrace F\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }} - F^{\prime }\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\right\rbrace \sum _{m =1}^{\infty }\,mC_m^2\right\rbrace \, \end{equation}

(53)

gives the collision integral explicitly as a functional of the DF.

5.2 Kinetic equation

Since all the (g, g΄) dependences have been eliminated from the collision integral of equation (53), we can now describe the RR of axisymmetric Keplerian discs by the new DF, f(I, L, τ) = 2πF(I, L, τ), which is a PDF in (I, L) space, normalized as⁹

\begin{equation} \int {\rm d}I\,{\rm d}L\, f(I, L, \tau ) \,=\, 1. \end{equation}

(54)

The self-gravitational potential is given in terms of f as

\begin{equation} \Phi (I, L, \tau ) \,=\, \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\, C_0(I, L, I^{\prime }, L^{\prime })f(I^{\prime }, L^{\prime }, \tau ). \end{equation}

(55)

The secular Hamiltonian is H(I, L, τ) = Φ(I, L, τ) + H^rel(I, L) + Φ^ext(I, L) with the apse precession rate Ω(I, L, τ) = ∂H/∂L . In the collision integral of equation (53) the net strength of the gravitational encounter between two rings is given by the sum over m. It is convenient to define the interaction coefficient (or kernel),

\begin{equation} K(I, L, I^{\prime }, L^{\prime }) \,=\, 2\pi \sum _{m =1}^{\infty }\,mC_m^2(I, L, I^{\prime }, L^{\prime }), \end{equation}

(56)

which will be regarded as a known function of its four arguments. K inherits some general properties from those of the C_m, given under F1– F4 in Section 5.1. These are:

K1: K is a non-negative function.
K2: K(I, L, I΄, L΄) is an even function of L and L΄.
K3: when either L = ±I or L΄ = ±I΄ or both, then K = 0.
K4: K(I, L, I΄, L΄) = K(I΄, L΄, I, L) is invariant under the interchange (I, L) ↔ (I΄, L΄).

From equations (41), (53) and (56) the kinetic equation for axisymmetric discs in ‘conservation’ form is

\begin{equation} \frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } \,+\, \frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,=\, 0, \end{equation}

(57)

where

\begin{equation} J(I, L, \tau ) \,=\, \frac{1}{N} \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K(I, L, I^{\prime }, L^{\prime }) \left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace \, \end{equation}

(58)

is the probability current density of RR in (I, L)-space. It is directed only along the L-direction because RR does not cause any change in the I of any ring. J(I, L, τ) is a functional of f, given as an integral over two-dimensional (I΄, L΄)-space which picks up contributions only from the one-dimensional apsidal-resonance curve defined by Ω(I΄, L΄, τ) = Ω(I, L, τ). Two useful general properties of J – that are valid for any non-singular DF f(I, L, τ) and either the ‘no-loss’ or ‘empty loss-cone’ boundary conditions – are the following:

(a) The maximum and minimum allowed values of L are I and −I, corresponding to prograde and retrograde circular orbits. Since these represent the upper and lower boundaries of the allowed regions of (I, L) space (see Figs 1 –3), we expect the current to vanish on the boundaries L = ±I – otherwise there would be a ‘leakage’ of rings into the forbidden regions, I ≤ |L| . This can be verified immediately: by property A3 we have K(I, ±I, I΄, L΄) = 0, so

\begin{equation} J(I, \pm I, \tau ) \,=\, 0. \end{equation}

(59)

(b) Let ξ(Ω) be an arbitrary non-singular function of Ω. Then the integral of ξ(Ω) J(I, L, τ) over (I, L)-space vanishes. To see this we note that

\begin{eqnarray*} \int {\rm d}I\,{\rm d}L\,\xi (\Omega )\,J &=& \frac{1}{N} \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K(I, L, I^{\prime }, L^{\prime })\,\xi (\Omega )\left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace \,\nonumber \\ &=& -\frac{1}{N} \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K(I, L, I^{\prime }, L^{\prime })\,\xi (\Omega ^{\prime })\left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace , \nonumber \end{eqnarray*}

because: under the interchange of integration variables, (I, L) ↔ (I΄, L΄) on the right side, both δ(Ω − Ω΄) and K are invariant whereas {f∂f΄/∂L΄ − f΄∂f/∂L} changes sign. Adding and dividing by two, we obtain the following general identity:

\begin{eqnarray} \int {\rm d}I\,{\rm d}L\,\xi (\Omega )\,J \,=\, \frac{1}{2N} \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, \left\lbrace \xi (\Omega ) - \xi (\Omega ^{\prime })\right\rbrace \left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace \,=\, 0, \end{eqnarray}

(60)

because of the factor δ(Ω − Ω΄){ξ(Ω) − ξ(Ω΄)} in the integrand. This is a general symmetry property of the current density J. Two particular cases we will use correspond to ξ(Ω) = 1 and ξ(Ω) = Ω; these give

\begin{equation} \int {\rm d}I\,{\rm d}L\,J \,=\, 0,\quad \int {\rm d}I\,{\rm d}L\,\Omega \,J \,=\, 0. \end{equation}

(61)

Equations (61) will be used in Section 6.1 to prove conservation of total angular momentum and energy of the disc.

Kinetic equation in the Fokker–Planck form: the quantities f and ∂f/∂L can be pulled out of the integrals in equation (58), and J written as a sum of advective and diffusive contributions:

\begin{eqnarray} J(I, L, \tau ) \,=\, J^{\rm adv}(I, L, \tau ) \,+\, J^{\rm dif}(I, L, \tau )\,; \end{eqnarray}

(62a)

\begin{eqnarray} J^{\rm adv} \,=\, Vf, \quad \mbox{where}\quad \,V(I, L, \tau ) \,=\, \frac{1}{N} \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, \frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }}\,, \end{eqnarray}

(62b)

\begin{eqnarray} J^{\rm dif} \,=\, -D\frac{\mathrm{\partial} f}{\mathrm{\partial} L}, \quad \mbox{where}\quad \, D(I, L, \tau ) \,=\, \frac{1}{N} \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, f^{\prime }\,. \end{eqnarray}

(62c)

The current J^adv represents an advective RR flux in L-space. It can be thought of as arising from the drag-torque exerted by the system on a ring at (I, L). Then V can be interpreted as a local velocity in L-space with which the DF is being advected. J^adv can be positive or negative because V can take any real value. The current J^dif expresses Fick's law in L-space, with positive diffusion because D is always non-negative. Both V(I, L, τ) and D(I, L, τ) are transport coefficients which are functionals of the DF. We can arrive at a physical interpretation of the kinetic equation by casting it in the Fokker–Planck form (Lifshitz & Pitaevskii 1981). Using equations (62a)–(62c) the kinetic equation (57) can be rewritten as

\begin{equation} \frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } \,=\, -\frac{\mathrm{\partial} }{\mathrm{\partial} L}\left(U f\right) \,+\, \frac{\mathrm{\partial} ^2 }{\mathrm{\partial} L^2}\left(D f\right), \end{equation}

(63)

where

\begin{equation} U(I, L, \tau ) \,=\, V \,+\, \frac{\mathrm{\partial} D}{\mathrm{\partial} L}\, \end{equation}

(64)

is given in terms of V and D defined in equations (62b) and (62c). Equation (63) is in the Fokker–Planck form, and U and D are diffusion coefficients that are expressible in terms of the average characteristics, per unit time, of the random changes in the specific angular momentum of a ring at (I, L):

\begin{equation} U(I, L, \tau ) \,=\, \frac{\left\langle \delta L\right\rangle }{\delta \tau }, \quad D(I, L, \tau ) \,=\, \frac{\left\langle \left(\delta L\right)^2\right\rangle }{2\,\delta \tau }. \end{equation}

(65)

The ring at (I, L) is viewed as a test ring that is embedded in a sea of other rings which exert stochastic torques on it. U is the mean change in L per unit time, and D equals one-half of the mean-squared change in L per unit time.

We now make an estimate of the RR time-scale implied by the Fokker–Planck equation (63). From equation (62c) we estimate that D ∼ K/NΩ , where we have used the fact that the DF is normalized as given in equation (54). From equation (56) for the kernel K in terms of the Fourier coefficients C_m, we have K ∼ Ψ² . Using the definition of Ψ given in equation (3b), K ∼ (GM_•/a)² where a is a typical semi-major axis. The precession frequency Ω ∼ (GM_•/a³)^1/2 because time has been measured in terms of the slow time variable τ = εt. Then the diffusion coefficient D ∼ (GM_•)^3/2/Na^1/2 . The RR time-scale, T_res , can be estimated by requiring that over times, δτ ∼ εT_res, the typical value of 〈(δL)²〉 ∼ L² ∼ GM_•a . Then equation (65) gives T_res ∼ 〈(δL)²〉/εD ∼ (N/ε)(a³/GM_•)^1/2 ∼ N(T_kep/ε) = NT_sec , as expected.

5.3 Mass, angular momentum and energy

We define the mass, angular momentum and energy of the stellar disc below. In Section 6 we will prove that these quantities are, as expected, conserved when there is no loss of stars to the MBH; hence the normalization of the DF f, given in equation (54), is valid for all τ. However, when stars are lost to the MBH this normalization is valid only at some initial time τ = τ₀. At later times τ > τ₀, we will have ∫ dI dL f(I, L, τ) < 1 , and the mass, angular momentum and energy of the stellar disc will be different; this fuelling of the black hole is treated in Section 7.

The probability that a ring lies in (I, I + dI) is

\begin{equation} P(I, \tau ) \,=\, \int {\rm d}L\, f(I, L, \tau ). \end{equation}

(66)

The disc mass is

\begin{equation} {\cal M}(\tau ) \,=\, M\int {\rm d}I\,P(I, \tau ) \,=\, M\int {\rm d}I\,{\rm d}L\, f(I, L, \tau ). \end{equation}

(67)

The disc angular momentum is

\begin{equation} {\cal L}(\tau ) \,=\, M\int {\rm d}I\,{\rm d}L\, L\,f(I, L, \tau ). \end{equation}

(68)

The disc energy is

\begin{eqnarray} {\cal E}(\tau ) \,=\, M\int {\rm d}I\,{\rm d}L\,E_{\rm k}(I)\,f(I, L, \tau ) \,+\, \frac{\varepsilon M}{2}\int {\rm d}I\,{\rm d}L\, \Phi (I, L)\,f(I, L, \tau ) \,+\, \varepsilon M\int {\rm d}I\,{\rm d}L\, \left\lbrace H^{\rm rel}(I, L) \,+\, \Phi ^{\rm ext}(I, L)\right\rbrace f(I, L, \tau ), \end{eqnarray}

(69)

where E_k(I) = −1/2(GM_•/I)² is the Kepler energy of a ring of unit mass. The second term is the mutual self-gravitational energy of all the rings, and the third term accounts for general relativistic corrections to 1.5 PN order and external gravitational sources. Using equation (55) the self-gravitational potential Φ can be expressed in terms of f, and the disc energy written as

\begin{eqnarray} {\cal E}(\tau ) = M\int {\rm d}I\,{\rm d}L\, \left\lbrace E_{\rm k}(I) \,+\,\varepsilon H^{\rm rel}(I, L) \,+\, \varepsilon \Phi ^{\rm ext}(I, L)\right\rbrace f(I, L, \tau )+ \frac{\varepsilon M}{2}\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\, C_0(I, L, I^{\prime }, L^{\prime })\,f(I, L, \tau )\,f(I^{\prime }, L^{\prime }, \tau ). \end{eqnarray}

(70)

We note that P(I, τ), |$\,{\cal M}(\tau )$| and |$\,{\cal L}(\tau )$| are linear functionals of f, whereas |$\,{\cal E}(\tau )$| is a quadratic functional of f.

6 PHYSICAL KINETICS OF A RAZOR-THIN LOSSLESS DISC

Here, we study the general properties of a stellar disc when there is no loss of stars to the MBH. The domain of the DF f(I, L, τ) is the maximum allowed: −I ≤ L ≤ I and I ≥ 0 , shown as the shaded region in Fig. 1 lying between the boundaries L = ±I that mark the prograde and retrograde circular orbits. In the shaded region f can take any non-negative value, subject to the normalization of equation (54). Below we verify conservation of disc mass, angular momentum and energy. We also prove an H-theorem on the increase of the Boltzmann entropy and discuss secular thermal equilibrium.

The region lying between the two diagonal lines (L = ±I) is the domain of the DF, f, when there is no loss of stars to the black hole. The upper and lower boundaries correspond to prograde and retrograde circular orbits, respectively; the current J vanishes on these boundaries.

Figure 1.

The region lying between the two diagonal lines (L = ±I) is the domain of the DF, f, when there is no loss of stars to the black hole. The upper and lower boundaries correspond to prograde and retrograde circular orbits, respectively; the current J vanishes on these boundaries.

Open in new tab Download slide

6.1 Conserved quantities

Since the stellar disc is lossless and the external sources of gravity are time-independent, we expect that the quantities P, |${\cal M}$|⁠, |${\cal L}$| and |${\cal E}$| are all constant in time. It is useful to verify that the kinetic equations (57)–(58) do indeed conserve these physical quantities:

\begin{equation} \frac{\mathrm{\partial} P}{\mathrm{\partial} \tau } \,=\, \int {\rm d}L\,\frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } \,=\, - \int _{-I}^{I}{\rm d}L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,=\, - \left\lbrace J(I, I, \tau ) \,-\, J(I, -I, \tau )\right\rbrace \,=\, 0, \end{equation}

(71)

because J(I, ±I, τ) = 0 from equation (59). The probability distribution of the rings in I space is frozen in time, because RR does not change the semi-major axis of any ring. Then it follows that the total disc mass is conserved:

\begin{equation} \frac{{\rm d}{\cal M}}{{\rm d}\tau } \,=\, M\int {\rm d}I\,\frac{\mathrm{\partial} P}{\mathrm{\partial} \tau } \,=\, 0\, \end{equation}

(72)

also verifies the consistency of the normalization equation (54) used for the DF. Similarly, using equation (59) and the first of equations (61), we have

\begin{equation*} \frac{{\rm d}{\cal L}}{{\rm d}\tau } \,=\, M\int {\rm d}I\,{\rm d}L\, L\,\frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } \,=\, - M\int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\,L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,=\, M\int {\rm d}I\,{\rm d}L\,J \,=\, 0. \end{equation*}

Proving energy conservation requires some manipulations because the self-gravity term in equation (70) has an integrand which is quadratic in the DF:

\begin{eqnarray*} \frac{{\rm d}{\cal E}}{{\rm d}\tau } &=& -M\int _0^\infty {\rm d}I\,E_{\rm k}(I)\int _{-I}^{I}{\rm d}L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,-\, \varepsilon M\int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\, \left\lbrace H^{\rm rel}(I, L) \,+\, \Phi ^{\rm ext}(I, L)\right\rbrace \,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \nonumber \\ &&-\, \frac{\varepsilon M}{2}\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\, C_0(I, L, I^{\prime }, L^{\prime })\, \left\lbrace f\,\frac{\mathrm{\partial} J^{\prime }}{\mathrm{\partial} L^{\prime }} \,+\, f^{\prime }\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L}\right\rbrace. \nonumber \end{eqnarray*}

The first term on the right side vanishes because J vanishes at the upper/lower boundaries L = ±I. Similarly the second term can be integrated by parts. The last term can be simplified because C₀(I, L, I΄, L΄) is symmetric under interchange of primed and unprimed variables:

\begin{eqnarray} \frac{1}{2}\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,C_0(I, L, I^{\prime }, L^{\prime }) \left\lbrace f\frac{\mathrm{\partial} J^{\prime }}{\mathrm{\partial} L^{\prime }} + f^{\prime }\frac{\mathrm{\partial} J}{\mathrm{\partial} L}\right\rbrace = \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,C_0(I, L, I^{\prime }, L^{\prime })f^{\prime }\frac{\mathrm{\partial} J}{\mathrm{\partial} L}= \int {\rm d}I\,{\rm d}L\,\Phi (I, L, \tau )\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L}, \end{eqnarray}

(73)

where we have used equation (55) relating Φ and the DF. Hence, we have

\begin{eqnarray} \frac{{\rm d}{\cal E}}{{\rm d}\tau } &=& -\varepsilon M\int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\, \left\lbrace \Phi (I, L, \tau ) \,+\, H^{\rm rel}(I, L) \,+\, \Phi ^{\rm ext}(I, L)\right\rbrace \,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} = -\varepsilon M\int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\, H(I, L, \tau )\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L}\nonumber \\ &=&\varepsilon M\int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\, \frac{\mathrm{\partial} H}{\mathrm{\partial} L}\,J = \varepsilon M \int {\rm d}I\,{\rm d}L\,\Omega \,J \,=\, 0, \end{eqnarray}

(74)

where the second identity of equation (61) has been used.

6.2 H-theorem

An H-function |$\,{\cal H}[f]\,$| is a functional of f, defined by

\begin{equation} {\cal H}[f] \,=\, -\int {\rm d}I\,{\rm d}L\,{{C}}(f), \end{equation}

(75)

where |${{C}}(f)$| is a convex function of f. We assume that |${{C}}(f)$| is at least twice-differentiable, so convexity implies that |$\,{\rm d}^2{{C}}/{\rm d}f^2 \ge 0\,$|⁠.

\begin{eqnarray*} \frac{{\rm d}{\cal H}}{{\rm d}\tau }= -\int {\rm d}I\,{\rm d}L\,\frac{{\rm d}{{C}}}{{\rm d}f}\, \frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } &=& \int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\,\frac{{\rm d}{{C}}}{{\rm d}f}\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L}= -\int _0^\infty {\rm d}I\int _{-I}^{I}{\rm d}L\,J\,\frac{\mathrm{\partial} }{\mathrm{\partial} L}\left\lbrace \frac{{\rm d}{{C}}}{{\rm d}f}\right\rbrace \,=\, -\int {\rm d}I\,{\rm d}L\,\frac{{\rm d}^2{{C}}}{{\rm d}f^2}\,\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\,J \nonumber \\ &=& -\frac{1}{N} \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, \frac{{\rm d}^2{{C}}}{{\rm d}f^2}\,\frac{\mathrm{\partial} f}{\mathrm{\partial} L} \left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace. \nonumber \end{eqnarray*}

Both K and the δ-function are invariant under the interchange of integration variables (I, L) ↔ (I΄, L΄), whereas the term in { } changes sign. Hence the integrand can be replaced by

\begin{equation*} -\delta (\Omega - \Omega ^{\prime })\,K\, \frac{{\rm d}^2{{C}}^{\prime }}{{\rm d}f^{\prime 2}}\,\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} \left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace , \end{equation*}

where |${{C}}^{\prime } = {{C}}(f^{\prime })$|⁠. Adding these two forms of the integrand and dividing by two, we obtain

\begin{eqnarray} \frac{{\rm d}{\cal H}}{{\rm d}\tau } &=& \frac{1}{2N} \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, \left\lbrace \frac{{\rm d}^2{{C}}^{\prime }}{{\rm d}f^{\prime 2}}\,\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - \frac{{\rm d}^2{{C}}}{{\rm d}f^2}\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace \left\lbrace f\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - f^{\prime }\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace \nonumber \\ &=& \frac{1}{2N} \int {\rm d}I{\rm d}L{\rm d}I^{\prime }{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\,ff^{\prime } \left\lbrace \frac{{\rm d}^2{{C}}^{\prime }}{{\rm d}f^{\prime 2}}\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - \frac{{\rm d}^2{{C}}}{{\rm d}f^2}\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace \left\lbrace \frac{1}{f^{\prime }}\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - \frac{1}{f}\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace , \end{eqnarray}

(76)

in a symmetric form. For a general convex function |${{C}}$|(f) the right side can be either positive or negative. Therefore, an H-function of general form is not a very useful probe of RR evolution. But for the convex function with |${\rm d}^2{{C}}/{\rm d}f^2 = k/f\,$|⁠, where k is a positive constant, the right side of equation (76),

\begin{equation*} \frac{k}{2N} \int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\,ff^{\prime } \left\lbrace \frac{\mathrm{\partial} \ln f^{\prime }}{\mathrm{\partial} L^{\prime }} - \frac{\mathrm{\partial} \ln f}{\mathrm{\partial} L}\right\rbrace ^{\!2}, \end{equation*}

is always non-negative. This happens for convex functions of the form |${{C}}(f) = k f\ln f + k^{\prime }f + k^{\prime \prime }$|⁠, where k΄ and k΄ are arbitrary constants. We choose k = 1 and k΄ = 0 and k΄ = 0 and define the (Boltzmann) entropy, S[f], as the H-function corresponding to |${{C}}(f) = f\ln f\,$|⁠. Then,

\begin{eqnarray} S[f] \,=\, -\int {\rm d}I\,{\rm d}L\,f\ln f, \end{eqnarray}

(77)

\begin{eqnarray} \frac{{\rm d}S}{{\rm d}\tau } \,=\, \frac{1}{2N}\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\,ff^{\prime }\left\lbrace \frac{\mathrm{\partial} \ln f^{\prime }}{\mathrm{\partial} L^{\prime }} - \frac{\mathrm{\partial} \ln f}{\mathrm{\partial} L}\right\rbrace ^{\!2} \,\,\ge \, 0. \end{eqnarray}

(78)

The entropy is an increasing function of time, with the equality occurring at RR equilibrium.

6.3 Secular thermal equilibrium

Let us consider DFs that extremize the Boltzmann entropy S, when the following three conserved quantities are held fixed: the probability DF P(I), the disc energy |${\cal E}$|⁠, and the disc angular momentum |${\cal L}$|⁠. Using Lagrange multipliers B(I), −β΄ and γ΄, the extremization problem is

\begin{equation} \delta S \,+\, \int {\rm d}I\,B(I)\, \delta P(I) \,-\, \beta ^{\prime }\,\delta {\cal E}\,+\, \gamma ^{\prime }\,\delta {\cal L}\,=\, 0, \end{equation}

(79)

where δS, δP(I), |$\delta {\cal E}$| and |$\delta {\cal L}$| are infinitesimal changes corresponding to an infinitesimal change, δf, in the DF. When the expressions in equations (77), (66), (70) and (68) are used for |$(S, \,P, \,{\cal E}, \,{\cal L})$|⁠, we get

\begin{equation} \int {\rm d}I\,{\rm d}L\, \left\lbrace -\ln f \,-\, 1 \,+\, B(I) \,-\, \beta ^{\prime }M\!\left(E_{\rm k} + \varepsilon H\right) \,+\, \gamma ^{\prime }ML\right\rbrace \,\delta f \,=\, 0, \end{equation}

(80)

for arbitrary δf. This implies that the DFs extremizing S must have the general form

\begin{equation} f_{\rm eq}(I, L) \,=\, A(I) \exp {\left\lbrace -\beta H(I, L) + \gamma L\right\rbrace }, \end{equation}

(81)

where A(I) is an arbitrary positive function, and β and γ are constant parameters.¹⁰ The equilibrium DF, f_eq, is given in terms of the secular Hamiltonian, H(I, L) = Φ(I, L) + H^rel(I, L) + Φ^ext(I, L) . Here H^rel and Φ^ext are known functions of (I, L), but the self-gravitational potential, Φ(I, L), is determined by the DF itself through the integral relation of equation (55). In order to determine f_eq explicitly as a function of (I, L), one must solve the integral equation

\begin{eqnarray} \Phi (I, L) = \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\, C_0(I, L, I^{\prime }, L^{\prime }) A(I^{\prime }) \exp {\left\lbrace -\beta \left[\Phi (I^{\prime }, L^{\prime }) + H^{\rm rel}(I^{\prime }, L^{\prime }) + \Phi ^{\rm ext}(I^{\prime }, L^{\prime }) \right] + \gamma L^{\prime }\right\rbrace }, \end{eqnarray}

(82)

for Φ(I, L) and substitute this into equation (81) to get f_eq(I, L). Once we have a normalized f_eq the total energy and angular momentum of the disc, |${\cal E}$| and |${\cal L}$|⁠, can be computed from equations (70) and (68). These will be given as functions, |${\cal E}(\beta , \gamma )$| and |${\cal L}(\beta , \gamma )$|⁠, of the constants β and γ. To the best of our knowledge secular equilibria have been explored (Touma & Tremaine 2014) only in the limit when all the rings have the same semi-major axis.

The DF of equation (81) is evidently a secular collisionless equilibrium because it is a function of I and L, which are two secular isolating integrals of motion in the axisymmetric disc geometry. It is also a secular thermal equilibrium, because f_eq(I, L) is a time-independent solution of the kinetic equation (57). In order to verify this we note that ∂f_eq/∂L = ( − βΩ + γ)f_eq , so that |$\left\lbrace f_{\rm eq}\mathrm{\partial} f^{\prime }_{\rm eq}/\mathrm{\partial} L^{\prime } - f^{\prime }_{\rm eq}\mathrm{\partial} f_{\rm eq}/\mathrm{\partial} L\right\rbrace = \beta f_{\rm eq}\,f^{\prime }_{\rm eq}\left(\Omega - \Omega ^{\prime }\right)$|⁠. Then equation (58) implies that the equilibrium RR current density vanishes everywhere:

\begin{equation} J_{\rm eq}(I, L) \,=\, \frac{\beta }{N}\int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\left(\Omega - \Omega ^{\prime }\right)f_{\rm eq}\,f^{\prime }_{\rm eq} \,=\, 0, \end{equation}

(83)

because of the factor (Ω − Ω΄) δ(Ω − Ω΄) in the integrand. This implies that the collision integral, |$\,{\cal C}[f_{\rm eq}] = -\mathrm{\partial} J_{\rm eq}/\mathrm{\partial} L\,$|⁠, also vanishes everywhere in (I, L)-space. Therefore, the f_eq of equation (81) are secular thermal equilibria where the parameter β is an inverse temperature. We can now calculate explicitly the two-ring correlation at equilibrium by substituting equation (81) for the equilibrium DF into equations (51a)–(51c) to get the irreducible parts. Then, the full equilibrium correlation function between two rings |${\cal R}=(I, L, g)$| and |${\cal R}^{\prime }=(I^{\prime }, L^{\prime }, g^{\prime })$| is

\begin{eqnarray} F^{(2)}_{\rm eq}({\cal R}, {\cal R}^{\prime }) = \frac{1}{4\pi ^2}f_{\rm eq}(I, L)\,f_{\rm eq}(I^{\prime }, L^{\prime }) \left\lbrace 1 +\frac{2\pi \beta }{N} (\Omega - \Omega ^{\prime })\delta (\Omega - \Omega ^{\prime }) \sum _{m =1}^{\infty }\,C_m\sin \lbrace m(g - g^{\prime })\rbrace - \frac{2\beta }{N} \sum _{m =1}^{\infty }\,C_m\cos \lbrace m(g - g^{\prime })\rbrace \right\rbrace , \end{eqnarray}

(84)

where the two terms proportional to β/N are the irreducible resonant and non-resonant contributions.

Since f_eq(I, L) does not evolve either dynamically or thermally all H-functions are also time-independent, as can be verified using equation (76):

\begin{eqnarray} \left.\frac{{\rm d}{\cal H}}{{\rm d}\tau }\right|_{\rm eq} \,=\, \frac{\beta }{2N}\!\!\int {\rm d}I{\rm d}L{\rm d}I^{\prime }{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\left(\Omega - \Omega ^{\prime }\right) \left\lbrace \frac{{\rm d}^2{{C}}^{\prime }}{{\rm d}f^{\prime 2}}\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} L^{\prime }} - \frac{{\rm d}^2{{C}}}{{\rm d}f^2}\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right\rbrace _{\rm eq} \!\!\!f_{\rm eq}\,f^{\prime }_{\rm eq} = 0, \end{eqnarray}

(85)

because of the factor (Ω − Ω΄) δ(Ω − Ω΄) in the integrand, and the fact that the term in { } with the convex functions |${{C}}$| and |${{C}}^{\prime }$| is well-defined and non-singular. During the approach to equilibrium, when f ≠ f_eq, a general H-function need not evolve monotonically in time. Hence it could have been either greater or lesser than |${\cal H}[f_{\rm eq}]$|⁠, its value at equilibrium. So we cannot immediately conclude whether, for DFs close to f_eq, the H-function is larger or smaller than |${\cal H}[f_{\rm eq}]$|⁠. But the situation is different for the Boltzmann entropy because, by equation (78), S is a non-decreasing function. We verify

\begin{equation} \frac{{\rm d}S}{{\rm d}\tau } \,=\, \frac{\beta ^2}{2N}\!\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\left(\Omega - \Omega ^{\prime }\right)^2 f_{\rm eq}\,f^{\prime }_{\rm eq} \,\,=\, 0, \end{equation}

(86)

because of the factor (Ω − Ω΄)² δ(Ω − Ω΄) in the integrand. Since S could not have been larger in the past, the f_eq(I, L) are also maximum entropy equilibria, for given |$(P(I), \,{\cal E}, \,{\cal L})$|⁠, when compared with nearby axisymmetric secular collisionless equilibria f(I, L). These axisymmetric f_eq(I, L), however, need not be maximum entropy states when the variations in the DF are allowed to be non-axisymmetric.

The high temperature limit β → 0 : In this limit the equilibrium DF

\begin{equation} f_{\rm eq} \,\rightarrow \, A(I)\exp {\lbrace \gamma L\rbrace }, \end{equation}

(87)

which is identical to the equilibrium DF of equation (48). This form for a secular thermal equilibrium was first discussed in Rauch & Tremaine (1996). The DF of equation (87) has two noteworthy properties:

The irreducible parts of the two-ring correlation function of equation (84) are proportional to β, and they vanish in the high temperature limit, which is consistent with the discussion following equation (48). Then the two-ring correlation function of equation (84) approaches its mean-field value, (4π²)⁻¹f_eq(I, L)f_eq(I΄, L΄).
By a stability criterion proved in Paper I (mentioned in Section 4.1) these very hot DFs are dynamically stable to all axisymmetric and non-axisymmetric secular perturbations, because (∂f_eq/∂L) has the same sign as the constant γ, for all I and L.

Thermal stability of f_eq to axisymmetric perturbations can be studied by linearizing the kinetic equation (57) about f_eq. The thermal stability problem for non-axisymmetric perturbations is an important problem that needs to be formulated by going back to the PRA kinetic equation (28).

7 BLACK HOLE FUELLING RATES

A star that comes close enough to the MBH will be either torn apart by its tidal field or swallowed whole by it. During RR the semi-major axis of every star is conserved whereas its eccentricity can change because of stochastic angular momentum exchanges with other stars. When a star loses enough angular momentum, for its pericentre to come closer than some minimum distance from the MBH, it is counted as lost to the MBH. In this section we derive formulae for the rates at which RR feeds mass, angular momentum and energy to the MBH. Since the stellar system loses mass, we can only require that at some initial time τ₀,

\begin{equation} \int {\rm d}I\,{\rm d}L\, f(I, L, \tau _0) \,=\, 1,\quad \mbox{initial normalization.} \end{equation}

(88)

At later times τ > τ₀, we will have ∫ dI dL f(I, L, τ) ≤ 1 . As discussed in Section 2.2, the loss of stars at small distances from the MBH can be simply modelled by an absorbing barrier: we assume that a star is lost to the MBH when its pericentre distance is smaller than some fixed value r_lc, the ‘loss-cone radius’. When the loss cone is empty, stars belonging to the cluster must necessarily have pericentre radii a(1 − e) larger than r_lc. Hence the domain of f(I, L, τ) is restricted to regions of (I, L)-space in which I and |L| are large enough:

\begin{eqnarray} I_{\rm lc} \,<\, I \quad \mbox{where}\quad I_{\rm lc} \,=\, \sqrt{GM_\bullet r_{\rm lc}};\qquad L_{\rm lc}(I) \,<\, |L| \,\le \, I, \quad \mbox{where}\quad L_{\rm lc}(I) \,=\, I_{\rm lc}\left[\,2 \,-\, \left(\frac{I_{\rm lc}}{I}\right)^2\,\right]^{1/2}. \end{eqnarray}

(89)

This domain consists of two disjoint wedge-like regions of (I, L) space, shown shaded in Fig. 2. The dashed lines marked ‘±lcb’ are the loss-cone boundaries for prograde (L > 0) and retrograde (L < 0) orbits, respectively, and are given by

\begin{equation} \mbox{{ $+$lcb:}}\quad L \,=\, L_{\rm lc}(I)\quad \mbox{for $I \,\ge \, I_{\rm lc}\,$;}\quad \mbox{{ $-$lcb:}}\quad L \,=\, -L_{\rm lc}(I)\quad \mbox{for $I \,\ge \, I_{\rm lc}\,$.} \end{equation}

(90)

We impose the empty loss-cone boundary condition on f by requiring that it vanishes on both ±lcb:

\begin{equation} f(I,\,\pm L_{\rm lc}(I), \,\tau ) \,=\, 0. \end{equation}

(91)

The domain of the DF, f, when stars are lost to the black hole through an empty loss cone. The domain is the union of the two disjoint wedge-like regions that are symmetric about the I-axis. The tip of the upper wedge is at I = L = Ilc and the tip of the lower wedge is at I = −L = Ilc. The DF vanishes on the dashed lines corresponding to prograde (+) and retrograde (−) loss-cone boundaries (‘lcb’).

Figure 2.

The domain of the DF, f, when stars are lost to the black hole through an empty loss cone. The domain is the union of the two disjoint wedge-like regions that are symmetric about the I-axis. The tip of the upper wedge is at I = L = I_lc and the tip of the lower wedge is at I = −L = I_lc. The DF vanishes on the dashed lines corresponding to prograde (+) and retrograde (−) loss-cone boundaries (‘lcb’).

Open in new tab Download slide

From equation (62a) we know that the probability current density J(I, L, τ) = J^adv(I, L, τ) + J^dif(I, L, τ), where J^adv = Vf and J^dif = −D(∂f/∂L) as given in equations (62b) and (62c). Since f = 0 on the ±lcb, J^adv = 0 on these boundaries. Hence the total current at the ±lcb is purely diffusive:

\begin{eqnarray} J_+(I, \tau ) \,\equiv \, J(I,\,+L_{\rm lc}(I), \,\tau ) \,=\, -D_+\left(\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right)_+, \end{eqnarray}

(92a)

\begin{eqnarray} J_-(I, \tau ) \,\equiv \, J(I,\,-L_{\rm lc}(I), \,\tau ) \,=\, -D_-\left(\frac{\mathrm{\partial} f}{\mathrm{\partial} L}\right)_-. \end{eqnarray}

(92b)

Here (∂f/∂L)_± are evaluated on the ±lcb, and the loss-cone diffusion coefficients, D_±, are given by equation (62c):

\begin{eqnarray} D_+(I, \tau ) \,=\, \frac{1}{N} \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, f(I^{\prime }, L^{\prime }, \tau )\,\quad {for L \,=\, +L_{\rm lc}(I);} \end{eqnarray}

(93a)

\begin{eqnarray} D_-(I, \tau ) \,=\, \frac{1}{N} \int {\rm d}I^{\prime }\,{\rm d}L^{\prime }\,\delta (\Omega - \Omega ^{\prime })\,K\, f(I^{\prime }, L^{\prime }, \tau )\,\quad \mbox{for $L \,=\, -L_{\rm lc}(I)$,} \end{eqnarray}

(93b)

where the integral, ∫ dI΄ dL΄ , is over the entire domain of f, consisting of the two disjoint wedges in Fig. 2. In the integrand both Ω(I, L) and K(I, L, I΄, L΄) should be evaluated on the ±lcb (L = ±L_lc), so the loss-cone diffusion coefficients have contributions only from stars that are in apsidal resonance, Ω(I΄, L΄) = Ω(I, ±L_lc(I)), with a star on one of the lcb. It is evident that both D₊ and D₋ are non-negative quantities, as is f everywhere in its domain of definition. Since f = 0 on the ± lcb, we must have (∂f/∂L)₊ ≥ 0 and (∂f/∂L)₋ ≤ 0. Therefore, from equations (92a) and (92b), we have

\begin{equation} J_+(I, \tau ) \,\le \, 0\,\quad \mbox{and}\quad J_-(I, \tau ) \,\ge \, 0\,\quad \mbox{for $I \,\ge \, I_{\rm lc}$,} \end{equation}

(94)

which is as we might have expected. Fig. 3 shows the current density on all the boundaries of the domain of f.

The probability current density J(I, L, τ) at the boundaries. J = 0 on the upper and lower boundaries marking the prograde and retrograde circular orbits. J = J+(I, τ) ≤ 0 on the +lcb, and J = J−(I, τ) ≥ 0 on the −lcb.

Figure 3.

The probability current density J(I, L, τ) at the boundaries. J = 0 on the upper and lower boundaries marking the prograde and retrograde circular orbits. J = J₊(I, τ) ≤ 0 on the +lcb, and J = J₋(I, τ) ≥ 0 on the −lcb.

Open in new tab Download slide

This loss of stars introduces time-dependence in physical quantities that were conserved for the case of a lossless disc (studied in the previous section). The probability distribution P(I, τ) varies as

\begin{eqnarray} \frac{\mathrm{\partial} P}{\mathrm{\partial} \tau } \,=\, \int {\rm d}L\,\frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } \,=\, -\int {\rm d}L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,=\, -\int _{-I}^{-L_{\rm lc}}{\rm d}L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,-\, \int _{L_{\rm lc}}^{I}{\rm d}L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} = J_+(I, \tau ) \,-\, J_-(I, \tau ), \end{eqnarray}

(95)

because J(I, ±I, τ) = 0 from equation (59). From equation (94) we know that J₊ is non-positive and J₋ is non-negative. Hence (∂P/∂τ) ≤ 0 , as expected. Then the rate of mass fuelling of the MBH, |$\,\dot{{\cal M}}_{\rm lc}\,$|⁠, is the negative of the rate at which the disc mass changes:

\begin{equation} \dot{{\cal M}}_{\rm lc} \,=\, -\frac{{\rm d}{\cal M}}{{\rm d}\tau } \,=\, -M\int _{I_{\rm lc}}^\infty {\rm d}I\,\frac{\mathrm{\partial} P}{\mathrm{\partial} \tau } \,=\, M\int _{I_{\rm lc}}^\infty {\rm d}I\left\lbrace J_-(I, \tau ) \,-\, J_+(I, \tau )\right\rbrace \,\ge \, 0. \end{equation}

(96)

Similarly the rate at which angular momentum is fed to the MBH, |$\,\dot{{\cal L}}_{\rm lc}\,$|⁠, is equal to the negative of the rate of change of the disc angular momentum:

\begin{eqnarray} \dot{{\cal L}}_{\rm lc} = -\frac{{\rm d}{\cal L}}{{\rm d}\tau } &=& -M\int {\rm d}I\,{\rm d}L\, L\,\frac{\mathrm{\partial} f}{\mathrm{\partial} \tau } \,=\, M\int _{I_{\rm lc}}^\infty {\rm d}I\left\lbrace \int _{-I}^{-L_{\rm lc}}{\rm d}L\,L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,+\, \int _{L_{\rm lc}}^{I}{\rm d}L\,L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L}\right\rbrace \nonumber \\ &=& M\int _{I_{\rm lc}}^\infty {\rm d}I\left\lbrace \, {[}\,LJ\, {]}_{-I}^{-L_{\rm lc}} \,\,+\,\, {[}\,LJ\,{]}_{L_{\rm lc}}^{I}\,\right\rbrace \,\,-\,\, M\int {\rm d}I\,{\rm d}L\, J. \end{eqnarray}

(97)

From the first identity of equation (61) the last term on the right side is zero. From equation (59) we know that J vanishes when L = ±I. Hence the rate at which angular momentum is fed to the MBH is

\begin{equation} \dot{{\cal L}}_{\rm lc} \,=\, -M\int _{I_{\rm lc}}^\infty {\rm d}I\,L_{\rm lc}(I) \left\lbrace J_-(I, \tau ) \,+\, J_+(I, \tau )\right\rbrace. \end{equation}

(98)

Each star falling into the MBH through the empty loss cone brings with it the specific angular momentum corresponding to the value at its lcb. In Fig. 3, the upper region (shaded blue) feeds the MBH a specific angular momentum of L_lc per star, and the lower region (shaded red) feeds the MBH a specific angular momentum of −L_lc per star. |$\dot{{\cal L}}_{\rm lc}$| can be positive or negative, depending on the net contribution from the prograde and retrograde stellar populations.

The rate at which energy is fed to the MBH, |$\,\dot{{\cal E}}_{\rm lc}\,$|⁠, is equal to the negative of the rate at which the disc energy changes. Using equation (70) we have

\begin{eqnarray} \dot{{\cal E}}_{\rm lc} &=& -\frac{{\rm d}{\cal E}}{{\rm d}\tau } = -M\int {\rm d}I\,{\rm d}L\, \left\lbrace E_{\rm k}(I) \,+\,\varepsilon H^{\rm rel}(I, L) \,+\, \varepsilon \Phi ^{\rm ext}(I, L)\right\rbrace \frac{\mathrm{\partial} f}{\mathrm{\partial} \tau }- \frac{\varepsilon M}{2}\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\, C_0(I, L, I^{\prime }, L^{\prime })\left\lbrace f\,\frac{\mathrm{\partial} f^{\prime }}{\mathrm{\partial} \tau } + f^{\prime }\,\frac{\mathrm{\partial} f}{\mathrm{\partial} \tau }\right\rbrace \nonumber \\ [3ex] &=& M\int {\rm d}I\,{\rm d}L\, \left\lbrace E_{\rm k}(I) \,+\,\varepsilon H^{\rm rel}(I, L) \,+\, \varepsilon \Phi ^{\rm ext}(I, L)\right\rbrace \frac{\mathrm{\partial} J}{\mathrm{\partial} L} + \frac{\varepsilon M}{2}\int {\rm d}I\,{\rm d}L\,{\rm d}I^{\prime }\,{\rm d}L^{\prime }\, C_0(I, L, I^{\prime }, L^{\prime })\left\lbrace f\,\frac{\mathrm{\partial} J^{\prime }}{\mathrm{\partial} L^{\prime }} + f^{\prime }\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L}\right\rbrace. \end{eqnarray}

(99)

Since C₀(I, L, I΄, L΄) is symmetric under the interchange (I, L) ↔ (I΄, L΄) of integration variables, we can replace f(∂J΄/∂L΄) by f΄(∂J/∂L). Then,

\begin{eqnarray} \dot{{\cal E}}_{\rm lc} &=& M\int {\rm d}I\,{\rm d}L\, \left\lbrace E_{\rm k}(I) \,+\,\varepsilon H^{\rm rel}(I, L) \,+\, \varepsilon \Phi ^{\rm ext}(I, L) \,+\, \varepsilon \Phi (I, L)\right\rbrace \frac{\mathrm{\partial} J}{\mathrm{\partial} L} = M\int {\rm d}I\,{\rm d}L\, \left\lbrace E_{\rm k}(I) \,+\,\varepsilon H(I, L)\right\rbrace \frac{\mathrm{\partial} J}{\mathrm{\partial} L} \nonumber \\ [1ex] &=& M\int {\rm d}I\,E_{\rm k}(I)\int {\rm d}L\,\frac{\mathrm{\partial} J}{\mathrm{\partial} L} \,+\, \varepsilon M\int {\rm d}I\,{\rm d}L\, H(I, L)\frac{\mathrm{\partial} J}{\mathrm{\partial} L}= M\int _{I_{\rm lc}}^\infty {\rm d}I\,E_{\rm k}(I)\left\lbrace J_-(I, \tau ) \,-\, J_+(I, \tau )\right\rbrace \nonumber \\ &&+\, \varepsilon M\int _{I_{\rm lc}}^\infty {\rm d}I\left\lbrace \, {[}\,HJ\, {]}_{-I}^{-L_{\rm lc}} \,\,+\,\, {[}\,HJ\,{]}_{L_{\rm lc}}^{I}\,\right\rbrace \,\,-\,\, \varepsilon M\int {\rm d}I\,{\rm d}L\, \Omega J. \end{eqnarray}

(100)

By the second identity of equation (61) the last integral vanishes. Then, we have

\begin{equation} \dot{{\cal E}}_{\rm lc} \,=\, M\int _{I_{\rm lc}}^\infty {\rm d}I\,E_{\rm k}(I)\left\lbrace J_- \,-\, J_+\right\rbrace \,\,+\,\, \varepsilon M\int _{I_{\rm lc}}^\infty {\rm d}I\left\lbrace H_-\,J_- \,-\, H_+\,J_+\right\rbrace. \end{equation}

(101)

The first term is the feeding rate of the Kepler orbital energy to the MBH. It is the dominant term and has a negative value because E_k(I) < 0, J₋ ≥ 0 and J₊ ≤ 0 . In the second term, H_± = H(I, ±L_lc, τ) are values of the ring Hamiltonian at the ±lcb. The second term is smaller by O(ε) and gives the feeding rate of the ring orbital energy.

Equations (96), (98) and (101) give the rates at which RR feeds mass, angular momentum and energy to the MBH in simple and readily interpretable forms in terms of the currents of stars that pass through the loss-cone boundaries in (I, L) space. The boundary currents are diffusive in nature and given by equations (92a) and (92b), with the diffusion coefficients determined self-consistently by the DF as given in equations (93a) and (93b).

8 CONCLUSIONS

Our goal in Papers I through III is to formulate a theoretical framework describing the secular dynamics and statistical mechanics of a low-mass (or ‘Keplerian’) stellar system orbiting a MBH. In Papers I and II we derived the equations governing the collisionless and collisional (‘RR’) evolution, respectively, of general Keplerian systems. Here, in Paper III, we have simplified the RR kinetic equation and applied it to understand the physical kinetics of razor-thin axisymmetric Keplerian stellar discs.

RR evolution is driven by |$F^{(2)}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )$|⁠, the ‘irreducible part of the two-ring correlator’, which measures the ring–ring correlations that develop over time, due to mutual stochastic torquing. This is conveniently calculated in terms of a related quantity, |$W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau )$|⁠, the ‘wake of |${\cal R}^{\prime }$| at |${\cal R}$|’, which has a nice physical interpretation as the linear collisionless response of the stellar system to the act of removing any one ring and re-introducing it into any specified ring orbit which arrives at |${\cal R}^{\prime }$| at time τ . Even though the equation governing W is simpler than the one governing |$F^{(2)}_{\rm irr}$|⁠, it is still a linear, partial integro-differential equation in six variables. Our first result in this paper is the introduction of the ‘PRA’, which is an iterative scheme by which successive approximations to the wake can be generated. The lowest-order PRA ignores the effects of the gravitational polarization of the wake on its own evolution. The wake is then expressed as a time integral over the mutual interactions between two rings, over their entire past orbital histories (equation 26a). This can be used to derive the PRA kinetic equation (i.e. 29a–29c), which applies to a general Keplerian system; this equation can be thought of as the secular analogue of equations (1)–(3) of Polyachenko & Shukhman (1982).

We then applied the PRA to razor-thin axisymmetric discs described by DFs of the form F(I, L, τ). Equations (44a) and (44b) show that the PRA wake takes a simple form: with |${\cal R}= (I, L, g)$| and |${\cal R}^{\prime } = (I^{\prime }, L^{\prime }, g^{\prime })$|⁠, we have |$W({\cal R}\,\vert \,{\cal R}^{\prime }, \tau ) = - (\mathrm{\partial} F/\mathrm{\partial} L)\,\Delta L$|⁠, where |$\Delta L({\cal R}, {\cal R}^{\prime })$| is the angular momentum change experienced by a ring at |${\cal R}$| due to the torque exerted on it by the ring at |${\cal R}^{\prime }$| during the entire past history of their interaction. This form has a simple physical interpretation: W is the perturbation to F, accounting for the gain and loss of rings in an elemental volume around |${\cal R}$|⁠, due to torquing by |${\cal R}^{\prime }$|⁠. The wake can be decomposed into a ‘resonant’ and a ‘non-resonant’ part (equations 50a–50c), where the resonance referred to is between the apsidal precession rates of pairs of stellar orbits. Only the resonant part of the wake contributes to the collision integral, and determines RR evolution. This is physically reasonable, because we expect the accumulated mutual torque between two rings to average to zero when they are not resonant.

We then derive the PRA kinetic equation governing the RR evolution of F. This is given in equations (57) and (58) in the ‘conservation’ form’, where the probability density current of the RR flow in (I, L) space is seen to get contributions only from the resonant orbits. Since the I of every ring is constant, the current density is directed only along the L direction. The PRA kinetic equation can also be put in a self-consistent Fokker–Planck form, equation (63), with the diffusion coefficients determined self-consistently by the DF, as given by equations (64), (62b) and (62c).¹¹ The nice thing about the Fokker–Planck form is the ready interpretation of the diffusion coefficients in terms of the average characteristics, per unit time, of the random changes in the specific angular momentum of a ring (see equation 65).

The broad features of the physical kinetics of RR were studied for the two major cases of interest, corresponding to the two different boundary conditions under which the kinetic equation is to be solved. These are (a) ‘Lossless’ discs in which the MBH is not a sink of stars (Section 6), and (b) discs which feed stars to the MBH through an ‘empty loss cone’ (Section 7). We obtained the following results.

Since the MBH does not consume stars, total disc mass, angular momentum and energy are all conserved. A general H-function can either increase or decrease during RR, but the Boltzmann entropy is (essentially) unique in being a non-decreasing function of time. It follows that maximum entropy states are also states of secular thermal equilibrium, and their DFs have the Boltzmann form, investigated by Touma & Tremaine (2014) for the limit in which all rings share the same semi-major axis. As a bonus we get a formula for the two-ring correlation function at equilibrium.
A star that comes close enough to the MBH will be either torn apart by its tidal field or swallowed whole by it. The stochastic evolution of the orbital angular momenta during RR evolution can make a star lose enough angular momentum for its pericentre to come closer than some minimum distance from the MBH; then the star is counted as lost to the MBH within a Kepler orbital period (which is essentially instantaneous when compared with secular or RR time-scales). This is modelled by an ‘Empty loss-cone’ boundary condition in (I, L) space, through which stars are lost to the MBH (see Figs 2 and 3). Equations (96), (98) and (101) give the rates at which RR feeds mass, angular momentum and energy to the MBH. These are given in terms of the currents of stars that pass through the loss-cone boundaries. The boundary currents are diffusive in nature and given by equations (92a) and (92b), with the diffusion coefficients determined self-consistently by the DF and given in equations (93a) and (93b).

All of these properties are valid for axisymmetric discs, allowing for 1 PN and 1.5 PN general relativistic effects on apsidal precession, as well as the gravity of an arbitrary axisymmetric external potential.

Having established an overall picture of the physical kinetics, it is necessary to compute and understand how RR affects orbits of different semi-major axes and eccentricities. Below are three interesting ways in which the present work can be taken forward:

Touma & Tremaine (2014) solved for maximum entropy equilibria of self-gravitating Keplerian discs, in the case when all stars have the same semi-major axis. Both axisymmetric and non-axisymmetric secular thermal equilibria were included. They found that a significant fraction of the axisymmetric equilibria were not really states of maximum entropy but saddle points, and are thermodynamically¹² unstable to bifurcations favouring lopsided discs at the same energy and angular momentum. The symmetry breaking equilibria were solved for numerically and found to be generically rotating. The axisymmetric equilibria which were thermodynamically unstable were also found to be dynamically unstable, in the sense of the linearized ring CBE (whereas thermodynamic stability implies dynamical stability, thermodynamic instability does not imply dynamical instability). When followed into the non-linear regimes through numerical simulations, the dynamical instabilities saturated and relaxed into lopsided, uniformly precessing configurations. To be able to describe complex processes such as these requires (i) using the methods of Paper I to first construct generic lopsided and rotating collisionless equilibria which are also dynamically stable, and (ii) generalizing the approach of the present work to derive a kinetic equation for the RR of these broken-symmetry equilibria.
The RR kinetic theory of Paper II applies to general stellar distributions, and needs to be developed for fully three-dimensional systems. The simplest are spherical systems around a non-spinning MBH, whose orbital structure is as simple as the axisymmetric discs we studied in this paper. The difference now is that we need to consider the gravitational interactions between pairs of mutually inclined stellar orbits. This adds to the technical complexity of Fourier expansions, but the physical kinetics of RR can be expected to be as simple as it is for two-dimensional axisymmetric discs: we expect to find a kinetic equation of the Fokker–Planck form, with an H-theorem and secular thermal equilibrium of the Boltzmann form for a lossless system, and diffusive loss-cone feeding of mass and energy to the MBH.
Three-dimensional axisymmetric systems are naturally rotating, and we expect RR to feed the MBH with mass, energy and angular momentum. It is also fortuitous that the secular dynamics of these systems is completely integrable with I, L_z and H(I, L, L_z, g) being three isolating integrals (Sridhar & Touma 1999; Sambhus & Sridhar 2000).

Acknowledgments

We are grateful to Jerome Perez, Stephane Colombi and the Institut Henri Poincaré for hosting us when a part of this work was done. We thank Scott Tremaine for comments on an earlier draft.

1

The typical angular momentum (per unit mass) of a test ring is |$L\sim \sqrt{GM_\bullet a}$| where a is its semi-major axis. The torque experienced by the test ring, due to all the other rings, is |${\cal N}\sim (Gm_\star /a)N^{1/2}$|⁠. Acting over a coherence time T_sec, this causes a change in L of the order of |${\cal N} T_{\rm sec}$|⁠. Over longer times t ≫ T_sec the net change in L is |${\sim\, }{\cal N} T_{\rm sec}(t/ T_{\rm sec})^{1/2}$|⁠, which is of the order of L when t ∼ NT_sec = T_res.

2

The BBGKY hierarchy applies to many-body systems and has been used to derive the well-known Boltzmann equation for dilute gases and the Balescu–Lenard equation for electrostatic plasmas (Liboff 2003). One begins with the N-particle DF obeying the Liouville equation in 6N-dimensional phase space and integrates this over various phase space coordinates. This gives a hierarchy of equations for ‘reduced’ DFs, where the equation for the s-particle DF (with 1 ≤ s ≤ N − 1) depends on the (s + 1)-particle DF. The aim is to derive a closed equation for the one-particle DF in six-dimensional phase space by closing the hierarchy of equations by exploiting a small parameter.

3

Here, |$\Psi ({\cal R}, {\cal R}^{\prime })$| is proportional to the gravitational interaction potential energy between two rings, |${\cal R}$| and |${\cal R}^{\prime }$|⁠, and has been scaled in a manner such that the equations of motion contain only quantities of order unity.

4

The typical time-scale of growth of a secular instability is T_sec, which is much shorter than the RR time-scale T_res = NT_sec. Hence, the collisional terms, as discussed in Section 2.2, may be neglected while considering the dynamical stability of a secular equilibrium.

5

In the course of RR evolution, it may happen that a system goes through a sequence of stable secular equilibria and arrives at one which is secularly unstable. Then RR evolution would have tipped the system into an instability which will grow and saturate (due to collisionless relaxation) on the short time-scale T_sec; clearly the collision terms can be ignored while describing this fast process. The new equilibrium state is evidently stable, and subsequent RR evolution is described as given above.

6

The decomposition of the collision integral into dissipative and fluctuating components is traditional, and follows from the theory for a general stellar system. In Gilbert's words (see section V of Gilbert 1968), ‘The first of these clearly represents the deceleration of a star by its associated polarization cloud. The second term represents the effect on a star of the rapidly fluctuating, random fields produced by its neighbors.’ We are dealing with the secular version which is valid for Keplerian stellar systems.

7

This is the convention for the symbols (L, g) we use in the rest of this paper.

8

Hence the theory does not apply to very cold discs violating the Toomre (1964) criterion, because they are unstable to axisymmetric modes growing over times ∼ T_kep. So the discs we study are understood to have velocity dispersions large enough to make them stable to fast modes.

9

As earlier the normalization of equation (54) is true for all τ only in the conservative case of no loss of stars to the MBH. When the MBH is a sink of stars, the above normalization holds only at some initial reference time τ = τ₀ when the disc mass is M.

10

These are related to the original Lagrange multipliers by ln A = ( − 1 + B − β΄ME_k), β = β΄Mε and γ = γ΄M .

11

This is a special case of a more general result derived in Polyachenko & Shukhman (1982): the lowest-order PRA kinetic equation for a general stellar system with integrable mean-field Hamiltonian can be expressed in a self-consistent Fokker–Planck form.

12

Here, ‘thermodynamically’ refers to the thermodynamics of a large number of Gaussian rings. Since the semi-major axis of every ring is a fixed constant, the phase space is finite and the problem is well-defined, unlike the case of a general self-gravitating system.

REFERENCES

Balescu

R.

1960

,

Phys. Fluids

,

3

,

52

Binney

J.

,

Tremaine

S.

2008

,

Galactic Dynamics

2nd edn

Princeton Univ. Press

Princeton, NJ

OpenURL Placeholder Text

Gilbert

I. H.

1968

,

ApJ

,

152

,

1043

Landau

L. D.

1936

,

Phys. Z. Sowj. Union

,

10

,

154

Lenard

A.

1960

,

Ann. Phys.

,

10

,

390

Liboff

R. L.

2003

,

Kinetic Theory: Classical, Quantum, and Relativistic Descriptions

Springer-Verlag

New York

OpenURL Placeholder Text

Lifshitz

E. M.

,

Pitaevskii

L. P.

1981

,

Physical Kinetics

1st edn

Pergamon Press

Oxford

OpenURL Placeholder Text

Murray

C. D.

,

Dermott

S. F.

1999

,

Solar system Dynamics

Cambridge University Press

Cambridge

OpenURL Placeholder Text

Polyachenko

V. L.

,

Shukhman

I. G.

1982

,

Sov. Astron.

,

26

,

140

Polyachenko

E. V.

,

Polyachenko

V. L.

,

Shukhman

I. G.

2007

,

MNRAS

,

379

,

573

Rauch

K. P.

,

Tremaine

S.

1996

,

New Astron.

,

1

,

149

Rostoker

N.

1964

,

Phys. Fluids

,

7

,

491

Sambhus

N.

,

Sridhar

S.

2000

,

ApJ

,

542

,

143

Sridhar

S.

,

Touma

J.

1999

,

MNRAS

,

303

,

483

Sridhar

S.

,

Touma

J. R.

2016a

,

MNRAS

,

458

,

4129

(Paper I)

Sridhar

S.

,

Touma

J. R.

2016b

,

MNRAS

,

458

,

4143

(Paper II)

Toomre

A.

1964

,

ApJ

,

139

,

1217

Touma

J.

,

Tremaine

S.

2014

,

J. Phys. A

,

47

,

292001

Tremaine

S.

2005

,

ApJ

,

625

,

143

APPENDIX A

Here, we prove that the non-resonant part of the wake does not contribute to the collision integral of equation (52). We first rewrite equation (51c) as

\begin{equation} F^{(2,{\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau ) \,=\, \Gamma (I, L, I^{\prime }, L^{\prime })\, \left\lbrace \Psi ({\cal R}, {\cal R}^{\prime }) - C_0\right\rbrace , \end{equation}

(A1)

where the Fourier coefficient C₀ = C₀(I, L, I΄, L΄), the ring–ring interaction potential Ψ is a function of |${\cal R}= (I, L, g)$| and |${\cal R}^{\prime } = (I^{\prime }, L^{\prime }, g^{\prime })$|⁠, and the function

\begin{equation} \Gamma (I, L, I^{\prime }, L^{\prime }) \,=\, -\frac{1}{\Omega - \Omega ^{\prime }}\left\lbrace F\frac{\mathrm{\partial} F^{\prime }}{\mathrm{\partial} L^{\prime }} \,-\, F^{\prime }\frac{\mathrm{\partial} F}{\mathrm{\partial} L}\right\rbrace. \end{equation}

(A2)

Then we manipulate the Poisson Bracket (over the |${\cal R}$| variables), |$\left[\Psi , F^{(2,{\rm n})}_{\rm irr}\right]$| as follows:

\begin{eqnarray} \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,{\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right] &=& \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,\Gamma \left\lbrace \Psi ({\cal R}, {\cal R}^{\prime }) - C_0\right\rbrace \,\right] =\left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,\Gamma \Psi ({\cal R}, {\cal R}^{\prime })\,\right] \,-\, \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,\Gamma C_0\,\right]\nonumber \\ &=& \Psi \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,\Gamma \,\right] \,-\, \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\,\Gamma C_0\,\right]= \Psi \,\frac{\mathrm{\partial} \Psi }{\mathrm{\partial} g}\,\frac{\mathrm{\partial} \Gamma }{\mathrm{\partial} L} \,-\, \frac{\mathrm{\partial} \Psi }{\mathrm{\partial} g}\,\frac{\mathrm{\partial} (\Gamma C_0)}{\mathrm{\partial} L}\nonumber\\ & =& \frac{\mathrm{\partial} }{\mathrm{\partial} g} \left\lbrace \frac{\Psi ^2}{2}\frac{\mathrm{\partial} \Gamma }{\mathrm{\partial} L} \,-\, \Psi \frac{\mathrm{\partial} (\Gamma C_0)}{\mathrm{\partial} L}\right\rbrace. \end{eqnarray}

(A3)

The expression in the curly bracket depends on g and g΄ only through the function |$\Psi ({\cal R}, {\cal R}^{\prime })$|⁠, in which the apsidal angles occur only in the combination (g − g΄). Hence,

\begin{equation} \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,{\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right] \,=\, -\frac{\mathrm{\partial} }{\mathrm{\partial} g^{\prime }} \left\lbrace \frac{\Psi ^2}{2}\frac{\mathrm{\partial} \Gamma }{\mathrm{\partial} L} \,-\, \Psi \frac{\mathrm{\partial} (\Gamma C_0)}{\mathrm{\partial} L}\right\rbrace , \end{equation}

(A4)

which implies that

\begin{equation} \int _0^{2\pi}{\rm d}g^{\prime } \left[\,\Psi ({\cal R}, {\cal R}^{\prime }),\, F^{(2,{\rm n})}_{\rm irr}({\cal R}, {\cal R}^{\prime }, \tau )\,\right] \,=\, 0. \end{equation}

(A5)

Therefore the non-resonant part of the wake does not contribute to the collision integral of equation (52), and plays no role in the kinetic evolution of the system.

© 2016 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society

Download all slides

Views

405

Altmetric

Total Views 405

262 Pageviews

143 PDF Downloads

Since 12/1/2016

Month:	Total Views:
December 2016	1
January 2017	6
February 2017	5
March 2017	14
May 2017	1
June 2017	7
July 2017	6
August 2017	4
September 2017	8
October 2017	5
November 2017	1
December 2017	1
January 2018	2
February 2018	5
March 2018	2
April 2018	4
May 2018	2
June 2018	5
July 2018	3
August 2018	3
December 2018	3
February 2019	2
March 2019	5
April 2019	1
June 2019	2
August 2019	1
October 2019	5
November 2019	1
December 2019	2
March 2020	1
April 2020	5
May 2020	3
June 2020	6
July 2020	9
August 2020	8
September 2020	3
October 2020	1
November 2020	1
December 2020	3
January 2021	1
February 2021	4
March 2021	4
April 2021	4
May 2021	2
June 2021	1
July 2021	1
August 2021	4
September 2021	2
October 2021	5
November 2021	7
December 2021	3
January 2022	4
February 2022	5
March 2022	1
April 2022	3
July 2022	4
August 2022	15
September 2022	4
October 2022	1
December 2022	11
January 2023	2
February 2023	1
March 2023	5
April 2023	4
May 2023	2
June 2023	1
July 2023	1
August 2023	18
September 2023	15
October 2023	14
November 2023	12
December 2023	5
January 2024	9
February 2024	14
March 2024	11
April 2024	5
May 2024	3
June 2024	5
July 2024	5
August 2024	4
September 2024	3
October 2024	3
November 2024	3
December 2024	2
January 2025	3
February 2025	7
March 2025	11
April 2025	1
May 2025	3