Stellar Dynamics around a Massive Black Hole I:
Secular Collisionless Theory
Abstract
We present a theory in three parts, of the secular dynamics of a (Keplerian) stellar system of mass M orbiting a black hole of mass M_{\bullet}\gg M. Here we describe the collisionless dynamics; Papers II and III are on the (collisional) theory of Resonant Relaxation. The mass ratio, \varepsilon=M/M_{\bullet}\ll 1, is a natural small parameter implying a separation of time scales between the short Kepler orbital periods and the longer orbital precessional periods. The collisionless Boltzmann equation (CBE) for the stellar distribution function (DF) is averaged over the fast Kepler orbital phase using the method of multiple scales. The orbit–averaged system is described by a secular DF, F, in a reduced phase space. F obeys a secular CBE that includes stellar self–gravity, general relativistic corrections up to 1.5 post–Newtonian order, and external sources varying over secular times. Secular dynamics, even with general time dependence, conserves the semi–major axis of every star. This additional integral of motion promotes extra regularity of the stellar orbits, and enables the construction of equilibria, F_{0}, through a secular Jeans theorem. A linearized secular CBE determines the response and stability of F_{0}. Spherical, non–rotating equilibria may support long–lived, warp–like distortions. We also prove that an axisymmetric, zero–thickness, flat disc is secularly stable to all in–plane perturbations, when its DF, F_{0}, is a monotonic function of the angular momentum at fixed energy.
keywords:
galaxies: kinematics and dynamics — galaxies: nuclei — Galaxy: center1 Introduction
Massive black holes (MBHs), consuming gas and stars, are believed to be responsible for the radiant power of quasars and other active galactic nuclei (Rees, 1984; Krolik, 1999). There is also dynamical evidence for MBHs with masses M_{\bullet} in the range 10^{5}10^{10}\,M_{\odot}\, in the nuclei of more nearby galaxies (Kormendy & Ho, 2013). Coevolution (or otherwise) of a MBH and its host galaxy has been discussed using many studies of correlations of M_{\bullet} with galaxy properties (Kormendy & Ho, 2013; Heckman & Best, 2014). Ground–based observations provided the first dynamical evidence for MBHs in nearby galaxies (Kormendy & Richstone, 1995). They also revealed nuclear star clusters where stellar number densities far exceeded those of globular clusters. Early work on the stellar dynamics of galactic nuclei considered several aspects of the interactions between the MBH and its stellar environment: feeding of stars to a MBH through two–body relaxation (Bahcall & Wolf, 1976; Cohn & Kulsrud, 1978); the effect of the adiabatic growth of a MBH on the formation of stellar density cusps (Young, 1980), and on stellar orbits (Goodman & Binney, 1984); destruction by a MBH of the regularity of stellar orbits in the inner parts of triaxial galaxy (Gerhard & Binney, 1985). With the Hubble Space Telescope spatially resolved kinematics, of the dense nuclear clusters of stars populating the MBH’s sphere of influence, became available. Recent ground–based observations with adaptive optics have resolved the sphere of influence of the MBH at our Galactic Centre so well, that the orbits of many individual stars are now being followed (Genzel, Eisenhauer & Gillessen, 2010).
Our interest is in these dense stellar systems in galactic nuclei lying within the sphere of influence of the MBH. Following Tremaine (2005) we may consider a stellar system orbiting a MBH to be Keplerian^{1}^{1}1Tremaine (2005) uses the term “near–Keplerian”; we have shortened it to Keplerian. when the gravitational force, over a large range of radii, is dominated by the MBH, and the dominant non–Keplerian force arises from the self–gravity of all the stars. In the absence of self–gravity the orbit of a star follows a closed Kepler ellipse around the MBH, because the radial and azimuthal frequencies are equal to each other (i.e. they are “degenerate” or “resonant”). When the small self–gravity of the cluster is taken into account, the orbit of a star is usefully represented by its instantaneous osculating Kepler ellipse, which evolves slowly over secular time scales that are much longer than the fast orbital times around the MBH. Progress in the understanding of the dynamics of Keplerian stellar systems can perhaps be traced to two papers: Tremaine (1995) on the structure of the (then) double nucleus of the Andromeda galaxy (M31), and Rauch & Tremaine (1996) on the role of resonant relaxation in stellar systems with degenerate (or resonant) frequencies. The former has stimulated much work on (collisionless) dynamics whereas the latter proposed a new and efficient (collisional) relaxation process called Resonant Relaxation, which is the statistical mechanics of Keplerian stellar systems. Two decades of research has paved the way for the construction of triaxial, even lopsided, equilibria around MBHs, and simplified the modeling, stability analysis, and numerical simulation of nuclear star clusters. We now have a better appreciation of the dynamics and the statistical mechanics of Keplerian star clusters, and of the additional role of the MBH as a consumer of stars. But the whole subject lacks a clearly formulated and general mathematical foundation. Our aim is to fill this gap and provide a secular theory of the structure and evolution of Keplerian stellar systems, described by the slow time evolution of a secular distribution function (DF), defined in a reduced phase space. Here (Paper I) we derive the (collisionless) theory describing their dynamics, beginning with the standard description of collisionless stellar systems given in Binney & Tremaine (2008). In Sridhar & Touma (2015; Paper II) we deal with their (collisional) statistical mechanics and present a first–principles theory of resonant relaxation based on Gilbert (1968). This is then applied in Sridhar & Touma (2016; Paper III) to the physical kinetics of axisymmetric, zero–thickness, flat discs.
The sketch of a model offered in Tremaine (1995) — for a lopsided Keplerian disc of stars around a MBH in the nucleus of M31 — has since been fleshed out in a number of numerical models (Statler, 1999; Bacon et al., 2001; Salow & Statler, 2001; Sambhus & Sridhar, 2002; Peiris & Tremaine, 2003; Salow & Statler, 2004; Kazandjian & Touma, 2013; Brown & Magorrian, 2013). This has brought to the forefront deep and interesting questions regarding the general nature of secular collisionless evolution. Below we give a brief time–line of developments in the secular collisionless dynamics of Keplerian stellar systems. Work in this field may be broadly divided into explorations of orbital dynamics, construction of stationary DFs, linear stability of stationary DFs, and the nonlinear evolution of the DF.
Orbital Dynamics: The orbital structure of model galactic–type potentials with a central MBH was studied in Sridhar & Touma (1997a, b); Merritt & Valluri (1999). Sridhar & Touma (1999) imported the Gaussian orbit–averaging methods of celestial mechanics into stellar dynamics within the sphere of influence of a MBH. In this procedure a stellar orbit is averaged over the (fast–varying) orbital phase of the corresponding osculating Kepler ellipse, a direct consequence of which is the near conservation of the semi–major axis of the stellar orbit. Since there arises an extra (approximate) integral of motion, they concluded that a MBH promotes extra order on stellar orbits within their sphere of influence. Orbit–averaging takes us to a reduced phase space describing secular dynamics: this is the theatre of the coupled evolution of the orbital angular momentum and Lenz vectors — giving the eccentricity and spatial orientation of the orbit — keeping fixed the semi–major axis. Studying orbits in planar, lopsided model disc potentials, they found a family of stable orbits that were lopsided in the same sense as the disc potential. The more complex orbital structure of triaxial systems was also studied in a number of papers: see e.g. Sambhus & Sridhar (2000); Poon & Merritt (2001); Merritt & Vasiliev (2011). Of particular interest is a family of centrophilic orbits that bring stars close to the MBH. These are reviewed in Merritt (2013), which also offers a treatment of the post–Newtonian (relativistic) effects of the MBH on stellar orbital dynamics.
Linear Stability of Stationary DFs: Stationary DFs are often taken as first approximations for describing the nuclear star clusters, and their stability to small perturbations is an important question. This is, in general, a difficult problem because the long–range nature of gravitational interactions leads to an integral eigenvalue problem. Two types of stellar systems have been studied: (a) Zero–thickness flat discs, and (b) Spherical non–rotating systems.
(a) Linear Stability of Discs: Sridhar, Syer & Touma (1999) used the Laplace–Lagrange secular theory to study slow, lopsided m=1 modes (where m is the azimuthal wavenumber) of cold stellar discs around MBHs. Lee & Goodman (1999) considered nonlinear m=1 spiral density waves in fluid discs, and used a variational method to derive the dispersion relation and angular momentum fluxes. The numerical simulations of Jacobs & Sellwood (2001) provided evidence for long–lived lopsided modes. Tremaine (2001) demonstrated that a variety (fluid, stellar and softened–gravity) of rotating discs were neutrally stable to m=1 modes in the Wentzel–Kramers–Brillouin (WKB) limit. He formulated a linear integral eigenvalue problem for softened–gravity discs, and demonstrated that all m=1 modes were neutrally stable; the problem was also solved numerically to determine the eigenfrequencies and eigenfunctions. Jalali & Tremaine (2012) derived the WKB dispersion relationship for slow modes in a collisionless disc with a Schwarzschild DF and showed that modes of all m were stable. Touma (2002) presented a Laplace–Lagrange theory for discs composed of softened Gaussian wires, confirmed the stability result of Tremaine (2001) for prograde discs, and also demonstrated that counter–rotating discs were unstable to slow, lopsided modes; a possible origin of the lopsided nuclear disc of M31 may be through such an instability. Sambhus & Sridhar (2002) built on this to include a retrograde population of stars in their planar model of the nuclear disc of M31, and suggested that this population could arise naturally through an accretion event. Tremaine (2005) considered stellar discs with DFs even in the angular momentum, and empty loss cones (i.e. DF is zero at zero angular momentum). Using the variational principle of Goodman (1988) he argued that many of them are likely to be unstable to m=1 modes. Polyachenko, Polyachenko & Shukhman (2007) considered mono–energetic discs dominated by nearly radial orbits, and found that they are prone to loss cone instabilities of all m, if there is some amount of counter–rotating stars. Sridhar & Saini (2010) derived the instability of m=1 modes of softened–gravity, counter–rotating discs in the WKB limit. Gulati, Saini & Sridhar (2012) studied the linear integral eigenvalue problem for these discs and determined eigenvalues, precession frequencies and growth rates.
(b) Linear Stability of Spherical Non–rotating Systems: Tremaine (2005) formulated the eigenvalue problem for the linear secular modes of spherical non–rotating systems and determined the following: (i) Modes with arbitrary l (where l is the latitudinal wavenumber) are of two types — either pairs of modes that are damped and growing in time, or pairs of oscillatory modes; (ii) Lopsided l=1 modes are stable when the DF is either a decreasing function of angular momentum (at fixed energy), or the DF is an increasing function of angular momentum and has a non–empty loss cone; (iii) DFs that are increasing functions of angular momentum with an empty loss cone are neutrally stable to an l=1 mode that corresponds to a uniform displacement. Polyachenko, Polyachenko & Shukhman (2007) considered a loss–cone instability for mono–energetic spherical stellar systems dominated by nearly radial orbits. They found that DFs which are non monotonic functions of the angular momentum may be unstable to modes with l\geq 3 linear perturbations; a necessary condition for instability is the retrograde precession of orbits, which obtains naturally in such clusters.
Nonlinear evolution of the DF: Numerical studies of the collisionless dynamics of Keplerian stellar systems have been more limited than their galactic counterparts. This is partly because observations are not yet precise enough to finely resolve the sphere of influence of central black holes (the Galactic center and M31’s nucleus being the exceptions), and partly because the kind of N–body simulations required to resolve the collisionless regime are computationally costly when a massive central body is present. Schwarzschild–type methods have been used to construct non–axisymmetric equilibria, both kinematic (Peiris & Tremaine, 2003) and selfconsistent (Poon & Merritt, 2001; Sambhus & Sridhar, 2002; Brown & Magorrian, 2013)). Early numerical simulations (Bacon et al., 2001; Jacobs & Sellwood, 2001) showed collisionless relaxation of perturbed stellar discs. Touma, Tremaine & Kazandjian (2009) built a numerical algorithm which is particularly suited to the numerical exploration of secular clusters. This is based on a generalization of an orbit–averaging algorithm, first introduced by Gauss, to systems with Plummer–softened interactions. The simulations proceed over secular time steps, in a vectorial framework which avoids singularities in the usual Keplerian elements. Averaging of force evaluations, which is costly, is ideal for parallel computations on a Beowulf cluster with or without GPU option. Their work explored various violent regimes of secular interactions, and provides the first evidence for collisionless relaxation of secularly unstable counter–rotating clusters into lopsided configurations. The (fully nonlinear) collisionless Boltzmann equation (CBE) for secular evolution was presented in Touma & Sridhar (2012), and used to explore the nonlinear, collisionless dynamics of counter–rotating discs. Kazandjian & Touma (2013) performed N–body simulations of unstable counter–rotating 3–dim clusters, and demonstrated that they relax collisionlessly into triaxial lopsided configurations; these appear to be promising models of the lopsided stellar nucleus of M31.
We are now in a position to bring all the work together under a general framework, describing the secular collisionless evolution of Keplerian stellar systems. We begin in § 2 by casting the standard collisionless Boltzmann equation for the phase space DF in the Delaunay action–angle variables (which are natural phase space coordinates for the pure Kepler problem). In § 3 the method of multiple scales is used to derive an explicit functional form for the DF. Complementing this is Appendix A on canonical perturbation theory applied to an individual stellar orbit. In § 4 self–gravity is augmented by relativistic corrections and external gravitational perturbations, to obtain the secular CBE governing the evolution of the secular DF. A secular Jeans theorem is proved in § 5, and used to discuss secular equilibria. Linear secular perturbations and stability are discussed in § 6, where some results are presented for two types of systems: (a) Spherical non–rotating systems; (b) Zero–thickness flat axisymmetric discs. Concluding remarks are offered in § 7.
2 Collisionless evolution of Keplerian Stellar Systems
In a Keplerian stellar system the gravitational force over a large range of radii is dominated by the MBH, and the dominant non–Keplerian force arises from the self–gravity of all the stars. The MBH also plays an additional role in destroying stars on orbits with near–zero angular momentum, through tidal disruption or being swallowed whole. Let us consider such a system composed of a large number (N\gg 1) stars, of total mass M and size R, orbiting a MBH of mass M_{\bullet}. Then we must have (1) \varepsilon\equiv(M/M_{\bullet})\ll 1\,, and (2) \,R\gg r_{\bullet}\equiv 2GM_{\bullet}/c^{2} such that (r_{\bullet}/R)\ll\varepsilon; this ensures that, over most of the cluster, stellar self–gravity dominates relativistic corrections to the MBH’s gravity. Then \varepsilon is the natural small parameter for the dynamics of non–relativistic star clusters.
The important time scales in the life of a Keplerian stellar system are: (i) The short Kepler orbital timescale, \mbox{$T_{\rm kep}$}=2\pi\left(R^{3}/GM_{\bullet}\right)^{1/2}\,; (ii) The longer secular timescale, \mbox{$T_{\rm sec}$}\,, over which the stars behave as if they were Keplerian Ellipses being torqued collisionlessly by the mean gravitational field of the cluster; (iii) The resonant relaxation timescale, T_{\rm res}\,, over which the angular momenta (but not the energies) of stellar orbits diffuse due to gravitational collisions between stars; (iv) The 2–body relaxation timescale, T_{\rm relax}, over which both the energies and angular momenta of the stellar orbits diffuse, due to gravitational collisions between stars. The different time scales are related to each other by, \,\mbox{$T_{\rm kep}$}:\mbox{$T_{\rm sec}$}:T_{\rm res}:T_{\rm relax}\,=\,1:% \varepsilon^{1}:N\varepsilon^{1}:N\varepsilon^{2}\,.^{2}^{2}2This estimate of T_{\rm res} applies to the diffusion of the magnitude of angular momentum, or scalar resonant relaxation. Logarithmic corrections have been neglected in the estimate for T_{\rm relax}. For galactic nuclei of interest T_{\rm relax} is usually greater than the Hubble time, so the usual 2–body relaxation process is not important. However, T_{\rm res} can be less than the Hubble time, so there are three regimes of interest:

Short times \left(t\,\lesssim\,\mbox{few}\;\mbox{$T_{\rm kep}$}\,\ll\mbox{$T_{\rm sec}$}\,\right) during which there are variations in the stellar distribution over Kepler orbital times. These have indeed been directly observed in our Galactic Centre. The general theory is undeveloped and could be important for understanding the first phase of the collisionless relaxation of a star cluster, were it to form in the sphere of influence of the MBH. This regime is not studied in this paper.

Intermediate times \left(\mbox{$T_{\rm kep}$}\,\ll\,\,t\,\sim\,\mbox{few}\;\mbox{$T_{\rm sec}$}\,% \,\ll\,T_{\rm res}\right) over which the cluster behaves collisionlessly. This secular regime is the focus of the present work (Paper I).

Long times \left(\mbox{$T_{\rm sec}$}\,\ll\,\,t\,\sim\,\mbox{few}\;T_{\rm res}\,\ll\,T_{% \rm relax}\right) over which the system evolves collisionally due to resonant relaxation. This is explored in Papers II and III.
The two Keplerian systems that have been studied best are also the closest ones, in which the sphere of influence of the MBH is well–resolved (Kormendy & Ho, 2013). The Galactic Centre MBH has an estimated mass of \sim 4.6\times 10^{6}~{}M_{\odot} (Yelda et al., 2014), and is surrounded by many young and old stars, with \varepsilon\sim 10^{2}10^{1} depending on what substructure is being studied. The ‘triple nucleus’ of M31 has a MBH of mass \sim 1.4\times 10^{8}~{}M_{\odot} (Bender et al., 2005), with a similar values of \varepsilon, depending on where the disc edge is located. One may ask if these values of \varepsilon are small enough. For M31 the answer seems to be in the affirmative because, as discussed earlier, the Keplerian disc model provides a natural explanation of the lopsidedness in terms of the trapping of orbits by self–gravity. The orbits of many individual stars close to the Galactic Centre MBH have been determined, and for these the dynamics is certainly Keplerian. Disc–like structures have also been observed, e.g. Yelda et al. (2014), but secular theory has not progressed to the stage of making clear predictions that can be tested.
2.1 Collisionless Boltzmann equation in the Delaunay variables
Let r be the position of a test star relative to the MBH, t be the time variable, and \mbox{\boldmath$u$}=({\rm d}\mbox{\boldmath$r$}/{\rm d}t) be the relative velocity. A large N stellar system can be described at time t by \,f(\mbox{\boldmath$r$},\mbox{\boldmath$u$},t), the distribution function (DF) in the 6–dim single–particle phase space \{\mbox{\boldmath$r$},\mbox{\boldmath$u$}\}. The DF is positive and, when the MBH is not considered a sink of stars, the total mass of the stellar system is conserved. Then the DF can be normalized for all time:
\int f(\mbox{\boldmath$r$},\mbox{\boldmath$u$},t)\,{\rm d}\mbox{\boldmath$r$}% \,{\rm d}\mbox{\boldmath$u$}\;=\;1,  (1) 
making f a (single–particle) probability distribution function. The gravitational force experienced by a star is due to the MBH and all the other stars. Per unit mass the former is \,GM_{\bullet}\hat{\mbox{\boldmath$r$}}/r^{2}\,, and the latter can be written as \,\mbox{$\partial$}(\varepsilon\varphi)/\mbox{$\partial$}\mbox{\boldmath$r$}\, where \,\varepsilon\varphi(\mbox{\boldmath$r$},t)\, is the mean–field self–gravitational potential given by,
\varphi(\mbox{\boldmath$r$},t)\;=\;\,GM_{\bullet}\int\frac{f(\mbox{\boldmath$% r$}^{\prime},\mbox{\boldmath$u$}^{\prime},t)}{\mbox{\boldmath$r$}\mbox{% \boldmath$r$}^{\prime}}\,{\rm d}\mbox{\boldmath$r$}^{\prime}\,{\rm d}\mbox{% \boldmath$u$}^{\prime}\,.  (2) 
We write the acceleration of the MBH due to all the stars as \varepsilon\mbox{\boldmath$A$}_{\bullet}(t), where
\mbox{\boldmath$A$}_{\bullet}(t)\;=\;GM_{\bullet}\int f(\mbox{\boldmath$r$},% \mbox{\boldmath$u$},t)\,\frac{\hat{\mbox{\boldmath$r$}}\;}{r^{2}}\,{\rm d}% \mbox{\boldmath$r$}\,{\rm d}\mbox{\boldmath$u$}\,.  (3) 
Note that \varphi and \mbox{\boldmath$A$}_{\bullet} have been defined such that they are O(1) quantities. Then the relative acceleration of a star (with respect to the MBH) is obtained by subtracting the acceleration of the MBH from the net gravitational force (per unit mass) on the star:
\frac{{\rm d}\mbox{\boldmath$u$}}{{\rm d}t}\;=\;\,\frac{GM_{\bullet}}{r^{2}}% \hat{\mbox{\boldmath$r$}}\;\;\varepsilon\frac{\mbox{$\partial$}\varphi}{\mbox% {$\partial$}\mbox{\boldmath$r$}}\;\;\varepsilon\mbox{\boldmath$A$}_{\bullet}\,.  (4) 
The stellar system may be regarded as collisionless over times shorter than the relaxation times for exchanging energy or angular momentum through stellar gravitational encounters. Then the time evolution of the DF is governed by the collisionless Boltzmann equation (CBE):
\displaystyle\frac{{\rm d}f}{{\rm d}t}  \displaystyle\;\equiv  \displaystyle\frac{\mbox{$\partial$}f}{\mbox{$\partial$}t}\;+\;\frac{{\rm d}% \mbox{\boldmath$r$}}{{\rm d}t}\mbox{\boldmath$\cdot$}\frac{\mbox{$\partial$}f}% {\mbox{$\partial$}\mbox{\boldmath$r$}}\;+\;\frac{{\rm d}\mbox{\boldmath$u$}}{{% \rm d}t}\mbox{\boldmath$\cdot$}\frac{\mbox{$\partial$}f}{\mbox{$\partial$}% \mbox{\boldmath$u$}}  
\displaystyle\;=  \displaystyle\frac{\mbox{$\partial$}f}{\mbox{$\partial$}t}\;+\;\mbox{\boldmath% $u$}\mbox{\boldmath$\cdot$}\frac{\mbox{$\partial$}f}{\mbox{$\partial$}\mbox{% \boldmath$r$}}\;+\;\left(\,\frac{GM_{\bullet}}{r^{2}}\hat{\mbox{\boldmath$r$}% }\;\;\varepsilon\frac{\mbox{$\partial$}\varphi}{\mbox{$\partial$}\mbox{% \boldmath$r$}}\;\;\varepsilon\mbox{\boldmath$A$}_{\bullet}\right)\mbox{% \boldmath$\cdot$}\frac{\mbox{$\partial$}f}{\mbox{$\partial$}\mbox{\boldmath$u$% }}\;=\;0\,. 
It can be verified that the CBE preserves the normalization of eqn.(1). Using Poisson Brackets (PBs), the CBE can be rewritten compactly as:
\frac{\mbox{$\partial$}f}{\mbox{$\partial$}t}\;+\;\left[\,f\,,\,H_{\rm org}\,% \right]_{(6)}\;=\;0\,,  (5) 
where
H_{\rm org}(\mbox{\boldmath$r$},\mbox{\boldmath$u$},t)\;=\;\frac{u^{2}}{2}\,% \,\frac{GM_{\bullet}}{r}\;+\;\varepsilon\varphi(\mbox{\boldmath$r$},t)\;+\;% \varepsilon\mbox{\boldmath$r$}\mbox{\boldmath$\cdot$}\mbox{\boldmath$A$}_{% \bullet}(t)  (6) 
is the Hamiltonian, and [\;,\;]_{(6)} is the 6–dim PB defined by:^{3}^{3}3We have used the subscript ‘6’ for the full PB because, in much of this paper, we will use a reduced PB (without subscript) — see eqns.(21) and (22) — that operates in a 4–dim subspace.
\left[\,f\,,\,H_{\rm org}\,\right]_{(6)}\;=\;\frac{\mbox{$\partial$}f}{\mbox{$% \partial$}\mbox{\boldmath$r$}}\mbox{\boldmath$\cdot$}\frac{\mbox{$\partial$}H_% {\rm org}}{\mbox{$\partial$}\mbox{\boldmath$u$}}\;\;\frac{\mbox{$\partial$}f}% {\mbox{$\partial$}\mbox{\boldmath$u$}}\mbox{\boldmath$\cdot$}\frac{\mbox{$% \partial$}H_{\rm org}}{\mbox{$\partial$}\mbox{\boldmath$r$}}\,.  (7) 
The Hamiltonian H_{\rm org} of eqn.(6) is the sum of an O(1) term \,E_{\rm k}=\left(u^{2}/2GM_{\bullet}/r\right)\, which is the Keplerian orbital energy, and an O(\varepsilon) term \,\varepsilon(\varphi+\mbox{\boldmath$r$}\mbox{\boldmath$\cdot$}\mbox{% \boldmath$A$}_{\bullet})\, which comes from the gravity of all the stars. Then the CBE can be written as:
\frac{\mbox{$\partial$}f}{\mbox{$\partial$}t}\;+\;\left[\,f\,,\,E_{\rm k}\,% \right]_{(6)}\;+\;\varepsilon\left[\,f\,,\,\varphi+\mbox{\boldmath$r$}\mbox{% \boldmath$\cdot$}\mbox{\boldmath$A$}_{\bullet}\,\right]_{(6)}\;=\;0\,,  (8) 
Since the Keplerian part is dominant, it is convenient to change the canonical variables from \{\mbox{\boldmath$r$},\mbox{\boldmath$u$}\} to new ones, the Delaunay variables in which the Kepler motion of every star appears simplest.
The Delaunay variables, {\cal D}\equiv\{I,L,L_{z};w,g,h\}, are a set of action–angle variables for the Kepler problem (Plummer, 1960; Murray & Dermott, 1999; Binney & Tremaine, 2008). The three actions are: I\,=\,\sqrt{GM_{\bullet}a\,}\,; L\,=\,I\sqrt{1e^{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. When expressed in terms of the {\cal D} variables, the Keplerian orbital energy assumes the simple and degenerate form E_{\rm k}(I)=1/2(GM_{\bullet}/I)^{2}. Hence for the Kepler problem all the Delaunay variables, excepting the orbital phase w, are constant in time; w itself advances at the rate
\Omega_{\rm k}(I)\;=\;\frac{{\rm d}E_{\rm k}}{{\rm d}I}\;=\;\frac{(GM_{\bullet% })^{2}}{I^{3}}\,.  (9) 
In terms of the {\cal D} variables, the 6–dim PB between two functions, \chi_{1}({\cal D}) and \chi_{2}({\cal D})\,, is:
\left[\,\chi_{1}\,,\,\chi_{2}\,\right]_{(6)}\;=\;\left(\frac{\mbox{$\partial$}% \chi_{1}}{\mbox{$\partial$}w}\frac{\mbox{$\partial$}\chi_{2}}{\mbox{$\partial$% }I}\frac{\mbox{$\partial$}\chi_{1}}{\mbox{$\partial$}I}\frac{\mbox{$\partial$% }\chi_{2}}{\mbox{$\partial$}w}\right)\,+\,\left(\frac{\mbox{$\partial$}\chi_{1% }}{\mbox{$\partial$}g}\frac{\mbox{$\partial$}\chi_{2}}{\mbox{$\partial$}L}% \frac{\mbox{$\partial$}\chi_{1}}{\mbox{$\partial$}L}\frac{\mbox{$\partial$}% \chi_{2}}{\mbox{$\partial$}g}\right)\,+\,\left(\frac{\mbox{$\partial$}\chi_{1}% }{\mbox{$\partial$}h}\frac{\mbox{$\partial$}\chi_{2}}{\mbox{$\partial$}L_{z}}% \frac{\mbox{$\partial$}\chi_{1}}{\mbox{$\partial$}L_{z}}\frac{\mbox{$\partial$% }\chi_{2}}{\mbox{$\partial$}h}\right)\,.  (10) 
The first PB in eqn.(8) is \,\left[\,f\,,\,E_{\rm k}\,\right]_{(6)}\,=\,\Omega_{\rm k}\left(\mbox{$% \partial$}f/\mbox{$\partial$}w\right)\,. In the second PB, \,\varphi and \,\mbox{\boldmath$A$}_{\bullet} have to be first expressed in terms of the {\cal D} variables. For this we need to rewrite the phase–space integrals of eqns.(2) and (3) in the {\cal D} variables. The DF is \,f({\cal D},t) and \,{\rm d}\mbox{\boldmath$r$}\,{\rm d}\mbox{\boldmath$u$}={\rm d}{\cal D}\equiv% {\rm d}I\,{\rm d}L\,{\rm d}L_{z}\,{\rm d}w\,{\rm d}g\,{\rm d}h\,. Lastly, the position vector \mbox{\boldmath$r$}({\cal D})=(x,y,z) is given by (Plummer, 1960; Murray & Dermott, 1999; Sambhus & Sridhar, 2000):
{\left(\begin{array}[]{c}x\\ \\ y\\ \\ z\end{array}\right)}={\left(\begin{array}[]{ccc}C_{g}C_{h}C_{i}S_{h}S_{g}&% \quadS_{g}C_{h}C_{i}S_{h}C_{g}&\quad S_{i}S_{h}\\ \\ C_{g}S_{h}+C_{i}C_{h}S_{g}&\quadS_{g}S_{h}+C_{i}C_{h}C_{g}&\quadS_{i}C_{h}\\ \\ S_{i}S_{g}&\quad S_{i}C_{g}&\quad C_{i}\end{array}\right)}{\left(\begin{array}% []{c}a(C_{\eta}e)\\ \\ a\sqrt{1e^{2}\,}\,S_{\eta}\\ \\ 0\end{array}\right)}  (11) 
where S and C are shorthand for sine and cosine of the angles given as subscript. Here a is the semi–major axis; \eta is the eccentric anomaly, related to the orbital phase through w=(\etae\sin\eta)\,; \,e=\sqrt{1L^{2}/I^{2}\,} is the eccentricity; and i is the inclination angle determined by \cos i=(L_{z}/L)\,. It is also useful to note that r=\sqrt{x^{2}+y^{2}+z^{2}}=a(1e\cos\eta)\,. Then eqns.(2) and (3) take the form:
\varphi({\cal D},t)\;=\;\,GM_{\bullet}\int\frac{\,f({\cal D}^{\prime},t)\,}{% \,\mbox{\boldmath$r$}\mbox{\boldmath$r$}^{\prime}\,}\,{\rm d}{\cal D}^{% \prime}\,,\qquad\qquad\mbox{\boldmath$A$}_{\bullet}(t)\;=\;GM_{\bullet}\int f(% {\cal D},t)\,\frac{\hat{\mbox{\boldmath$r$}}\;}{r^{2}}\,{\rm d}{\cal D}\,,  (12) 
where \mbox{\boldmath$r$}=\mbox{\boldmath$r$}({\cal D}) and \mbox{\boldmath$r$}^{\prime}=\mbox{\boldmath$r$}({\cal D}^{\prime}) are given in eqn.(11). Then the CBE of eqn.(8) in the {\cal D} variables is:
\frac{\mbox{$\partial$}f}{\mbox{$\partial$}t}\;+\;\Omega_{\rm k}\frac{\mbox{$% \partial$}f}{\mbox{$\partial$}w}\;+\;\varepsilon\left[\,f\,,\,\varphi\,+\,% \mbox{\boldmath$r$}\mbox{\boldmath$\cdot$}\mbox{\boldmath$A$}_{\bullet}\,% \right]_{(6)}\;=\;0\,.  (13) 
This form of the CBE makes apparent two natural time scales in the problem: (i) the Kepler orbital period \mbox{$T_{\rm kep}$}=2\pi/\Omega_{\rm k}\,, on which the phase w varies, and (ii) the longer secular time scale \mbox{$T_{\rm sec}$}=\varepsilon^{1}\mbox{$T_{\rm kep}$}\,. In the next section we orbit–average the CBE over the fast orbital phase w, and obtain a reduced CBE describing secular evolution.
3 Orbit–averaging via the Method of Multiple Scales
Since there are two well–separated time scales \,\left(\mbox{$T_{\rm kep}$}\;\,\mbox{and}\,\;\mbox{$T_{\rm sec}$}=\varepsilon% ^{1}\mbox{$T_{\rm kep}$}\gg\mbox{$T_{\rm kep}$}\right)\, in the problem, it is convenient to analyze the CBE of eqn.(13) using the method of multiple scales (MMS) — see e.g. Bender & Orszag (1978) for a general introduction. The MMS proceeds by assuming that all relevant quantities — like the DF, potential etc — are functions of two time variables, t and \tau=\varepsilon t, whose variations are treated as independent. This is an artifice whose efficacy is proved when, on obtaining an MMS solution, we set \tau=\varepsilon t\, to get the desired solution to the full CBE eqn.(13); thus we write f=f({\cal D},t,\tau)\,, \,\varphi=\varphi({\cal D},t,\tau)\,, \,\mbox{\boldmath$A$}_{\bullet}=\mbox{\boldmath$A$}_{\bullet}(t,\tau)\,. Then, replacing the time derivative \,(\mbox{$\partial$}/\mbox{$\partial$}t)\, in the full CBE of eqn.(13) by \,\left(\mbox{$\partial$}/\mbox{$\partial$}t+\varepsilon\,\mbox{$\partial$}/% \mbox{$\partial$}\tau\right)\, we obtain the MMS–CBE:
\frac{\mbox{$\partial$}f}{\mbox{$\partial$}t}\;+\;\varepsilon\frac{\mbox{$% \partial$}f}{\mbox{$\partial$}\tau}\;+\;\Omega_{\rm k}\frac{\mbox{$\partial$}f% }{\mbox{$\partial$}w}\;+\;\varepsilon\left[\,f\,,\,\varphi\,+\,\mbox{\boldmath% $r$}\mbox{\boldmath$\cdot$}\mbox{\boldmath$A$}_{\bullet}\,\right]_{(6)}\;=\;0\,.  (14) 
The next step in the MMS is to seek a solution to eqn.(14) that describes secular evolution plus small fluctuations. Then the dominant part of the DF must be independent of the fast time t, as well as the fast orbital phase w. Let us now gather the remaining 5 Delaunay variables and call them {\cal R}\equiv\{I,L,L_{z};g,h\}\,. A point in 5–dim {\cal R}–space represents a Keplerian ellipse of given semi–major axis, eccentricity, inclination, periapse angle and nodal longitude. Henceforth such an ellipse will be referred to as a Gaussian Ring or simply Ring. We want to imagine secular evolution of N\gg 1 stars as the evolution of N\gg 1 Gaussian Rings over times \sim\mbox{$T_{\rm sec}$}\,. Hence, to O(1), the DF we seek must be a function only of the variables ({\cal R},\tau). With these considerations we begin with the following ansatz for the DF, as a Fourier expansion in w :
f({\cal D},t,\tau)\;=\;\frac{1}{2\pi}\,F({\cal R},\tau)\;+\;\frac{\varepsilon}% {2\pi}\sum_{n=\infty}^{\infty}f_{n}({\cal R},t,\tau)\exp{[{\rm i}nw]}\;+\;O(% \varepsilon^{2})\,,  (15) 
where \,F\geq 0\, is the O(1) secular DF, and the O(\varepsilon) fluctuation is given by its Fourier–coefficients \,f_{n}=f_{n}^{\star}\,.
We will derive equations for F and all the f_{n} by substituting the expansion of eqn.(15) in the MMS–CBE eqn.(14), and ignoring O(\varepsilon^{2}) terms. To do this we take the following preparatory steps:

The quantities, \varphi({\cal D},t,\tau) and \mbox{\boldmath$A$}_{\bullet}(t,\tau), are needed only to O(1), which implies that we can set f\to(2\pi)^{1}\,F({\cal R},\tau) in eqns.(12).

To O(1) the self–consistent potential is independent of t, and is given by:
\varphi({\cal D},\tau)\;=\;\,GM_{\bullet}\int F({\cal R}^{\prime},\tau)\,{\rm d% }{\cal R}^{\prime}\oint\frac{{\rm d}w^{\prime}}{2\pi}\frac{1}{\,\mbox{% \boldmath$r$}\mbox{\boldmath$r$}^{\prime}\,}\,, (16) where on the left side we have dropped the dependence on the fast time. Fourier–expanding the right side in w, the potential can be written as:
\varphi({\cal D},\tau)\;=\;\Phi({\cal R},\tau)\;+\;\sum_{n\neq 0}\varphi_{n}({% \cal R},\tau)\exp{[{\rm i}nw]}\,, (17) where
\displaystyle\Phi({\cal R},\tau) \displaystyle\;= \displaystyle\int{\rm d}{\cal R}^{\prime}\,F({\cal R}^{\prime},\tau)\,\Psi({% \cal R},{\cal R}^{\prime})\,; (18a) \displaystyle\varphi_{n}({\cal R},\tau) \displaystyle\;= \displaystyle\int{\rm d}{\cal R}^{\prime}\,F({\cal R}^{\prime},\tau)\,\psi_{n}% ({\cal R},{\cal R}^{\prime})\,, (18b) with the Ring–Ring interaction potential functions given by
\displaystyle\Psi({\cal R},{\cal R}^{\prime}) \displaystyle\;= \displaystyleGM_{\bullet}\oint\oint\frac{{\rm d}w}{2\pi}\,\frac{{\rm d}w^{% \prime}}{2\pi}\,\frac{1}{\left\mbox{\boldmath$r$}\mbox{\boldmath$r$}^{\prime% }\right}\,, (19a) \displaystyle\psi_{n}({\cal R},{\cal R}^{\prime}) \displaystyle\;= \displaystyleGM_{\bullet}\oint\oint\frac{{\rm d}w}{2\pi}\,\frac{{\rm d}w^{% \prime}}{2\pi}\,\frac{\exp{[{\rm i}nw]}}{\left\mbox{\boldmath$r$}\mbox{% \boldmath$r$}^{\prime}\right}\,,\qquad n\neq 0\,. (19b) Note that \,\varphi_{n}=\varphi_{n}^{\star}\, and \,\psi_{n}=\psi_{n}^{\star}\,, because \Phi and \Psi are real quantities.

Using {\rm d}{\cal D}={\rm d}{\cal R}\,{\rm d}w\, in the second of eqns.(12), we get:
\mbox{\boldmath$A$}_{\bullet}\;=\;GM_{\bullet}\int F({\cal R}^{\prime},\tau)\,% {\rm d}{\cal R}^{\prime}\oint\frac{{\rm d}w}{2\pi}\,\frac{\hat{\mbox{\boldmath% $r$}}\;}{r^{2}}\;=\;{\bf 0}\,, (20) because \oint{\rm d}w\,\hat{\mbox{\boldmath$r$}}/r^{2}={\bf 0}\,, by the conservation of angular momentum along a Kepler orbit. So we can drop the \mbox{\boldmath$r$}\mbox{\boldmath$\cdot$}\mbox{\boldmath$A$}_{\bullet} term in the PB of eqn.(14).

The 6–dim PB of two functions, \chi_{1}({\cal D}) and \chi_{2}({\cal D}), can be expanded as:
\left[\,\chi_{1}\,,\,\chi_{2}\,\right]_{(6)}\;=\;\left(\frac{\mbox{$\partial$}% \chi_{1}}{\mbox{$\partial$}w}\frac{\mbox{$\partial$}\chi_{2}}{\mbox{$\partial$% }I}\frac{\mbox{$\partial$}\chi_{1}}{\mbox{$\partial$}I}\frac{\mbox{$\partial$% }\chi_{2}}{\mbox{$\partial$}w}\right)\;+\;\left[\,\chi_{1}\,,\,\chi_{2}\,% \right]\,, (21) where the PB on the right side (without subscript),
\left[\,\chi_{1}\,,\,\chi_{2}\,\right]\;=\;\left(\frac{\mbox{$\partial$}\chi_{% 1}}{\mbox{$\partial$}g}\frac{\mbox{$\partial$}\chi_{2}}{\mbox{$\partial$}L}% \frac{\mbox{$\partial$}\chi_{1}}{\mbox{$\partial$}L}\frac{\mbox{$\partial$}% \chi_{2}}{\mbox{$\partial$}g}\right)\,+\,\left(\frac{\mbox{$\partial$}\chi_{1}% }{\mbox{$\partial$}h}\frac{\mbox{$\partial$}\chi_{2}}{\mbox{$\partial$}L_{z}}% \frac{\mbox{$\partial$}\chi_{1}}{\mbox{$\partial$}L_{z}}\frac{\mbox{$\partial$% }\chi_{2}}{\mbox{$\partial$}h}\right)\,, (22) is the 4–dim PB, whose action is restricted to the 4 dimensional I=\mbox{constant}\; surfaces in the 5 dimensional Ring–space. In particular, for two functions \chi_{1}({\cal R}) and \chi_{2}({\cal R})\,, the 6–dim PB equals the 4–dim PB: \;\left[\,\chi_{1}\,,\,\chi_{2}\,\right]_{(6)}\,=\,\left[\,\chi_{1}\,,\,\chi_{% 2}\,\right]\,.
When eqns.(15) and (17) are substituted in the MMS–CBE eqn.(14), all the O(1) terms are zero. Collecting the O(\varepsilon) terms we get:
\displaystyle\frac{\mbox{$\partial$}f_{0}}{\mbox{$\partial$}t}\;+\;\frac{\mbox% {$\partial$}F}{\mbox{$\partial$}\tau}\;+\;\left[\,F\,,\,\Phi\,\right]\;=\;0\,;  (23a)  
\displaystyle\frac{\mbox{$\partial$}f_{n}}{\mbox{$\partial$}t}\;+\;{\rm i}n% \Omega_{\rm k}\,f_{n}\;=\;{\rm i}n\left(\frac{\mbox{$\partial$}F}{\mbox{$% \partial$}I}\varphi_{n}\;+\;\frac{{\rm i}}{n}\left[\,F\,,\,\varphi_{n}\,\right% ]\right)\,,\qquad\quad n\neq 0\,.  (23b) 
The second and third terms on the left side of eqn.(23a) are independent of the fast time t, so it is necessary that \mbox{$\partial$}f_{0}/\mbox{$\partial$}t also be independent of t. Then the general solution is f_{0}({\cal R},t,\tau)=\lambda_{0}({\cal R},\tau)+b({\cal R},\tau)\,t, where \lambda_{0} and b are arbitrary real functions of {\cal R} and \tau. However, f_{0} is an increasing function of t when b\neq 0. A cornerstone of the MMS is the requirement that physically admissible solutions must not grow on the fast time scale. This implies that b=0, and the MMS solution of eqn.(23a) is:
\displaystyle\frac{\mbox{$\partial$}F}{\mbox{$\partial$}\tau}\;+\;\left[\,F\,,% \,\Phi\,\right]\;=\;0\,,  (24a)  
\displaystyle f_{0}({\cal R},t,\tau)\;=\;\lambda_{0}({\cal R},\tau)\,.  (24b) 
Eqn.(24a) is the desired secular CBE for F({\cal R},\tau), and is studied in more detail in the next two sections. From eqn.(24b) we see that our O(\varepsilon) MMS theory only requires f_{0} to independent of t, while letting it be an arbitrary real function of {\cal R} and \tau. The right side of eqn.(23b) is independent of t, and it is straightforward to write down the general solution:
f_{n}({\cal R},t,\tau)\;=\;\frac{1}{\Omega_{\rm k}}\left(\frac{\mbox{$\partial% $}F}{\mbox{$\partial$}I}\varphi_{n}\;+\;\frac{{\rm i}}{n}\left[\,F\,,\,\varphi% _{n}\,\right]\right)\;+\;\lambda_{n}({\cal R},\tau)\,\exp{[{\rm i}n\Omega_{% \rm k}t]}\,,\qquad\quad n\neq 0\,,  (25) 
where the \lambda_{n} are arbitrary functions of {\cal R} and \tau, with \lambda_{n}=\lambda_{n}^{\star}\,. Substituting eqns.(24b) and (25) in eqn.(15), the general solution to the MMS–CBE eqn.(14) can be written as:
\displaystyle f({\cal D},t,\tau)  \displaystyle\;=  \displaystyle\frac{1}{2\pi}F({\cal R},\tau)\;+\;\frac{\varepsilon}{2\pi\Omega_% {\rm k}}\,\sum_{n\neq 0}\left(\frac{\mbox{$\partial$}F}{\mbox{$\partial$}I}% \varphi_{n}\;+\;\frac{{\rm i}}{n}\left[\,F\,,\,\varphi_{n}\,\right]\right)\exp% {[{\rm i}nw]}  (26)  
\displaystyle\qquad\;+\;\frac{\varepsilon}{2\pi}\,\sum_{n=\infty}^{\infty}% \lambda_{n}({\cal R},\tau)\exp{[{\rm i}n(w\Omega_{\rm k}t)]}\;+\;O(% \varepsilon^{2})\,. 
Define two t–independent functions, \Delta F({\cal D},\tau) and \Lambda({\cal R},w,\tau), by the following Fourier series:
\displaystyle\Delta F({\cal D},\tau)  \displaystyle\;=  \displaystyle\frac{1}{\Omega_{\rm k}}\,\sum_{n\neq 0}\left(\frac{\mbox{$% \partial$}F}{\mbox{$\partial$}I}\varphi_{n}\;+\;\frac{{\rm i}}{n}\left[\,F\,,% \,\varphi_{n}\,\right]\right)\exp{[{\rm i}nw]}\,,  (27a)  
\displaystyle\Lambda({\cal R},w,\tau)  \displaystyle\;=  \displaystyle\sum_{n=\infty}^{\infty}\,\lambda_{n}({\cal R},\tau)\exp{[{\rm i% }nw]}\,.  (27b) 
The Fourier coefficients on the right side of eqn.(27a) are completely determined by F\,; they are well–defined because n and \Omega_{\rm k} are both non zero. We do not attempt to prove convergence of the Fourier series, and assume that this can be checked in any particular case. In eqn.(27b) the \lambda_{n} are arbitrary functions of {\cal R} and \tau; they can be chosen such that the Fourier series converges and \Lambda({\cal R},w,\tau) is a welldefined function of {\cal D} and \tau. Having obtained the MMS solution, we go back to using a single time variable t, because we want a solution to the original problem of the full CBE eqn.(13). In other words, \tau is no longer regarded as an independent variable, but is simply \tau=\varepsilon t, which is a function of the true time variable t, and the small parameter \varepsilon. Setting \tau=\varepsilon t in eqn.(26), we obtain
f({\cal D},t)\;=\;\frac{1}{2\pi}F({\cal R},\varepsilon t)\;+\;\frac{% \varepsilon}{2\pi}\Delta F({\cal D},\varepsilon t)\;+\;\frac{\varepsilon}{2\pi% }\Lambda\left({\cal R},\,w\Omega_{\rm k}t,\,\varepsilon t\right)\;+\;O(% \varepsilon^{2})\,.  (28) 
which is a solution of the full CBE eqn.(13), accurate to O(\varepsilon) over times \sim\mbox{$T_{\rm sec}$}; this can be verified by direct substitution of eqn.(28) in eqn.(13).
The DF of eqn.(28) is the main result of this section. It is the sum of three terms: the first is O(1) and proportional to F\,; the two other terms are O(\varepsilon), and can be thought of as the small corrections to F. Some salient points to be noted are:

The secular DF F is O(1), independent of the orbital phase w, and evolves over times \mbox{$T_{\rm sec}$}\,, as determined by the CBE eqn.(24a).

The first O(\varepsilon) correction is proportional to \Delta F which is a fluctuating function of the orbital phase w with zero mean. It is completely determined by F and, like F, evolves over times \mbox{$T_{\rm sec}$}\,.

The second O(\varepsilon) correction is proportional to \Lambda which is a function of \,{\cal R}, \,(w\Omega_{\rm k}t)\, and \,\varepsilon t\,. Its functional form is undetermined in our O(\varepsilon) MMS theory. However, it is important to note that the fast time t occurs only in the form (w\Omega_{\rm k}t), which contributes only in a periodic manner, with no growth or decay.
The derivation of the full DF of eqn.(28) given above makes it natural to think of the orbits followed by individual stars, to O(1) accuracy, as slowly deforming Gaussian Rings. However, it is also useful to see this directly from stellar orbital dynamics. This is done in Appendix A using the standard apparatus of secular perturbation theory. There it is proved: Every stellar orbit can be described to O(1) accuracy as a Gaussian Ring whose dynamics is governed by the Hamiltonian, \Phi({\cal R},\tau). Superimposed on this are O(\varepsilon) oscillations on the fast time T_{\rm kep}, so we may think of a stellar orbit as a slowly evolving, ‘noisy’ Gaussian Ring.
MMS analysis is very useful because it deals directly with the DF — bypassing consideration of individual stellar orbits — and offers the following: for a given secular evolution described by the O(1) DF F({\cal R},\tau) in 5–dim {\cal R} space, there is a full DF f({\cal D},t) in the 6–dim {\cal D} space for which the O(\varepsilon) noise contributed by all the individual stellar orbits does not act cooperatively to give rise to growth on the fast time scale T_{\rm kep}. Hence the fluctuations remain O(\varepsilon) over the long secular time \mbox{$T_{\rm sec}$}\,, and the O(1) dynamics described by F is well–defined. Thus MMS provides a consistent framework for secular dynamics, but it would be well to note that this very consistency has been achieved at a certain cost. MMS filters out all O(1) changes over fast orbital times T_{\rm kep}, and hence the resulting description is insensitive to either fast external perturbations or instabilities. It is generally expected that these fast phenomena would have a period of growth followed by saturation. Further evolution — including any instability — unfolds over times \sim\mbox{$T_{\rm sec}$}. This is the regime that is accurately described by the MMS analysis of this Section.^{4}^{4}4An example of a fast instability is the well–known axisymmetric instability (Toomre, 1964) applied to a Keplerian disc. All but the most unrealistically cold Keplerian discs would be stable by the Toomre criterion. But, as a matter of principle, we can indeed consider a Keplerian disc cold enough to be unstable. Then a secular description would make sense only after the instability has grown and saturated over times T_{\rm kep}, leaving the disc in a warmer state that is stable to all fast modes.
4 Collisionless Boltzmann equation for Gaussian Rings
Henceforth we will study the O(1) slow evolution described by the secular DF F({\cal R},\tau). To this order the normalization of the full DF eqn.(1) implies that
\int F({\cal R},\tau)\,{\rm d}{\cal R}\;=\;1,  (29) 
so we may regard F({\cal R},\tau) as a probability distribution function in 5–dim Ring space. We shall refer to F as the Ring distribution function. The time evolution of F is described by the secular CBE eqn.(24a), \;(\mbox{$\partial$}F/\mbox{$\partial$}\tau)+[F,\Phi]=0\,. We recall that \Phi({\cal R},\tau) is the (scaled) orbit–averaged self–gravitational potential of the stars which is determined self–consistently by F itself:
\Phi({\cal R},\tau)\;=\;\int{\rm d}{\cal R}^{\prime}\,F({\cal R}^{\prime},\tau% )\,\Psi({\cal R},{\cal R}^{\prime})\,,  (30) 
where
\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\mbox{\boldmath$r$}\mbox% {\boldmath$r$}^{\prime}\right}  (31) 
is the (scaled) Ring–Ring interaction potential energy. Note that \Phi is independent of \varepsilon, so the secular evolution of eqn.(24a) has the property of mass–scale invariance: i.e. the behaviour of F as a function of {\cal R} and \tau is independent of \varepsilon when \varepsilon\ll 1\,. This means that, for given MBH mass M_{\bullet}, the evolution of the stellar system appears independent of its mass M, so long as M\ll M_{\bullet}\, and the time variable used is \tau=\varepsilon t. This invariance holds when every star is acted upon by the Keplerian potentials of the MBH and all the other stars. We now add two corrections to \Phi\,; these are the orbit–averaged effects of (a) relativistic corrections and (b) external perturbations.
(a) Relativistic Corrections to the gravity of the MBH introduce post–Newtonian (PN) corrections, the most important ones being the 1 PN Schwarzschild precession of the apses, and the 1.5 PN Lense–Thirring precession of the apses and nodes for a spinning MBH. Both effects are described by the relativistic (secular) Hamiltonian \,\varepsilon H^{\rm rel}\,:
\displaystyle H^{\rm rel}(I,L,L_{z})  \displaystyle\;=  \displaystyle H^{PN1}(I,L)\;+\;H^{PN1.5}(I,L,L_{z})\,,  
\displaystyle H^{PN1}(I,L)  \displaystyle\;=  \displaystyle\,B_{1}\,\frac{1}{I^{3}L}\qquad\mbox{with}\qquad B_{1}\;=\;\frac% {3(GM_{\bullet})^{4}}{c^{2}}\,\frac{M_{\bullet}}{M}\,,  
\displaystyle H^{PN1.5}(I,L,L_{z})  \displaystyle\;=  \displaystyle B_{1.5}\,\frac{L_{z}}{I^{3}L^{3}}\,,\qquad\mbox{with}\qquad B_{1% .5}\;=\;\frac{2(GM_{\bullet})^{5}}{c^{3}}\,\frac{M_{\bullet}}{M}\,\sigma\,,  (32) 
We have assumed that the spin of the MBH points along the z–axis; \,0\leq\sigma\leq 1\, is the spin parameter. Note that the definitions of the constants B_{1} and B_{1.5}\, include the factor (M_{\bullet}/M)=\varepsilon^{1}.
(b) External Perturbations on the stars as well as the MBH — due to nuclear density cusps and/or slowly moving distant masses — can also be included in the secular description. Let \,\varepsilon\Phi^{\rm ext}({\cal R},\tau)\, be the orbit–averaged contribution and \,\varepsilon\mbox{\boldmath$A$}^{\rm ext}(\tau) be the acceleration of the MBH. These will contribute an orbit–averaged tidal potential \,\varepsilon\Phi^{\rm tid}({\cal R},\tau)\, to the secular Hamiltonian for a star, where
\Phi^{\rm tid}({\cal R},\tau)\;=\;\Phi^{\rm ext}({\cal R},\tau)\;+\;\mbox{% \boldmath$X$}({\cal R})\cdot\mbox{\boldmath$A$}^{\rm ext}_{\bullet}(\tau)\,,  (33) 
with
\mbox{\boldmath$X$}({\cal R})\;=\;\oint\frac{{\rm d}w}{2\pi}\,\mbox{\boldmath$% r$}({\cal D})  (34) 
equal to the orbit–averaged position vector. Note that, similar to the case of relativistic perturbations, \,\Phi^{\rm tid}\, has a factor \varepsilon^{1} implicit in its definition.
Putting the pieces together, the general Hamiltonian for a Gaussian Ring is:
H({\cal R},\tau)\;=\;\Phi({\cal R},\tau)\;+\;H^{\rm rel}(I,L,L_{z})\;+\;\Phi^{% \rm tid}({\cal R},\tau)\,,  (35) 
with the Ring equations of motion:
\displaystyle I  \displaystyle\;=  \displaystyle\sqrt{GM_{\bullet}a}\;=\;\mbox{constant}\,,  
\displaystyle\frac{{\rm d}L}{{\rm d}\tau}  \displaystyle\;=  \displaystyle\,\frac{\mbox{$\partial$}H}{\mbox{$\partial$}g}\,,\qquad\frac{{% \rm d}g}{{\rm d}\tau}\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L}\,;% \qquad\frac{{\rm d}L_{z}}{{\rm d}\tau}\;=\;\,\frac{\mbox{$\partial$}H}{\mbox{% $\partial$}h}\,,\qquad\frac{{\rm d}h}{{\rm d}\tau}\;=\;\frac{\mbox{$\partial$}% H}{\mbox{$\partial$}L_{z}}\,.  (36) 
A secular integral of motion is a function \,{\cal I}({\cal R},\tau) that is constant along the orbit of a Ring:^{5}^{5}5Henceforth we drop the word ‘secular’ when referring to an integral of motion.
\frac{{\rm d}}{{\rm d}\tau}{\cal I}\left({\cal R}(\tau),\tau\right)\;=\;0\,.  (37) 
It is a feature of secular dynamics that I=\sqrt{GM_{\bullet}a}\, is an integral of motion for quite general and time–dependent H({\cal R},\tau). When the Hamiltonian is time–independent, it provides a second integral H({\cal R}).
The equation governing the evolution of the distribution of Gaussian Rings is obtained by replacing \Phi by H in eqn.(24a):
\frac{\mbox{$\partial$}F}{\mbox{$\partial$}\tau}\;+\;\left[\,F\,,\,H\,\right]% \;=\;0\,.  (38) 
The time evolution is collisionless: using eqns.(36) and (38) we can verify that
\displaystyle\frac{{\rm d}F}{{\rm d}\tau}  \displaystyle\;\equiv  \displaystyle\frac{\mbox{$\partial$}F}{\mbox{$\partial$}t}\;+\;\frac{{\rm d}L}% {{\rm d}\tau}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}L}\;+\;\frac{{\rm d}g}% {{\rm d}\tau}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}g}\;+\;\frac{{\rm d}L_% {z}}{{\rm d}\tau}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}L_{z}}\;+\;\frac{{% \rm d}h}{{\rm d}\tau}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}h}  (39)  
\displaystyle\;=  \displaystyle\frac{\mbox{$\partial$}F}{\mbox{$\partial$}\tau}\;\;\frac{\mbox{% $\partial$}H}{\mbox{$\partial$}g}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}L}% \;+\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L}\frac{\mbox{$\partial$}F}{% \mbox{$\partial$}g}\;\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}h}\frac{% \mbox{$\partial$}F}{\mbox{$\partial$}L_{z}}\;+\;\frac{\mbox{$\partial$}H}{% \mbox{$\partial$}L_{z}}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}h}  
\displaystyle\;=  \displaystyle\frac{\mbox{$\partial$}F}{\mbox{$\partial$}\tau}\;+\;\left[\,F\,,% \,H\,\right]\;=\;0\,, 
Therefore eqn.(38) describes a system of Gaussian Rings evolving collisionlessly under the combined actions of self–gravity, relativistic corrections to the MBH’s Newtonian gravity and slow external perturbations. It will be referred to as the Ring CBE.
The Ring Hamiltonian H of eqn.(35) is in general not independent of \varepsilon, because both H^{\rm rel} and \Phi^{\rm tid} are proportional to \varepsilon^{1}. Therefore the property of mass–scale invariance of the Ring CBE is in general broken. Another important general property, valid for quite general F({\cal R},\tau) and H({\cal R},\tau), is this: The probability for a Ring to be in (I,\,I+{\rm d}I) is a conserved quantity. In other words the PDF in 1–dim I–space, defined by
P(I)\;=\;\int\,dL\,dL_{z}\,dg\,dh\,F(I,\,L,\,L_{z},\,g,\,h,\,\tau)\,,  (40) 
is independent of \tau, as can be verified directly using the Ring CBE eqn.(38). The physical reason for this is that the semi–major axis of every Ring is conserved, so the function P(I) must stay frozen as F({\cal R},\tau) evolves over secular times. Therefore secular evolution is driven by gravitational torques between Gaussian Rings, with no exchange of Keplerian energies.
4.1 Collisionless Boltzmann equation for Zero–thickness Flat Discs
It is useful to study the idealized case when all the Rings are confined to the xy plane. In these zero–thickness discs, the angular momentum of every Ring points along \pm\hat{z}, so that a Ring is specified by its semi–major axis, angular momentum and apsidal longitude. The Ring phase space is now 3–dim: we write {\cal R}=\{I,L,g\}, where L stands for the angular momentum which can now take both positive and negative values \left(I\leq L\leq I\,\right), and g is the longitude of the periapse.^{6}^{6}6This is the convention for the symbols (L,g) we will always use when dealing with Ring dynamics that is confined to the xy plane. Rings with L>0 will be referred to as prograde and Rings with L<0 will be referred to as retrograde.
The general time–dependent Ring DF is F(I,L,g,\tau). The (scaled) orbit–averaged self–gravitational potential, \Phi(I,L,g,\tau), is given in terms of F and \Psi({\cal R},{\cal R}^{\prime}) by eqns.(30) and (31). To calculate \Psi({\cal R},{\cal R}^{\prime}), we express \mbox{\boldmath$r$}=(x,y) in terms of the 2–dim Delaunay variables:
\displaystyle\qquad\qquad x  \displaystyle\;=  \displaystyle a(\cos\etae)\cos g\,\;\mp\;\,a\sqrt{1e^{2}\,}\sin\eta\,\sin g\,,  
\displaystyle y  \displaystyle\;=  \displaystyle a(\cos\etae)\sin g\,\;\pm\;\,a\sqrt{1e^{2}\,}\sin\eta\,\cos g\,,  (41) 
where the upper(lower) signs in the second terms on the right hand side correspond to prograde(retrograde) orbits. As earlier, \eta is related to the mean anomaly by w=\etae\sin\eta. Using corresponding expressions for \mbox{\boldmath$r$}^{\prime}=(x^{\prime},y^{\prime}), we see that \left\mbox{\boldmath$r$}\mbox{\boldmath$r$}^{\prime}\right^{1} is an explicitly known function of the Delaunay variables, {\cal D} and {\cal D}^{\prime}, and can be averaged over w and w^{\prime}. Then the secular Hamiltonian governing Ring dynamics is:
H(I,L,g,\tau)\;=\;\Phi(I,L,g,\tau)\;+\;H^{\rm rel}(I,L)\;+\;\Phi^{\rm tid}(I,L% ,g,\tau)\,,  (42) 
where the relativistic correction,
H^{\rm rel}(I,L)\;=\;B_{1}\frac{1}{I^{3}L}\;+\;B_{1.5}\frac{{\rm Sgn}(L)}{I% ^{3}L^{2}}\,,  (43) 
takes a slightly different form from eqn.(32), because we have allowed L to take both positive and negative values for zero–thickness flat discs. The Ring equations of motion are:
I\;=\;\sqrt{GM_{\bullet}a}\;=\;\mbox{constant}\,,\qquad\quad\frac{{\rm d}L}{{% \rm d}\tau}\;=\;\,\frac{\mbox{$\partial$}H}{\mbox{$\partial$}g}\,,\qquad\quad% \frac{{\rm d}g}{{\rm d}\tau}\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L}\,.  (44) 
The DF satisfies the Ring CBE,
\frac{{\rm d}F}{{\rm d}\tau}\;=\;\frac{\mbox{$\partial$}F}{\mbox{$\partial$}% \tau}\;\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}g}\frac{\mbox{$\partial$}% F}{\mbox{$\partial$}L}\;+\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L}\frac{% \mbox{$\partial$}F}{\mbox{$\partial$}g}\;=\;\frac{\mbox{$\partial$}F}{\mbox{$% \partial$}\tau}\;+\;\left[\,F\,,\,H\,\right]\;=\;0\,,  (45) 
where the Poisson Bracket is now 2–dim, because it involves only the canonically conjugate pair (L,g). As earlier, the probability for a Ring to be in (I,\,I+{\rm d}I) is a conserved quantity: i.e. the PDF in 1–dim I–space, defined by
P(I)\;=\;\int\,dL\,dg\,\,F(I,\,L,\,g,\,\tau)\,,  (46) 
is independent of \tau, as can be verified directly using the Ring CBE eqn.(45).
5 Secular Collisionless Equilibria
Secular collisionless equilibria are described by DFs F that are independent of \tau; these DFs will also be referred to as stationary. From eqn.(38) a stationary DF must satisfy \left[F,H\right]=0\,. These can be constructed using a secular version of Jeans theorem.
Secular Jeans theorem: Any stationary solution of the Ring CBE eqn.(38) is a function of {\cal R} only through the time–independent integrals of motion of H({\cal R}), and any function of the time–independent integrals of H({\cal R}) is a stationary solution of eqn.(38).
Proof: This is similar to the one given in Binney & Tremaine (2008): If F is a stationary solution of eqn.(38), then it is also an integral of motion; so the first part of the theorem is proved. Conversely, if \;{\cal I}_{1}({\cal R}),\ldots,{\cal I}_{n}({\cal R})\, be n time–independent integrals, then any F({\cal I}_{1},\ldots,{\cal I}_{n}) satisfies the Ring CBE:
\frac{{\rm d}}{{\rm d}\tau}F\left({\cal I}_{1}({\cal R}),\ldots,{\cal I}_{n}({% \cal R})\right)\;=\;\sum_{m=1}^{n}\frac{\mbox{$\partial$}F}{\mbox{$\partial$}{% \cal I}_{m}}\frac{{\rm d}{\cal I}_{m}}{{\rm d}\tau}\;=\;0\,,  (47) 
because ({\rm d}{\cal I}_{m}/{\rm d}\tau)=0 for all m=1,\ldots,n\,.\quad\square
In practice a stronger version of the Secular Jeans theorem is used, by requiring a stationary DF to be a function of the isolating integrals of motion. In contrast with non–secular dynamics two isolating integrals, I and H({\cal R}), are always available in secular dynamics. Additional symmetries can promote more isolating integrals. Hence a richer variety of stationary DFs can be constructed, with greater ease than in non secular dynamics. Below we discuss briefly the simplest cases.
5.1 Spherical Equilibria
For a spherically symmetric stellar system around a MBH, the orbit–averaged self–gravitational potential \Phi is a function only of I and L. An external potential can be included if it is also spherically symmetric with respect to the MBH; this causes no acceleration of the MBH, so \Phi^{\rm tid}=\Phi^{\rm ext} are also functions only of I and L. Therefore the Ring Hamiltonian of eqn.(35),
H(I,L,L_{z})\;=\;\Phi(I,L)\;+\;H^{\rm rel}(I,L,L_{z})\;+\;\Phi^{\rm ext}(I,L)\,,  (48) 
is independent of the angle variables. Using eqn.(36) the orbit of a Ring is given by:
\displaystyle I  \displaystyle\;=  \displaystyle\mbox{constant}\,,\qquad L\;=\;\mbox{constant}\,,\qquad L_{z}\;=% \;\mbox{constant}\,,  
\displaystyle\frac{{\rm d}g}{{\rm d}\tau}  \displaystyle\;\equiv  \displaystyle\Omega(I,L,L_{z})\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}% L}\;;\qquad\mbox{apse precession frequency}\,,  
\displaystyle\frac{{\rm d}h}{{\rm d}\tau}  \displaystyle\;\equiv  \displaystyle\nu(I,L)\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L_{z}}\;=% \;\frac{B_{1.5}}{(IL)^{3}}\;;\qquad\mbox{node precession frequency}\,.  (49) 
Therefore, every Ring in a spherical system has constant semi–major axis, eccentricity and inclination; its periapse precesses at a constant rate \Omega(I,L,L_{z}); its ascending node precesses at a constant, inclination–independent rate \nu(I,L), due solely to 1.5 PN relativistic effects for a spinning MBH. The Ring orbital dynamics is integrable, with the Delaunay actions (I,L,L_{z}) serving as three isolating integrals of motion.
The Secular Jeans theorem implies that a general stationary spherical DF can be a function of I, L and L_{z}\,. We will restrict attention to DFs which have complete spherical symmetry which these are of the form F(I,L). It should be noted that the stellar distribution in phase space is given explicitly to O(1), without having to solve an integral equation. Indeed, the functional dependence on the six coordinates \{\mbox{\boldmath$r$},\mbox{\boldmath$u$}\} is obtained by setting
I\;=\;GM_{\bullet}\left[\frac{2GM_{\bullet}}{r}\,\,u^{2}\right]^{1/2}\,,% \qquad\qquad L\;=\;\mbox{\boldmath$r$}\mbox{\boldmath$\times$}\mbox{\boldmath% $u$}\,,  (50) 
in the DF F(I,L). Then, observationally interesting quantities — such as the radial profiles of the surface density, the velocity dispersion tensor and the line–of–sight velocity distribution (LOSVD) — can be calculated by direct integration of the DF. It is readily seen that the DFs F(I,L) describe systems with zero streaming motions. The velocity dispersion tensor is in general anisotropic.
5.2 Axisymmetric Equilibria
When the stellar system is axisymmetric about the MBH with symmetry axis along \hat{z}, the orbit–averaged self–gravitational potential \Phi is a function of (I,\,L,\,L_{z},\,g). An external potential can be included if it is also axisymmetric with respect to the MBH with the same symmetry axis; this causes no acceleration of the MBH, so \Phi^{\rm tid}=\Phi^{\rm ext} are also functions only of (I,\,L,\,L_{z},\,g). Therefore the Ring Hamiltonian of eqn.(35) simplifies to
H(I,L,L_{z},g)\;=\;\Phi(I,L,L_{z},g)\;+\;H^{\rm rel}(I,L,L_{z})\;+\;\Phi^{\rm ext% }(I,L,L_{z},g)\,.  (51) 
Using eqn.(36), the orbit of a Ring is given by
I\;=\;\mbox{constant}\,,\qquad L_{z}\;=\;\mbox{constant}\,,\qquad\frac{{\rm d}% g}{{\rm d}\tau}\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L}\,,\qquad% \frac{{\rm d}h}{{\rm d}\tau}\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial$}L_% {z}}\,.  (52) 
As a Ring evolves with constant I and L_{z}, its eccentricity, inclination, apse and node can all vary with time. From the Secular Jeans theorem, a general stationary axisymmetric DF must be of the form F(I,L_{z},H). Since H=\Phi+H^{\rm rel}+\Phi^{\rm ext} is a function of (I,L,L_{z},g), so is the DF. However, to express the DF explicitly as a function of the phase space coordinates, we need to determine \Phi(I,L,L_{z},g) by solving the self–consistent problem of eqns.(30) and (51):
\Phi(I,L,L_{z},g)\;=\;\int{\rm d}I^{\prime}\,{\rm d}L^{\prime}\,{\rm d}L_{z}^{% \prime}\,{\rm d}g^{\prime}\,\,F(I^{\prime},L^{\prime}_{z},H^{\prime})\oint{\rm d% }h^{\prime}\,\Psi({\cal R},{\cal R}^{\prime})\,,  (53) 
where H^{\prime}=H(I^{\prime},L^{\prime},L_{z}^{\prime},g^{\prime}).^{7}^{7}7The left side of eqn.(53) indicates that \Phi is independent of h, as must be the case for any axisymmetric system. We can check that the right side of eqn.(53) is also independent of h: the nodal longitudes occur in \Psi({\cal R},{\cal R}^{\prime}) only in the combination (hh^{\prime}). So, in the integral of \Phi over h^{\prime}, when we change the integration variable from h^{\prime} to (h^{\prime}h), we get a quantity that is independent of h.
The problem simplifies considerably for DFs that are of the form F(I,L_{z}). By replacing I and L_{z} by
I\;=\;GM_{\bullet}\left[\frac{2GM_{\bullet}}{r}\,\,u^{2}\right]^{1/2}\,,% \qquad\qquad L_{z}\;=\;\hat{z}\mbox{\boldmath$\cdot$}\left(\mbox{\boldmath$r$}% \mbox{\boldmath$\times$}\mbox{\boldmath$u$}\right)\,,  (54) 
we see that, just like in the spherical case, the DF is known explicitly to O(1) as a function of \{\mbox{\boldmath$r$},\mbox{\boldmath$u$}\}, without having to solve an integral equation. Then observationally interesting quantities — such as the sky–maps of the surface density, the velocity dispersion tensor and the line–of–sight velocity distribution (LOSVD) — can be calculated by direct integration of the DF.
(a) Zero–thickness Flat Discs: We recall from § 4.1 that the Ring space is now 3–dim; {\cal R}=(I,L,g), where I\leq L\leq I is the angular momentum and g is the apsidal longitude. For an axisymmetric disc around a MBH, the (scaled) orbit–averaged self–gravitational potential \Phi is a function only of I and L. An external potential can be included if it is also axisymmetric with respect to the MBH; this causes no acceleration of the MBH, so \,\Phi^{\rm tid}=\Phi^{\rm ext} are also functions only of I and L. Therefore the Ring Hamiltonian of eqn.(42),
H(I,L)\;=\;\Phi(I,L)\;+\;H^{\rm rel}(I,L)\;+\;\Phi^{\rm ext}(I,L)\,,  (55) 
is independent of the angle variables. Using eqn.(44) the orbit of a Ring is given by:
I\;=\;\mbox{constant}\,,\qquad L\;=\;\mbox{constant}\,,\qquad\frac{{\rm d}g}{{% \rm d}\tau}\;\equiv\;\Omega(I,L)\;=\;\frac{\mbox{$\partial$}H}{\mbox{$\partial% $}L}\,.  (56) 
Therefore, every Ring in an axisymmetric system has constant semi–major axis and eccentricity, with its apsidal longitude precessing at the constant rate \Omega(I,L). The Ring orbital dynamics is integrable, with the Delaunay actions (I,L) serving as two isolating integrals of motion. The Secular Jeans theorem implies that a general stationary axisymmetric DF has the form F(I,L). The stellar distribution in phase space is given explicitly to O(1), without having to solve a self–consistent problem. Indeed, the functional dependence on the four coordinates \{x,y,u_{x},u_{y}\} is obtained by setting
I\;=\;GM_{\bullet}\left[\frac{2GM_{\bullet}}{r}\,\,u^{2}\right]^{1/2}\,,% \qquad\qquad L\;=\;xu_{y}\,\,yu_{x}\,,  (57) 
in the DF F(I,L). Similar to 3–dim axisymmetric systems, F(I,L_{z}) discussed above, observationally interesting quantities — such as the velocity dispersion tensor or the LOSVD — can be calculated by direct integration of the DF. Since L can take both positive and negative values, a general F(I,L) describes discs with or without rotation, and generally anisotropic velocity dispersions. These DFs are relevant to studies such as Touma & Tremaine (2014), who solved for maximum entropy equilibria of self–gravitating Keplerian disks (in the case where all particles have the same semi–major axis). When F(I,L) is an even function of L, local streaming motions are absent and the disc is non–rotating with anisotropic velocity dispersions. A special case consists of DFs of the form F(I), which are non–rotating and have isotropic velocity dispersions.
6 Perturbations and Stability of Secular Equilibria
If F_{0}({\cal R}) be a stationary DF, the Hamiltonian governing the orbits of Gaussian Rings in the unperturbed system is H_{0}({\cal R})=\Phi_{0}({\cal R})+H^{\rm rel}(I,L,L_{z})+\Phi_{0}^{\rm tid}({% \cal R}), and we must have [F_{0}\,,\,H_{0}]=0. Let F_{1}({\cal R},\tau) be an infinitesimal perturbation about F_{0}. This gives rise to an infinitesimal perturbation of the self–gravitational potential,
\Phi_{1}({\cal R},\tau)\;=\;\int{\rm d}{\cal R}^{\prime}\,F_{1}({\cal R}^{% \prime},\tau)\,\Psi({\cal R},{\cal R}^{\prime})\,.  (58) 
Let the external potential also be allowed to change infinitesimally, so it contributes an additional \Phi_{1}^{\rm tid}({\cal R},\tau)\,. Then the infinitesimal change in the Ring Hamiltonian is
H_{1}({\cal R},\tau)\;=\;\Phi_{1}({\cal R},\tau)\;+\;\Phi_{1}^{\rm tid}({\cal R% },\tau)\,.  (59) 
Substituting F=F_{0}+F_{1} and H=H_{0}+H_{1} in the Ring CBE eqn.(38), and keeping only the linear terms, gives the linearized collisionless Boltzmann equation (LCBE):
\frac{\mbox{$\partial$}F_{1}}{\mbox{$\partial$}\tau}\;+\;\left[\,F_{1}\,,\,H_{% 0}\,\right]\;+\;\left[\,F_{0}\,,\,H_{1}\,\right]\;=\;0\,.  (60) 
Eqns.(58)—(60) define the linear integral equation for the perturbation F_{1}({\cal R},\tau). When \Phi_{1}^{\rm tid}=0, the evolution of an infinitesimal perturbation is determined by only its self–gravity \Phi_{1}({\cal R},\tau) — and the unperturbed flows due to H_{0}({\cal R}) in which we have allowed for an external contribution \Phi_{0}^{\rm tid}({\cal R}). The unperturbed DF F_{0}({\cal R}) is said to be secularly stable if there is no growing solution F_{1}. Below we discuss two problems:

Long–lived, warp–like perturbations of spherically symmetric DFs of the form F_{0}(I,L).

Axisymmetric, zero–thickness, flat discs with DFs of the form F_{0}(I,L).
6.1 Perturbations of Spherical Equilibria
For a completely spherical unperturbed system the DF is F_{0}(I,L), with the unperturbed Ring Hamiltonian given by eqn.(48):
H_{0}(I,L,L_{z})\;=\;\Phi_{0}(I,L)\;+\;H^{\rm rel}(I,L,L_{z})\;+\;\Phi_{0}^{% \rm ext}(I,L)\,.  (61) 
From eqn.(49), the apse precession frequency is \Omega(I,L,L_{z})=\mbox{$\partial$}H_{0}/\mbox{$\partial$}L; and the inclination–independent node precession frequency, \nu(I,L)=B_{1.5}(IL)^{3}\,, is due only to the 1.5 PN relativistic correction for a spinning MBH. Then the LCBE of eqn.(60) can be written as:
\frac{\mbox{$\partial$}F_{1}}{\mbox{$\partial$}\tau}\;+\;\Omega\frac{\mbox{$% \partial$}F_{1}}{\mbox{$\partial$}g}\;+\;\nu\frac{\mbox{$\partial$}F_{1}}{% \mbox{$\partial$}h}\;=\;\frac{\mbox{$\partial$}F_{0}}{\mbox{$\partial$}L}\frac% {\mbox{$\partial$}H_{1}}{\mbox{$\partial$}g}\,.  (62) 
where H_{1}=\Phi_{1}+\Phi_{1}^{{\rm tid}} contains the effects of the perturbed self–gravity and arbitrary external (slow) perturbation. From this equation for the linear (secular) perturbation F_{1}({\cal R},\tau), we can arrive at the following general conclusion: Any growing or decaying behaviour of F_{1} must necessarily be g–dependent, even in the presence of slow external perturbations of arbitrary form. To see this let us define the g–averaged perturbed DF:
\bar{F}_{1}(I,L,L_{z},h,\tau)\;=\;\oint\frac{{\rm d}g}{2\pi}\,F_{1}(I,L,L_{z},% g,h,\tau)\,,  (63) 
which can be thought of as the DF for stellar orbits, each of which is smeared into an axisymmetric annulus whose inner/outer radii are equal to the peri/apo centre distances. Averaging eqn.(62) over g, we get
\frac{\mbox{$\partial$}\bar{F}_{1}}{\mbox{$\partial$}\tau}\;+\;\nu(I,L)\frac{% \mbox{$\partial$}\bar{F}_{1}}{\mbox{$\partial$}h}\;=\;0\,.  (64) 
The general solution is
\bar{F}_{1}(I,L,L_{z},h,\tau)\;=\;Q(I,L,L_{z},h\nu\tau)\,,  (65) 
where Q(I,L,L_{z},h) is an arbitrary initial g–independent perturbation. According to eqn.(65) such an initial condition evolves through precession of the nodes at the rate \nu(I,L) without either growing or decaying. The nodal precession around the spin axis of the MBH will be negligible in non–relativistic stellar systems. Then the perturbations are quasi–static, and could describe long–lived warps in spherical nuclear clusters. Since the time evolution — in the linear regime — is independent of arbitrary external tidal fields \Phi_{1}^{{\rm tid}}({\cal R},\tau), the perturbed DF retains memory of its origins.
A linear perturbation must be g–dependent for it to either grow or decay. We do not investigate these, and refer the reader to Tremaine (2005); Polyachenko, Polyachenko & Shukhman (2007) which were discussed briefly in the Introduction. Spherical isotropic DFs of the form F_{0}(I) are special because the right hand side of eqn.(62) vanishes. Then the LCBE is:
\frac{\mbox{$\partial$}F_{1}}{\mbox{$\partial$}\tau}\;+\;\Omega\frac{\mbox{$% \partial$}F_{1}}{\mbox{$\partial$}g}\;+\;\nu\frac{\mbox{$\partial$}F_{1}}{% \mbox{$\partial$}h}\;=\;0\,,  (66) 
Therefore an arbitrary infinitesimal perturbation of an isotropic distribution function F_{0}(I) flows collisionlessly along the unperturbed trajectories in {\cal R}–space. The general solution of eqn.(66) is:
F_{1}({\cal R},\tau)\;=\;Q(I,L,L_{z},\,g\Omega\tau,\,h\nu\tau)\,,\qquad% \qquad\mbox{for $\,\tau\geq 0\,$.}  (67) 
where Q({\cal R})=Q(I,L,L_{z},g,h) is an arbitrary initial perturbation that is assumed to be given. This implies: (a) All isotropic DFs F_{0}(I) are secularly stable, because the perturbation evidently is non growing; (b) The evolution of an initial perturbation is independent of H_{1}({\cal R},\tau).
6.2 Stability of Axisymmetric Zero–thickness Flat Discs
We recall from § 4.1 that, for zero–thickness flat discs, the Ring phase space is 3–dim: {\cal R}=(I,L,g), where I\leq L\leq I is the angular momentum, and g is the apsidal longitude. Stationary, axisymmetric DFs are of the form F_{0}(I,L), with the corresponding Hamiltonian of eqn.(55):
H_{0}(I,L)\;=\;\Phi_{0}(I,L)\;\;B_{1}\frac{1}{I^{3}L}\;+\;B_{1.5}\frac{{\rm Sgn% }(L)}{I^{3}L^{2}}\;+\;\Phi_{0}^{\rm ext}(I,L)\,,  (68) 
From eqn.(56) we see that, in the unperturbed system, every Ring precesses rigidly at the apse precession rate \Omega(I,L)=\mbox{$\partial$}H_{0}/\mbox{$\partial$}L\,. A linear perturbation, F_{1}({\cal R},\tau), obeys the LCBE:
\frac{\mbox{$\partial$}F_{1}}{\mbox{$\partial$}\tau}\;+\;\Omega\frac{\mbox{$% \partial$}F_{1}}{\mbox{$\partial$}g}\;=\;\frac{\mbox{$\partial$}F_{0}}{\mbox{$% \partial$}L}\frac{\mbox{$\partial$}H_{1}}{\mbox{$\partial$}g}\,,  (69) 
where
\displaystyle H_{1}(I,L,g,\tau)  \displaystyle\;=  \displaystyle\Phi_{1}(I,L,g,\tau)\;+\;\Phi_{1}^{\rm tid}(I,L,g,\tau)\,,  (70a)  
\displaystyle\Phi_{1}(I,L,g,\tau)  \displaystyle\;=  \displaystyle\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)\,.  (70b) 
The linear problem posed by eqns.(69) and (70b) is still difficult to solve, because of the integral relationship between \Phi_{1} and F_{1} in eqn.(70b).
Below we derive a stability result, by exploiting the symmetry properties of the Ring–Ring interaction potential function \Psi({\cal R},{\cal R}^{\prime}). We write:
\Psi(I,L,g,I^{\prime},L^{\prime},g^{\prime})\;=\;GM_{\bullet}\oint\oint\frac{% {\rm d}w}{2\pi}\,\frac{{\rm d}w^{\prime}}{2\pi}\,\frac{1}{\left\mbox{% \boldmath$r$}\mbox{\boldmath$r$}^{\prime}\right}\,,  (71) 
where \mbox{\boldmath$r$}=(x,y) and \mbox{\boldmath$r$}^{\prime}=(x^{\prime},y^{\prime}) can be written in terms of (I,L,\eta,g) and (I^{\prime},L^{\prime},\eta^{\prime},g^{\prime}), respectively, using eqn.(41); we also recall that {\rm d}w=(1e\cos\eta){\rm d}\eta, and {\rm d}w^{\prime}=(1e^{\prime}\cos\eta^{\prime}){\rm d}\eta^{\prime}. The following symmetries of \Psi({\cal R},{\cal R}^{\prime}) can be verified either by writing the integral in eqn.(71) in explicit form, or more simply by thinking about the gravitational energy between two Gaussian Rings:

\Psi is a real function which is even in both L\, and L^{\prime}\,.

The apsidal longitudes g and g^{\prime} occur in \Psi only in the combination gg^{\prime}.

\Psi is independent of both g and g^{\prime} when one of the two Rings is circular (i.e. when L=\pm I\, or \,L^{\prime}=\pm I^{\prime} or both).

\Psi is invariant under the interchange of the two Rings. This can be achieved by any of the transformations: \{I,\,L\}\leftrightarrow\{I^{\prime},\,L^{\prime}\}\,, or \,g\leftrightarrow g^{\prime}\, or both. In explicit form we have: \,\Psi(I,L,g,I^{\prime},L^{\prime},g^{\prime})=\Psi(I^{\prime},L^{\prime},g,I,% L,g^{\prime}), \;\Psi(I,L,g,I^{\prime},L^{\prime},g^{\prime})=\Psi(I,L,g^{\prime},I^{\prime},% L^{\prime},g), which implies that \Psi(I,L,g,I^{\prime},L^{\prime},g^{\prime})=\Psi(I^{\prime},L^{\prime},g^{% \prime},I,L,g).
We now focus attention on disc eigenmodes and hence set \Phi_{1}^{\rm tid}=0 in eqn.(70a), so that the perturbed Hamiltonian H_{1}=\Phi_{1} is due only to the perturbed self–gravity.^{8}^{8}8However, in the unperturbed Hamiltonian H_{0} of eqn.(68), the external axisymmetric contribution, \Phi_{0}^{\rm ext}(I,L), need not be zero. Since the coefficients in the LCBE of eqn.(69) are independent of the variables, g and \tau, it is natural to develop the perturbations in a Fourier series. We begin with the Fourier–expansion of the Ring–Ring interaction potential function \Psi({\cal R},{\cal R}^{\prime}) in g and g^{\prime}. Using the fact that \Psi depends only on the difference of the apsidal longitudes, we write:
\Psi(I,L,g,I^{\prime},L^{\prime},g^{\prime})\;=\;\sum_{m=\infty}^{\infty}\,C_% {m}(I,L,I^{\prime},L^{\prime})\exp{\left[{\rm i}m(gg^{\prime})\right]}\,.  (72) 
The symmetry properties of \Psi, given above in P1 to P4, translate to the following properties of its Fourier coefficients:

The C_{m} are real functions that are even in L and L^{\prime}.

The C_{m} are even in m: i.e. \,C_{m}=C_{m}\,.

When either L=\pm I\, or \,L^{\prime}=\pm I^{\prime}\, or both, then C_{m}=0 for all m\neq 0.

C_{m}(I,L,I^{\prime},L^{\prime})=C_{m}(I^{\prime},L^{\prime},I,L).
The perturbed DF for an eigenmode with azimuthal wavenumber m and eigenfrequency \omega can be written as:
F_{1}(I,L,g,\tau)\;=\;{\rm Re}\left\{F_{1m}(I,L)\exp{\left[{\rm i}\left(mg% \omega\tau\right)\right]}\right\}\,.  (73) 
where F_{1,m}=F_{1m}^{*}\,. To compute the corresponding perturbation in the self–gravity, we substitute eqns.(73) and (72) in (70b). Then the perturbed Hamiltonian is
\displaystyle H_{1}(I,L,g,\tau)  \displaystyle\;=  \displaystyle\Phi_{1}(I,L,g,\tau)\;=\;{\rm Re}\left\{\Phi_{1m}(I,L)\exp{\left[% {\rm i}\left(mg\omega\tau\right)\right]}\right\}\,,  
\displaystyle\Phi_{1m}(I,L)  \displaystyle\;=  \displaystyle 2\pi\int\,{\rm d}I^{\prime}\,{\rm d}L^{\prime}\,C_{m}(I,L,I^{% \prime},L^{\prime})\,F_{1m}(I^{\prime},L^{\prime})\,.  (74) 
Substituting eqns.(73) and (74) in the LCBE eqn.(69), we obtain the following linear integral eigenvalue problem:
\left[m\,\Omega(I,L)\;\;\omega\,\right]\,F_{1m}\;=\;2\pi m\,\left(\frac{\mbox% {$\partial$}F_{0}}{\mbox{$\partial$}L}\right)\,\int\,{\rm d}I^{\prime}\,{\rm d% }L^{\prime}\,\,C_{m}(I,L,I^{\prime},L^{\prime})\,F_{1m}(I^{\prime},L^{\prime})\,,  (75) 
which must be solved to determine the eigenvalues \omega and eigenfunctions F_{1m}(I,L).
When m=0, eqn.(75) implies that \omega\,=\,0\,. Therefore all DFs of the form F_{0}(I,L) are secularly stable to axisymmetric perturbations. Moreover these perturbations are also time–independent, giving rise to nearby axisymmetric equilibria. It should be noted that this does not imply these discs are stable to fast axisymmetric instabilities: very cold discs violating the Toomre (1964) criterion will be unstable to axisymmetric modes growing over times \sim\mbox{$T_{\rm kep}$}. As we discussed at the end of § 3, MMS analysis is insensitive to these fast modes. Therefore when we consider Keplerian stellar discs, we understand that they must have velocity dispersions that are large enough to make them stable to fast modes.
For non–axisymmetric modes with m\neq 0, the integral eigenvalue problem is, in general, difficult to solve. However, we can prove a stability result for a class of stationary DFs. Let F_{0}(I,L) be a strictly monotonic function of L at fixed I. Then we can write
\frac{\mbox{$\partial$}F_{0}}{\mbox{$\partial$}L}\;=\;\sigma\left\frac{\mbox{% $\partial$}F_{0}}{\mbox{$\partial$}L}\right\;\neq\;0\,,  (76) 
where \sigma=+1 when F_{0} is an increasing function of L, and \sigma=1 when F_{0} is a decreasing function of L. Let us define a new eigenfunction, X_{m}(I,L), and a new kernel A_{m}(I,L,I^{\prime},L^{\prime}) by
\displaystyle X_{m}(I,L)  \displaystyle\;=  \displaystyle\left\frac{\mbox{$\partial$}F_{0}}{\mbox{$\partial$}L}\right^{% 1/2}F_{1m}(I,L)\,,  (77a)  
\displaystyle A_{m}(I,L,I^{\prime},L^{\prime})  \displaystyle\;=  \displaystyle 2\pi\left\frac{\mbox{$\partial$}F_{0}}{\mbox{$\partial$}L}% \right^{1/2}C_{m}(I,L,I^{\prime},L^{\prime})\left\frac{\mbox{$\partial$}F^{% \prime}_{0}}{\mbox{$\partial$}L^{\prime}}\right^{1/2}\,,  (77b) 
where F^{\prime}_{0}=F_{0}(I^{\prime},L^{\prime}). Then, for m\neq 0, the eigenvalue eqn.(75) becomes
\left[\,\Omega(I,L)\;\;\frac{\omega}{m}\,\right]X_{m}(I,L)\;=\;\sigma\int\,{% \rm d}I^{\prime}\,{\rm d}L^{\prime}\,\,A_{m}(I,L,I^{\prime},L^{\prime})\,X_{m}% (I^{\prime},L^{\prime})\,.  (78) 
From F2 and F4 we see that the kernel, A_{m}(I,L,I^{\prime},L^{\prime}) is real, even in m, and symmetric under the interchange \{I,\,L\}\leftrightarrow\{I^{\prime},\,L^{\prime}\}\,. Hence A_{m} is Hermitian, and all the eigenvalues \omega are real for all m\,. Therefore we have proved:
\bullet\, Stationary, axisymmetric discs with DFs F_{0}(I,L) are neutrally stable to all secular perturbations when \,\mbox{$\partial$}F_{0}/\mbox{$\partial$}L\, is of the same sign (either positive or negative) everywhere in its domain of support, \,I\leq L\leq I\, and I_{\rm min}\leq I\leq I_{\rm max}\,.
These secularly stable DFs can have both prograde and retrograde populations of stars because I\leq L\leq I\,. Hence the stability result applies to counter–rotating discs of stars with DFs that are strictly monotonic (either increasing or decreasing) functions of L, at fixed I. Depending on the sign of \mbox{$\partial$}F_{0}/\mbox{$\partial$}L, either prograde or retrograde stars dominate at every value of I\,. These secularly stable stellar discs have net rotation and include physically interesting cases, such as a secular analogue of the well–known Schwarzschild DF. We may also restate the above stability result as:
\bullet\; A necessary condition for F_{0}(I,L) to be secularly unstable is that \mbox{$\partial$}F_{0}/\mbox{$\partial$}L must vanish somewhere in its domain of support, \,I\leq L\leq I\, and I_{\rm min}\leq I\leq I_{\rm max}\,.
Jalali & Tremaine (2012) derived the dispersion relationship for slow modes in a collisionless disc with a Schwarzschild DF — using the epicyclic approximation and the WKB limit — and showed that modes of all m were stable. The stability result presented above is more general, and not restricted to a particular DF, or the epicyclic approximation or the WKB limit. Secular instabilities have been discussed in Tremaine (2005); Polyachenko, Polyachenko & Shukhman (2007). Tremaine (2005) considered non–rotating stellar discs described by DFs F_{0}(I,L) that are even functions of L\,, and empty loss cones (F_{0}/L nonsingular when L\to 0\,). Using a variational principle (Goodman, 1988) he argued that many of them are likely to be unstable to m=1 modes. Polyachenko, Polyachenko & Shukhman (2007) considered mono–energetic counter–rotating discs dominated by nearly radial orbits, with DFs equivalent to F_{0}\propto\delta(II_{0})A(L)\,, and found that they are prone to loss cone instabilities of all m\,.
7 Discussion
The main results of this paper are:

The full DF in 6–dim phase space, of a secularly evolving, collisionless Keplerian stellar system, is given by eqn.(28). This the sum of an O(1) secular DF in a reduced 5–dim phase space, and small fluctuations that remain of O(\varepsilon) over secular times \mbox{$T_{\rm sec}$}=\varepsilon^{1}\mbox{$T_{\rm kep}$}\,.

The slow time evolution of the secular (or Ring) DF F({\cal R},\tau) in 5–dim {\cal R}–space is governed by the Ring CBE eqn.(38), which includes the combined secular effects of the gravity of the MBH upto 1.5 PN order, the self–gravity of all the stars, and slowly varying external sources of arbitrary form.
The lower dimensional phase space results from the secular conservation of the semi–major axis of every stellar orbit. From this fully nonlinear, self–consistent formulation of secular dynamics follows a Secular Jeans theorem enabling the construction of Secular Equilibria, and problems concerning the Linear Response and Stability of these equilibria. Secular equilibria are easier to construct than full stellar dynamical equilibria. For instance, Jeans theorem implies that the exact DF of a stationary, spherical, non–rotating system is of the form f(E,L). In order to determine f as a function of the phase space coordinates, \{\mbox{\boldmath$r$},\mbox{\boldmath$u$}\}, it is necessary to solve a self–consistent problem. However, the secular DF for the same system, F(I,L), which differs from f only by O(\varepsilon), is known explicitly as a function of \{\mbox{\boldmath$r$},\mbox{\boldmath$u$}\}.
Linear secular stability of disc–like and spherical configurations has been studied earlier (Tremaine, 2005; Polyachenko, Polyachenko & Shukhman, 2007). Both begin with equilibria and linear perturbations of the full (non–secular) problem. Having formulated an eigenvalue equation for normal modes, they then take the slow mode — i.e. \omega\sim O(\varepsilon) — limit. In contrast we have formulated a fully nonlinear secular theory, described by Gaussian Ring dynamics. From this theory followed secular equilibria and development of secular linear perturbation theory. Broad conclusions were reached about evolutions of small perturbations in spherical non–rotating systems. Tremaine (2005) has noted that flat systems are less simple than spherical systems due to a difficulty with defining inner–products, and the more complicated relation between mass distribution and potential. We considered perturbations of DFs describing axisymmetric, zero–thickness flat discs, and derived a secular integral eigenvalue problem. Even though the potential is related to the surface density in a more complicated manner, its symmetry properties were sufficient to let us prove a stability result for a class of physically interesting rotating DFs.
The fully nonlinear secular description of the Ring CBE is needed to go beyond linear theory. Touma & Sridhar (2012) used a fully nonlinear secular CBE — equivalent to the Ring CBE — to discuss the counter–rotating instability, but they did not derive the secular equation from the full CBE, as we have done in this paper. Touma, Tremaine & Kazandjian (2009) performed numerical simulations of secular dynamics, wherein each star is replaced by a Gaussian Ring that responds to the torques of all the other stars. They found that the counter–rotating instability in unstable discs saturates over times \sim\mbox{$T_{\rm sec}$}, and the discs settle into uniformly precessing equilibria. The in–plane instability unfolds into a triaxial one, which couples radial (or eccentricity) and vertical (or inclination) degrees–of–freedom, resulting in the dispersal of the lighter of the two counter–rotating components into a triaxial halo in which the other, more massive, precessing, lopsided disc is embedded. This process is the secular analogue of “violent relaxation” (Lynden–Bell, 1967) and forms a central problem in collisionless secular dynamics. The results obtained in this paper are used directly in the papers to follow. In Paper II we formulate a theory of the resonant relaxation of general Keplerian stellar systems, in terms of the quasi–static evolution (due to gravitational collisions between Gaussian Rings) of the Ring DF through a sequence of collisionless equilibria, each of which must be dynamically stable. In Paper III we apply this theory to study the physical kinetics of the resonant relaxation of axisymmetric discs.
Acknowledgments
This work has benefited from a series of visits by S. Sridhar to Ras Beirut, made possible by the generous support of the Dean’s Office and the hospitality of the Department of Physics, at the American University of Beirut. We thank Karamveer Kaur for a careful reading of an earlier draft.
References
 Bacon et al. (2001) Bacon R., Emsellem E., Combes F., Copin Y., Monnet G., Martin P., 2001, A&A, 371, 409
 Bahcall & Wolf (1976) Bahcall J. N., Wolf R. A., 1976, ApJ, 209, 214
 Bender & Orszag (1978) Bender C. M., Orszag S. A., 1978, Advanced Mathematical Methods for Scientists and Engineers, McGrawHill, New York
 Bender et al. (2005) Bender, R., Kormendy, J., Bower, G., et al. 2005, ApJ, 631, 280
 Binney & Tremaine (2008) Binney, J., & Tremaine, S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton, NJ
 Brown & Magorrian (2013) Brown C. K., Magorrian, J., 2013, MNRAS, 431, 80
 Cohn & Kulsrud (1978) Cohn H., Kulsrud R. M., 1978, ApJ, 226, 1087
 Genzel, Eisenhauer & Gillessen (2010) Genzel R., Eisenhauer F., Gillessen, S., 2010, Reviews of Modern Physics, 82, 3121
 Gerhard & Binney (1985) Gerhard O. E., Binney, J., 1985, MNRAS, 216, 467
 Gilbert (1968) Gilbert, I. H., 1968, The Astrophysical Journal, 152, 1043.
 Goodman & Binney (1984) Goodman J., Binney, J., 1984, MNRAS, 207, 511
 Goodman (1988) Goodman J., 1988, ApJ, 329, 612
 Gulati, Saini & Sridhar (2012) Gulati M., Saini T. D., Sridhar S., 2012, MNRAS, 424, 348
 Heckman & Best (2014) Heckman T. M., Best P. N., 2014, ARA&A, 52, 589
 Jacobs & Sellwood (2001) Jacobs V., Sellwood J. A., 2001, ApJ, 555, L25
 Jalali & Tremaine (2012) Jalali M. A., & Tremaine S., 2012, MNRAS, 421, 2368
 Kazandjian & Touma (2013) Kazandjian M. V., Touma J. R., 2013, MNRAS, 430, 2732
 Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
 Kormendy & Richstone (1995) Kormendy J., Richstone, D., 1995, ARA&A, 33, 581
 Krolik (1999) Krolik J. H., 1999, Active Galactic Nuclei: from the central black hole to the galactic environment, Princeton Univ. Press, Princeton, NJ
 Lee & Goodman (1999) Lee E., Goodman J., 1999, MNRAS, 308, 984
 Lynden–Bell (1967) Lynden–Bell D., 1967, MNRAS, 136, 101
 Merritt & Valluri (1999) Merritt D., Valluri, M., 1999, AJ, 118, 1177
 Merritt & Vasiliev (2011) Merritt D., Vasiliev, E., 2011, ApJ, 726, 61
 Merritt (2013) Merritt D., 2013, Dynamics and Evolution of Galactic Nuclei, Princeton Univ. Press, Princeton, NJ
 Murray & Dermott (1999) Murray C. D., Dermott, S. F., 1999, Solar System Dynamics, Cambridge University Press, Cambridge
 Peiris & Tremaine (2003) Peiris H. V., Tremaine, S., 2003, ApJ, 599, 237
 Plummer (1960) Plummer H. C., 1960, An Introductory Treatise on Dynamical Astronomy, Dover, New York
 Polyachenko, Polyachenko & Shukhman (2007) Polyachenko E. V., Polyachenko V. L., Shukhman I. G., 2007, MNRAS, 379, 573
 Poon & Merritt (2001) Poon M. Y., Merritt, D., 2001, ApJ, 549, 192
 Rauch & Tremaine (1996) Rauch K. P., Tremaine S., 1996, New Astronomy, 1, 149
 Rees (1984) Rees M. J. 1984, ARA&A, 22, 471
 Salow & Statler (2001) Salow R. M., & Statler T. S., 2001, ApJ, 551, L49
 Salow & Statler (2004) Salow R. M., Statler, T. S., 2004, ApJ, 611, 245
 Sambhus & Sridhar (2000) Sambhus N., Sridhar, S., 2000, ApJ, 542, 143
 Sambhus & Sridhar (2002) Sambhus N., Sridhar, S., 2002, A&A, 388, 766
 Sridhar & Saini (2010) Sridhar S., Saini T. D., 2010, MNRAS, 404, 527
 Sridhar, Syer & Touma (1999) Sridhar S., Syer D., Touma J., 1999, in ASP Conf. Ser. 160, Astrophysical Discs – an EC Summer School, ed. J. A. Sellwood & J. Goodman, San Francisco, 307
 Sridhar & Touma (1997a) Sridhar S., Touma, J., 1997, MNRAS, 287, L1
 Sridhar & Touma (1997b) Sridhar S., Touma, J., 1997, MNRAS, 292, 657
 Sridhar & Touma (1999) Sridhar S. & Touma J., 1999, MNRAS, 303, 483
 Sridhar & Touma (2015; Paper II) Sridhar S., Touma J. R., 2015, arXiv:1509.02401 (Paper II)
 Sridhar & Touma (2016; Paper III) Sridhar S., Touma J. R., 2016, arXiv:1602.05763 (Paper III)
 Statler (1999) Statler T. S., 1999, ApJ, 524, L87
 Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
 Touma (2002) Touma J. R., 2002, MNRAS, 333, 583
 Touma & Sridhar (2012) Touma J. R., Sridhar, S., 2012, MNRAS, 423, 2083
 Touma, Tremaine & Kazandjian (2009) Touma J. R., Tremaine S., Kazandjian M. V., 2009, MNRAS, 394, 1085
 Touma & Tremaine (2014) Touma J., Tremaine S., 2014, J. Phys. A, 47, 292001
 Tremaine (1995) Tremaine S., 1995, AJ, 110, 628
 Tremaine (2001) Tremaine S., 2001, AJ, 121, 1776
 Tremaine (2005) Tremaine S., 2005, ApJ, 625, 143
 Yelda et al. (2014) Yelda, S., Ghez, A. M., Lu, J. R., et al. 2014, ApJ, 783, 131
 Young (1980) Young P., 1980, ApJ, 242, 1232
Appendix A Secular Dynamics of Stellar Orbits
The full DF, f({\cal D},t) of eqn.(28), describes the collisionless behaviour of N\gg 1 stars around a MBH. The orbit of a test star that belongs to this DF is governed by the Hamiltonian H_{\rm org} of eqn.(6). H_{\rm org} is the sum of three terms: of these the one containing \mbox{\boldmath$A$}_{\bullet} vanishes to O(\varepsilon) accuracy (see eqn.20); the Kepler energy is E_{\rm k}(I)=1/2\left(GM_{\bullet}/I\right)^{2}; and the self–gravitational potential is
\displaystyle\varepsilon\,GM_{\bullet}\int\frac{f({\cal D}^{\prime},t)}{\,% \mbox{\boldmath$r$}\mbox{\boldmath$r$}^{\prime}\,}\,{\rm d}{\cal D}^{\prime}  \displaystyle\;=  \displaystyle\varepsilon\,GM_{\bullet}\int F({\cal R}^{\prime},\varepsilon t)% \,{\rm d}{\cal R}^{\prime}\oint\frac{{\rm d}w^{\prime}}{2\pi}\frac{1}{\,\mbox% {\boldmath$r$}\mbox{\boldmath$r$}^{\prime}\,}\;+\;O(\varepsilon^{2})  (79)  
\displaystyle\;=  \displaystyle\varepsilon\Phi({\cal R},\varepsilon t)\;+\;\varepsilon\sum_{n% \neq 0}\varphi_{n}({\cal R},\varepsilon t)\exp{[{\rm i}nw]}\;+\;O(\varepsilon^% {2})\,, 
where we have used eqns.(16) and (17). Hence, to O(\varepsilon) accuracy,
H_{\rm org}({\cal D},t)\;=\;E_{\rm k}(I)\;+\;\varepsilon\Phi({\cal R},% \varepsilon t)\;+\;\varepsilon\sum_{n\neq 0}\varphi_{n}({\cal R},\varepsilon t% )\exp{[{\rm i}nw]}\,,  (80) 
depends on all six Delaunay variables, but has only slow variations with time. Therefore it is natural to want to use \tau=\varepsilon t as the time variable, instead of t. The equations of motion are invariant under the rescaling, \{t\,,\;H_{\rm org}\}\to\{\tau=\varepsilon t\,,\;H^{\prime}_{\rm org}=% \varepsilon^{1}H_{\rm org}\}. We have:
H^{\prime}_{\rm org}({\cal D},\tau)\;=\;\frac{E_{\rm k}(I)}{\varepsilon}\;+\;% \Phi({\cal R},\tau)\;+\;\sum_{n\neq 0}\varphi_{n}({\cal R},\tau)\exp{[{\rm i}% nw]}\,.  (81) 
The next step makes use of a near–identity canonical transformation from old to new Delaunay variables, {\cal D}\to{\cal D}^{\prime}\equiv\{I^{\prime},L^{\prime},L^{\prime}_{z};w^{% \prime},g^{\prime},h^{\prime}\}, through the mixed–variable generating function,
{\cal S}_{0}(I^{\prime},L^{\prime},L^{\prime}_{z};w,g,h;\tau)\;=\;I^{\prime}w% \,+\,L^{\prime}g\,+\,L^{\prime}_{z}h\;+\;\varepsilon{\cal S}(I^{\prime},L^{% \prime},L^{\prime}_{z};w,g,h;\tau)\,.  (82) 
The old actions (I,L,L_{z}) and new angles (w^{\prime},g^{\prime},h^{\prime}) are given by:
\displaystyle I  \displaystyle\;=  \displaystyle I^{\prime}\;+\;\varepsilon\frac{\mbox{$\partial$}{\cal S}}{\mbox% {$\partial$}w}\,,\qquad\qquad L\;=\;L^{\prime}\;+\;\varepsilon\frac{\mbox{$% \partial$}{\cal S}}{\mbox{$\partial$}g}\,,\qquad\qquad L_{z}\;=\;L^{\prime}_{z% }\;+\;\varepsilon\frac{\mbox{$\partial$}{\cal S}}{\mbox{$\partial$}h}\,,  
\displaystyle w^{\prime}  \displaystyle\;=  \displaystyle w\;+\;\varepsilon\frac{\mbox{$\partial$}{\cal S}}{\mbox{$% \partial$}I^{\prime}}\,,\qquad\qquad g^{\prime}\;=\;g\;+\;\varepsilon\frac{% \mbox{$\partial$}{\cal S}}{\mbox{$\partial$}L^{\prime}}\,,\qquad\qquad h^{% \prime}\;=\;h\;+\;\varepsilon\frac{\mbox{$\partial$}{\cal S}}{\mbox{$\partial$% }L^{\prime}_{z}}\,.  (83) 
Hence {\cal D} and {\cal D}^{\prime}, the old and new variables, differ only by O(\varepsilon). The new Hamiltonian governing the dynamics of {\cal D}^{\prime} is:
\displaystyle H^{{}^{\prime\prime}}_{\rm org}({\cal D}^{\prime},\tau)  \displaystyle\;=  \displaystyle H^{\prime}_{0}({\cal D},\tau)\;+\;\varepsilon\frac{\mbox{$% \partial$}{\cal S}}{\mbox{$\partial$}\tau}  (84)  
\displaystyle\;=  \displaystyle\frac{E_{\rm k}(I)}{\varepsilon}\;+\;\Phi({\cal R},\tau)\;+\;\sum% _{n\neq 0}\varphi_{n}({\cal R},\tau)\exp{[{\rm i}nw]}\;+\;\varepsilon\frac{% \mbox{$\partial$}{\cal S}}{\mbox{$\partial$}\tau}\,. 
We need to express the right hand side in terms of the new variables {\cal D}^{\prime}, and keep only the O(1) terms. The first term can be expanded as:
\frac{E_{\rm k}(I)}{\varepsilon}\;=\;\frac{1}{\varepsilon}E_{\rm k}\!\left(I^{% \prime}+\varepsilon\,\frac{\mbox{$\partial$}{\cal S}}{\mbox{$\partial$}w}% \right)\;=\;\frac{E_{\rm k}(I^{\prime})}{\varepsilon}\;+\;\Omega_{\rm k}(I^{% \prime})\frac{\mbox{$\partial$}{\cal S^{\prime}}}{\mbox{$\partial$}w^{\prime}}% \;+\;O(\varepsilon)\,, 
where {\cal S}^{\prime}({\cal D}^{\prime})={\cal S}(I^{\prime},L^{\prime},L^{\prime}% _{z};w^{\prime},g^{\prime},h^{\prime};\tau) is a function only of the new variables, obtained by replacing the old angles (w,g,h) in {\cal S} by the new ones (w^{\prime},g^{\prime},h^{\prime}); this is legitimate, because, by eqns.(83), the old and new variables differ only by O(\varepsilon). For the same reason, in the two O(1) terms, we can replace {\cal D} by {\cal D}^{\prime}. Then
H^{{}^{\prime\prime}}_{\rm org}({\cal D}^{\prime},\tau)\;=\;\frac{E_{\rm k}(I^% {\prime})}{\varepsilon}\;+\;\Omega_{\rm k}(I^{\prime})\frac{\mbox{$\partial$}{% \cal S}^{\prime}}{\mbox{$\partial$}w^{\prime}}\;+\;\Phi({\cal R}^{\prime},\tau% )\;+\;\sum_{n\neq 0}\varphi_{n}({\cal R}^{\prime},\tau)\exp{[{\rm i}nw^{\prime% }]}\;+\;O(\varepsilon)\,.  (85) 
We now choose {\cal S}^{\prime} such that all the w^{\prime}–dependent terms on the right hand side cancel. This amounts to choosing the mixed variable generating function,
{\cal S}(I^{\prime},L^{\prime},L^{\prime}_{z};w,g,h;\tau)\;=\;\frac{1}{\Omega% _{\rm k}(I^{\prime})}\,\sum_{n\neq 0}\frac{1}{{\rm i}n}\varphi_{n}(I^{\prime},% L^{\prime},L^{\prime}_{z},g,h,\tau)\exp{[{\rm i}nw]}\,.  (86) 
{\cal S} is well–defined, since \Omega_{\rm k}(I^{\prime}) and n are non–zero. Note that the degeneracy of the Kepler problem has prevented the occurrence of small denominators. Then the Hamiltonian for the dynamics of {\cal D}^{\prime} is:
H^{{}^{\prime\prime}}_{\rm org}({\cal D}^{\prime},\tau)\;=\;\frac{E_{\rm k}(I^% {\prime})}{\varepsilon}\;+\;\Phi({\cal R}^{\prime},\tau)\;+\;O(\varepsilon)\,.  (87) 
with the corresponding equations of motion:
\displaystyle\frac{{\rm d}I^{\prime}}{{\rm d}\tau}  \displaystyle\;=  \displaystyle O(\varepsilon)\,,\qquad\frac{{\rm d}w^{\prime}}{{\rm d}\tau}\;=% \;\frac{\Omega_{\rm k}(I^{\prime})}{\varepsilon}\;+\;\frac{\mbox{$\partial$}% \Phi}{\mbox{$\partial$}I^{\prime}}\;+\;O(\varepsilon)\,;  
\displaystyle\frac{{\rm d}L^{\prime}}{{\rm d}\tau}  \displaystyle\;=  \displaystyle\,\frac{\mbox{$\partial$}\Phi}{\mbox{$\partial$}g^{\prime}}\;+\;% O(\varepsilon)\,,\qquad\frac{{\rm d}g^{\prime}}{{\rm d}\tau}\;=\;\frac{\mbox{$% \partial$}\Phi}{\mbox{$\partial$}L^{\prime}}\;+\;O(\varepsilon)\,;  
\displaystyle\frac{{\rm d}L^{\prime}_{z}}{{\rm d}\tau}  \displaystyle\;=  \displaystyle\,\frac{\mbox{$\partial$}\Phi}{\mbox{$\partial$}h^{\prime}}\;+\;% O(\varepsilon)\,,\qquad\frac{{\rm d}h^{\prime}}{{\rm d}\tau}\;=\;\frac{\mbox{$% \partial$}\Phi}{\mbox{$\partial$}L^{\prime}_{z}}\;+\;O(\varepsilon)\,.  (88) 
The dominant behaviour of the orbital phase w^{\prime} is rapid circulation at the Kepler orbital rate, (\Omega_{\rm k}/\varepsilon). Superposed on this are two kinds of modulations; a slow O(1) variation and a fast O(\varepsilon) oscillation. The other five variables, {\cal R}^{\prime}=(I^{\prime},L^{\prime},L^{\prime}_{z},g^{\prime},h^{\prime}), behave like Gaussian Rings with additional fast O(\varepsilon) oscillations: to O(1) accuracy, I^{\prime} is constant, while the other four variables, (L^{\prime},L^{\prime}_{z};g^{\prime},h^{\prime}), vary over secular time scales with \Phi({\cal R}^{\prime},\tau) acting as the Hamiltonian for their dynamics. This is just the dynamics of {\cal R}^{\prime} in the orbit–averaged self–gravitational potential. Since, by eqn.(83), {\cal R} and {\cal R}^{\prime} differ only by O(\varepsilon), the old {\cal R} variables also have O(1) slow dynamics governed by the Hamiltonian \Phi({\cal R},\tau) plus O(\varepsilon) fast oscillations. Therefore a stellar orbit in a Keplerian stellar system can be thought of as a secularly evolving, “noisy” Gaussian Ring.\hfill\square