The centers of most galaxies contain massive black holes surrounded by dense star clusters. The structure of these clusters determines the rate and properties of observable transient events, such as flares from tidally disrupted stars and gravitational-wave signals from stars spiraling into the black hole. Most estimates of these rates enforce spherical symmetry on the cluster. Here we show that, in the course of generic evolutionary processes, a star cluster surrounding a black hole can undergo a robust phase transition from a spherical thermal equilibrium to a lopsided equilibrium, in which most stars are on high-eccentricity orbits with aligned orientations. The rate of transient events is expected to be much higher in the ordered phase. Better models of cluster formation and evolution are needed to determine whether clusters should be found in the ordered or disordered phase.

A relatively massive and moderately eccentric disk of trans-Neptunian objects (TNOs) can effectively counteract apse precession induced by the outer planets, and in the process shepherd highly eccentric members of its population into nearly-stationary configurations which are anti-aligned with the disk itself. We were sufficiently intrigued by this remarkable feature to embark on an extensive exploration of the full spatial dynamics sustained by the combined action of giant planets and a massive trans-Neptunian debris disk. In the process, we identified ranges of disk mass, eccentricity and precession rate which allow apse-clustered populations that faithfully reproduce key orbital properties of the much discussed TNO population. The shepherding disk hypothesis is to be sure complementary to any potential ninth member of the Solar System pantheon, and could obviate the need for it altogether. We discuss its essential ingredients in the context of Solar System formation and evolution, and argue for their naturalness in view of the growing body of observational and theoretical knowledge about self-gravitating disks around massive bodies, extra-solar debris disks included.

We present idealized models of a razor-thin, axisymmetric, Keplerian stellar disc around a massive black hole, and study non-axisymmetric secular instabilities in the absence of either counter-rotation or loss cones. These discs are prograde mono-energetic waterbags, whose phase space distribution functions are constant for orbits within a range of eccentricities (e) and zero outside this range. The linear normal modes of waterbags are composed of sinusoidal disturbances of the edges of distribution function in phase space. Waterbags which include circular orbits (polarcaps) have one stable linear normal mode for each azimuthal wavenumber m. The m = 1 mode always has positive pattern speed and, for polarcaps consisting of orbits with e < 0.9428, only the m = 1 mode has positive pattern speed. Waterbags excluding circular orbits (bands) have two linear normal modes for each m, which can be stable or unstable. We derive analytical expressions for the instability condition, pattern speeds, growth rates and normal mode structure. Narrow bands are unstable to modes with a wide range in m. Numerical simulations confirm linear theory and follow the non-linear evolution of instabilities. Long-time integration suggests that instabilities of different m grow, interact non-linearly and relax collisionlessly to a coarse-grained equilibrium with a wide range of eccentricities.

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.

We present a first-principles theory of Resonant Relaxation (RR) of a low mass stellar system orbiting a more massive black hole (MBH). We first extend the kinetic theory of Gilbert (1968) to include the Keplerian field of a black hole of mass $M_\bullet$. Specializing to a Keplerian stellar system of mass $M \ll M_\bullet$, we use the orbit-averaging method of Sridhar & Touma (2015; Paper I) to derive a kinetic equation for RR. This describes the collisional evolution of a system of $N \gg 1$ Gaussian Rings in a reduced 5-dim space, under the combined actions of self-gravity, 1 PN and 1.5 PN relativistic effects of the MBH and an arbitrary external potential. In general geometries RR is driven by both apsidal and nodal resonances, so the distinction between scalar-RR and vector-RR disappears. The system passes through a sequence of quasi-steady secular collisionless equilibria, driven by irreversible 2-Ring correlations that accrue through gravitational interactions, both direct and collective. This correlation function is related to a `wake function', which is the linear response of the system to the perturbation of a chosen Ring. The wake function is easier to appreciate, and satisfies a simpler equation, than the correlation function. We discuss general implications for the interplay of secular dynamics and non-equilibrium statistical mechanics in the evolution of Keplerian stellar systems toward secular thermodynamic equilibria, and set the stage for applications to the RR of axisymmetric discs in Paper III.

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

Most exoplanetary systems in binary stars are of S--type, and consist of one or more planets orbiting a primary star with a wide binary stellar companion. Gravitational forcing of a single planet by a sufficiently inclined binary orbit can induce large amplitude oscillations of the planet's eccentricity and inclination through the Kozai-Lidov (KL) instability. KL cycling was invoked to explain: the large eccentricities of planetary orbits; the family of close--in hot Jupiters; and the retrograde planetary orbits in eccentric binary systems. However, several kinds of perturbations can quench the KL instability, by inducing fast periapse precessions which stabilize circular orbits of all inclinations: these could be a Jupiter--mass planet, a massive remnant disc or general relativistic precession. Indeed, mutual gravitational perturbations in multiplanet S--type systems can be strong enough to lend a certain dynamical rigidity to their orbital planes. Here we present a new and faster process that is driven by this very agent inhibiting KL cycling. Planetary perturbations enable secular oscillations of planetary eccentricities and inclinations, also called Laplace--Lagrange (LL) eigenmodes. Interactions with a remnant disc of planetesimals can make planets migrate, causing a drift of LL mode periods which can bring one or more LL modes into resonance with binary orbital motion. The results can be dramatic, ranging from excitation of large eccentricities and mutual inclinations to total disruption. Not requiring special physical or initial conditions, binary resonant driving is generic and could have profoundly altered the architecture of many S--type multiplanet systems. It can also weaken the multiplanet occurrence rate in wide binaries, and affect planet formation in close binaries.

We describe the dynamics and thermodynamics of collisionless particle disks orbiting a massive central body, in the case where the disk mass is small compared to the central mass, the self-gravity of the disk dominates the non-Keplerian force, and the spread in semi-major axes is small. We show that with plausible approximations such disks have logarithmic two-body interactions and a compact phase space, and therefore exhibit thermodynamics that are simpler than most other gravitating systems, which require a confining box and artificial softening of the potential at small scales to be thermodynamically well-behaved. We solve for the microcanonical axisymmetric thermal equilibria and demonstrate the existence of a symmetry-breaking bifurcation into lopsided equilibria. We discuss the relation between thermal and dynamical instability in these systems and draw connections to astrophysical settings, as well as to the wider subject of the statistical mechanics of particles with logarithmic long-range interactions, such as point vortices in two-dimensional fluids.

We revisit in one dimension the waterbag method to solve numerically Vlasov-Poisson equations. In this approach, the phase-space distribution function $f(x,v)$ is initially sampled by an ensemble of patches, the waterbags, where $f$ is assumed to be constant. As a consequence of Liouville theorem it is only needed to follow the evolution of the border of these waterbags, which can be done by employing an orientated, self-adaptive polygon tracing isocontours of $f$. This method, which is entropy conserving in essence, is very accurate and can trace very well non linear instabilities as illustrated by specific examples. As an application of the method, we generate an ensemble of single waterbag simulations with decreasing thickness, to perform a convergence study to the cold case. Our measurements show that the system relaxes to a steady state where the gravitational potential profile is a power-law of slowly varying index $\beta$, with $\beta$ close to $3/2$ as found in the literature. However, detailed analysis of the properties of the gravitational potential shows that at the center, $\beta > 1.54$. Moreover, our measurements are consistent with the value $\beta=8/5=1.6$ that can be analytically derived by assuming that the average of the phase-space density per energy level obtained at crossing times is conserved during the mixing phase. These results are incompatible with the logarithmic slope of the projected density profile $\beta-2 \simeq -0.47$ obtained recently by Schulz et al. (2013) using a $N$-body technique. This sheds again strong doubts on the capability of $N$-body techniques to converge to the correct steady state expected in the continuous limit.

It is strongly believed that Andromeda's double nucleus signals a disk of stars revolving around its central super-massive black hole on eccentric Keplerian orbits with nearly aligned apsides. A self-consistent stellar dynamical origin for such apparently long-lived alignment has so far been lacking, with indications that cluster self-gravity is capable of sustaining such lopsided configurations if and when stimulated by external perturbations. Here, we present results of N-body simulations which show unstable counter-rotating stellar clusters around super-massive black holes saturating into uniformly precessing lopsided nuclei. The double nucleus in our featured experiment decomposes naturally into a thick eccentric disk of apo-apse aligned stars which is embedded in a lighter triaxial cluster. The eccentric disk reproduces key features of Keplerian disk models of Andromeda's double nucleus; the triaxial cluster has a distinctive kinematic signature which is evident in HST observations of Andromeda's double nucleus, and has been difficult to reproduce with Keplerian disks alone. Our simulations demonstrate how the combination of eccentric disk and triaxial cluster arises naturally when a star cluster accreted over a pre-existing and counter-rotating disk of stars, drives disk and cluster into a mutually destabilizing dance. Such accretion events are inherent to standard galaxy formation scenarios. They are here shown to double stellar black hole nuclei as they feed them.

We formulate the collisionless Boltzmann equation (CBE) for dense star clusters that lie within the radius of influence of a massive black hole in galactic nuclei. Our approach to these nearly Keplerian systems follows that of Sridhar and Touma (1999): Delaunay canonical variables are used to describe stellar orbits and we average over the fast Keplerian orbital phases. The stellar distribution function (DF) evolves on the longer time scale of precessional motions, whose dynamics is governed by a Hamiltonian, given by the orbit-averaged self-gravitational potential of the cluster. We specialize to razor-thin, planar discs and consider two counter-rotating ("$\pm$") populations of stars. To describe discs of small eccentricities, we expand the $\pm$ Hamiltonian to fourth order in the eccentricities, with coefficients that depend self-consistently on the $\pm$ DFs. We construct approximate $\pm$ dynamical invariants and use Jeans' theorem to construct time-dependent $\pm$ DFs, which are completely described by their centroid coordinates and shape matrices. When the centroid eccentricities are larger than the dispersion in eccentricities, the $\pm$ centroids obey a set of 4 autonomous ordinary differential equations. We show that these can be cast as a two-degree of freedom Hamiltonian system which is nonlinear, yet integrable. We study the linear instability of initially circular discs and derive a criterion for the counter-rotating instability. We then explore the rich nonlinear dynamics of counter-rotating discs, with focus on the variety of steadily precessing eccentric configurations that are allowed. The stability and properties of these configurations are studied as functions of parameters such as the disc mass ratios and angular momentum.

We show that the algorithm proposed by Gauss to compute the secular evolution of gravitationally interacting Keplerian rings extends naturally to softened gravitational interactions. The resulting tool is ideal for the study of the secular dynamical evolution of nearly Keplerian systems such as stellar clusters surrounding black holes in galactic nuclei, cometary clouds, or planetesimal discs. We illustrate its accuracy, efficiency and versatility on a variety of configurations. In particular, we examine a secularly unstable unstable system of counter-rotating disks, and follow the unfolding and saturation of the instability into a global, uniformly precessing, lopsided (m=1) mode.

The orbital dynamics of most planetary satellites is governed by the quadrupole moment from the equatorial bulge of the host planet and the tidal field from the Sun. On the Laplace surface, the long-term orbital evolution driven by the combined effects of these forces is zero, so that orbits have a fixed orientation and shape. The "classical" Laplace surface is defined for circular orbits, and coincides with the planet's equator at small planetocentric distances and with its orbital plane at large distances. A dissipative circumplanetary disk should settle to this surface, and hence satellites formed from such a disk are likely to orbit in or near the classical Laplace surface. This paper studies the properties of Laplace surfaces. Our principal results are: (i) if the planetary obliquity exceeds 68.875 deg there is a range of semimajor axes in which the classical Laplace surface is unstable; (ii) at some obliquities and planetocentric distances there is a distinct Laplace surface consisting of nested eccentric orbits, which bifurcates from the classical Laplace surface at the point where instability sets in; (iii) there is also a "polar" Laplace surface perpendicular to the line of nodes of the planetary equator on the planetary orbit; (iv) for circular orbits, the polar Laplace surface is stable at small planetocentric distances and unstable at large distances; (v) at the onset of instability this polar Laplace surface bifurcates into two polar Laplace surfaces composed of nested eccentric orbits.

Caustics are formally singular structures, with infinite density, that form in collisionless media. The non-negligible velocity dispersion of dark matter particles renders their density finite. We evaluate the maximum density of the caustics within the framework of secondary infall model of formation of dark matter haloes. The result is then used to demonstrate that caustics can be probed by properly stacking the weak-lensing signal of about 600 haloes. CFHTLS accompanied by X-ray observations and the space-based experiments like SNAP or DUNE can provide us with the required statistics and hence a way of distinguishing between the viable dark matter particle candidates. The extension of our results to more realistic models including the effects of mergers of haloes is briefly outlined.

We classify orbits of stars that are bound to central black holes in galactic nuclei. The stars move under the combined gravitational influences of the black hole and the central star cluster. Within the sphere of influence of the black hole, the orbital periods of the stars are much shorter than the periods of precession. We average over the orbital motion and end up with a simpler problem and an extra integral of motion: the product of the black hole mass and the semimajor axis of the orbit. Thus the black hole enforces some degree of regularity in its neighborhood. Well within the sphere of influence, (i) planar, as well as three dimensional, axisymmetric configurations-both of which could be lopsided-are integrable, (ii) fully three dimensional clusters with no spatial symmetry whatsover must have semi-regular dynamics with two integrals of motion. Similar considerations apply to stellar orbits when the black hole grows adiabatically. We introduce a family of planar, non-axisymmetric potential perturbations, and study the orbital structure for the harmonic case in some detail. In the centered potentials there are essentially two main families of orbits: the familiar loops and lenses, which were discussed in Sridhar and Touma (1997, MNRAS, 287, L1-L4). We study the effect of lopsidedness, and identify a family of loop orbits, whose orientation reinforces the lopsidedness, an encouraging sign for the construction of self-consistent models of eccentric, discs around black holes, such as in M31 and NGC 4486B.

We construct a simple symplectic map to study the dynamics of eccentric orbits in non-spherical potentials. The map offers a dramatic improvement in speed over traditional integration methods, while accurately representing the qualitative details of the dynamics. We focus attention on planar, non-axisymmetric power-law potentials, in particular the logarithmic potential. We confirm the presence of resonant orbit families (``boxlets'') in this potential and uncover new dynamics such as the emergence of a stochastic web in nearly axisymmetric logarithmic potentials. The map can also be applied to triaxial, lopsided, non-power-law and rotating potentials.

We construct three dimensional axisymmetric, cuspy density distributions, whose potentials are of St\"ackel form in parabolic coordinates. As in Sridhar and Touma (1997), a black hole of arbitrary mass may be added at the centre, without destroying the St\"ackel form of the potentials. The construction uses a classic method, originally due to Kuzmin (1956), which is here extended to parabolic coordinates. The models are highly oblate, and the cusps are "weak", with the density, $\rho \propto 1/r^k$, where $0<k<1$.

We present cuspy, non-axisymmetric, scale-free mass models of discs, whose gravitational potentials are of St\"ackel form in parabolic coordinates. A black hole may be added at the centre, without in any way affecting the St\"ackel form; the dynamics in these potentials is, of course, fully integrable. The surface density, $\Sigma_{disc}\propto 1/r^{\gamma}$, where $0 <\gamma < 1$ corresponds to steep cusps for which the central force diverges. Thus cusps, black holes, and non-axisymmetry are not a sure recipe for chaos, as is generally assumed. A new family of orbits, lens orbits, emerges to replace the box orbits of models of elliptical galaxies that have constant-density cores. Loop orbits are conspicuous by their absence. Both lenses and boxlets (the other family of orbits), can be elongated in the direction of the density distribution, a property that is favourable for the construction of non-axisymmetric, self-consistent equilibrium models of elliptical galaxies.