The statistical mechanics of self-gravitating Keplerian disks

The statistical mechanics of self-gravitating Keplerian disks

Jihad Touma11affiliation: Department of Physics, American University of Beirut, PO Box 11–0236, Riad El-Solh, Beirut 1107 2020, Lebanon; 22affiliation: School of Natural Sciences, Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA; and Scott Tremaine22affiliation: School of Natural Sciences, Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA;

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.

1 Introduction

The thermodynamics of isolated, bound, self-gravitating stellar systems (systems of point particles interacting only through gravitational forces) is notoriously pathological. The infinite volume of physical space forbids maximum-entropy states and the short-range singularity in the potential makes for an unbounded energy, hence undermining the construction of microcanonical ensembles (for reviews see Padmanabhan 1990 and Katz 2003). To make progress, the usual approach is to introduce artificial cutoffs by confining the N-body system in a spherical box, and/or regularizing or “softening” the short-range singularity. With these artifices, canonical and microcanonical equilibria of the self-gravitating gas can be constructed (Lynden-Bell & Wood, 1968; Thirring, 1970; Lynden-Bell & Lynden-Bell, 1977). There are regimes with negative heat capacity in the microcanonical equilibria, and these are associated with phase transitions of the corresponding canonical equilibrium (Thirring 1970; Lynden-Bell & Lynden-Bell 1977; Katz 1978; see Chavanis 2006 for a relatively recent review). Even simpler toy models are constructed with a view to isolating and analyzing what is thought to be generic behavior, e.g., the Hamiltonian Mean Field (HMF) (Antoni & Ruffo, 1995), the Self-Gravitating Ring (SGR) model (Sota et al., 2001), slab models in which point particles are replaced by infinite sheets (Rybicki, 1971; Joyce & Worrakitpoonpon, 2010; Schulz et al., 2013), and cylindrical models in which point particles are replaced by infinite wires (Stodólkiewicz, 1963; Ostriker, 1964; Aly, 1994; Aly & Perez, 1999). A large body of literature has evolved around these models, studying their equilibria, phase transitions, dynamical stability and metastability, and connections to the actual systems they are meant to model (e.g., Campa et al., 2009). Although instructive and elegant, these models leave one with the nagging question of what all of this has to do with actual self-gravitating systems in the real world. What remains at the end of the day are robust results on the thermodynamics of artificially imprisoned and mutilated self-gravitating systems, more tentative and largely numerical results on the evolution of realistic systems (Binney & Tremaine, 2008), and heuristic rules relating the properties of the former to the latter.

2 The Keplerian disk and ring

Here, we bridge the gap between tractable and realistic self-gravitating systems by examining a system that arises naturally (in the study of protoplanetary disks, stellar disks around supermassive black holes, etc.), and in which the pathologies of systems with long-range interactions are naturally resolved. We start with an infinitesimally thin disk composed of identical point particles, each of mass . The particles orbit a central point mass . The disk is flat; nevertheless particles may orbit in the prograde or retrograde direction (i.e., the inclinations are 0 or )111It might seem more natural to model a disk containing only prograde orbits. However, the orbit-averaged gravitational torque on an eccentric orbit does not approach zero as the eccentricity approaches unity. Thus if the phase space is restricted to prograde orbits, there will be a loss of particles through the boundary at zero angular momentum or .. Since , the particle orbits are nearly Keplerian. In such disks the dominant relaxation process is resonant relaxation (Rauch & Tremaine, 1996), involving secular interactions between particles that cause the orbits to evolve on time-scales of order times the orbital period. Relaxation can be studied by averaging over the fast orbital time-scale (the orbital period), that is, replacing each particle by a so-called Gaussian wire (a closed wire following the Keplerian orbit, with linear mass density inversely proportional to velocity). In interactions of this kind the angular momenta or eccentricities of the wires relax, but their energies or semi-major axes do not. Since each wire has a constant semi-major axis, it is completely specified by its mass , sense of rotation ( for prograde and for retrograde), eccentricity , and azimuth of periapsis , or instead of the last two the eccentricity vector , which points towards periapsis222For retrograde particles, our definition of differs from the usual convention (because is always measured counter-clockwise from the origin of azimuth rather than in the direction of orbital motion); the advantage of our convention is that prograde and retrograde orbits with the same eccentricity vector occupy the same locus in space.. The eccentricity vector rotates slowly due to the orbit-averaged force field of the other wires, and varies stochastically due to the even slower effect of resonant relaxation. Note that the conservation of Keplerian energy (semi-major axis) leaves us with a compact phase space for the wires to relax in, hence removing the need for artificial confinement.

Alternatively, a particle orbit can be specified by the Poincaré variables ; these are canonical coordinate-momentum pairs when multiplied by (see Appendix A.2). Note that and both range from 0 to 1, with as and as . We shall sometimes call the Poincaré eccentricity, and shift between eccentricity and Poincaré eccentricity as needed to keep the formulae as simple as possible. In the models described in this paper, which we call Keplerian rings, all particles are further assumed to share a common (and conserved) semi-major axis . Such a limit is reasonable in disks where the fractional spread in semi-major axes is smaller than the typical orbital eccentricity, but is chosen here mostly for simplicity, as the methods we describe are applicable to disks with any distribution of semi-major axes.

Last but not least, we require the orbit-averaged gravitational potential energy between two particles in the disk: , where denotes a time average over both orbits (see Appendix A.1). When eccentricities are small, the averaged potential can be evaluated analytically (Borderies et al., 1983), plus terms that are . For eccentricities that are not small, must be evaluated numerically by a double integral over the two orbital phases. Most of the calculations described below have been carried out both with the logarithmic potential and with a numerical evaluation of on a grid, and the main conclusions are qualitatively and quantitatively unaffected. Therefore for simplicity we present mostly the results with the logarithmic potential333An alternative approximation is that the potential is logarithmic in the distance between the Poincaré eccentricities, . We have experimented with this approximation and find that the rich behavior described here—bifurcation points, lopsided equilibria, etc.—is present with the exact potential and the approximate potential but not with ., except for a brief discussion associated with Figure 4.

In the continuum limit, let be the number of prograde or retrograde particles on orbits in the eccentricity range or . The total number of prograde and retrograde particles is or . In transforming between these we use the relation between phase-space area elements, , to write . We define a dimensionless mean-field potential of the disk by


with . This potential is the mean-field Hamiltonian, in the sense that (see eq. A12)


The disk’s entropy is defined by


Our aim is to extremize the entropy subject to the conservation of the number of particles ; the energy ; and the angular momentum . We denote the resulting distribution functions and potential and . Using Lagrange multipliers this optimization problem can be solved to give


where , , are dimensionless constants and as usual for prograde or retrograde particles. We must have (the distribution function cannot be negative). The parameter is an inverse temperature, which can be either positive or negative since the phase space is compact. Setting , Poisson’s equation (1) may now be written


with . This is a nonlinear integral equation for the dimensionless potential , whose solution depends on the parameters and . The inverse temperature is determined from the solution of (5) by substituting equation (4) into the relation :


Throughout this paper we shall approximate the potential by the logarithmic potential , and since the integral equation can be replaced by a differential one444Apart from the factor , when this is the equation for the self-gravitating isothermal cylinder (Stodólkiewicz, 1963; Ostriker, 1964; Katz & Lynden-Bell, 1978; Aly, 1994).,


Rather than total angular momentum or energy, we shall work with the dimensionless quantities


where . These, together with the integral equation (5), determine the potential and the parameters and , given the conserved quantities and ; thus, all thermodynamic equilibria can be parametrized by their dimensionless energy and angular momentum. The dimensionless energy cannot exceed , corresponding to particles uniformly distributed on the circle (see Appendix B.1); the absolute value of the dimensionless angular momentum cannot exceed unity, and without loss of generality we can restrict to the range .

States that are entropy extrema according to equation (7) can be either axisymmetric (i.e., depending on only through ) or non-axisymmetric. If they are non-axisymmetric the figure is stationary in a frame rotating at the pattern speed given by equation (B12).

3 Thermodynamics of the Keplerian ring

We first study axisymmetric entropy extrema, which we construct by solving the differential equation (7) (see Appendix B.1). We plot the results in Figure 1. The four colored curves show solutions with dimensionless angular momentum , as functions of the dimensionless energy . The four panels show the mean eccentricity , fraction of prograde particles, inverse temperature , and entropy (the last of these is for the normalization ; more generally ). Each constant angular-momentum sequence terminates at a point marked by a cross. Sequences of models with non-zero angular momentum terminate at a mean eccentricity less than unity (as they must, since orbits with have zero angular momentum). The small open triangles in the left part of each panel show the predictions of the low-eccentricity analytic limit (eqs. B5B7), which agree well with the numerical solutions. Remarkably, the curves of mean eccentricity versus energy (top left panel) almost coincide for the whole range of angular momenta shown.

The axisymmetric entropy extrema become increasingly prograde with increasing energy and mean eccentricity, as they should to maintain a constant angular momentum. The family of axisymmetric solutions includes regions of negative heat capacity () and negative temperature . The sequence of models with zero angular momentum has as ; in this limit the distribution function approaches a singular form in which all the particles have . For a given non-zero angular momentum the sequence terminates at finite entropy.

We now investigate the response of these equilibria to small non-axisymmetric perturbations, , , where defines the potential of the unperturbed axisymmetric system (see Appendices B.2 and B.4). We substitute this form into the differential equation (7) and linearize in the small parameter . The existence of a solution to the linearized equation implies a bifurcation to a sequence of non-axisymmetric disks that initially have -fold symmetry. We find numerically that (i) bifurcations exist for only; (ii) there is one and only one bifurcation point along the sequence of axisymmetric equilibria at fixed angular momentum for (see derivation at the end of Appendix B.4), and none for ; (iii) these bifurcations are associated with a transition from entropy maxima, hence thermally stable equilibria, to entropy saddle points which are thermally unstable. In Figure 1, we distinguish the regions in which each sequence is stable or unstable by solid and dotted lines, respectively, and mark the locus of bifurcation points by a heavy solid line555The properties of the system at the bifurcation are continuous functions of the energy in a microcanonical setting. Our preliminary exploration of the thermodynamics of our model disks in the canonical ensemble reveals a richer behavior, including a first-order phase transition at zero angular momentum, which transitions into a second-order transition with increasing angular momentum. The study of the canonical ensemble will take us too far afield in an already lengthy exploration of the microcanonical states and is relegated to future work.. The axisymmetric systems are thermally unstable at small mean eccentricity and stable at large mean eccentricity. This result is surprising, since in the limit of small mean eccentricity the equilibrium disks are identical to the isothermal cylinder (eq. B4), which is known to be an entropy maximum, and therefore stable (Katz & Lynden-Bell 1978; Aly 1994; see Aly & Perez 1999 for a generalization to unbounded systems with angular momentum constraint). The explanation is that any isolated system such as the isothermal cylinder is neutrally stable to displacements—in other words, the differential equation governing the density distribution is autonomous—whereas the differential equation (7) governing the eccentricity distribution contains terms involving that can make the neutral mode slightly unstable, no matter how small the mean eccentricity.

Figure 1: The properties of axisymmetric disks with a logarithmic potential. Each panel shows systems with dimensionless angular momentum (magenta), 0.5 (blue), 0.8 (green), and 0.95 (red). The four panels show the mean eccentricity (top left), fraction of prograde particles (top right), inverse temperature (bottom left), and entropy when (bottom right). Also shown as open triangles are the analytic predictions for low-eccentricity disks. Each sequence of models terminates at the point marked by a cross. Unstable parts of the constant angular-momentum sequences are represented by dotted lines, while stable parts are shown by solid lines. The locus of bifurcation points, which separates the stable and unstable regions, is shown by a heavy solid line (except in the top left panel, to avoid obscuring the equilibrium sequences with which it almost coincides). Sequences with have no bifurcation.

We next construct non-axisymmetric disks using the nonlinear optimization methods described in Appendix C. The results are shown in Figure 2. The top left and bottom right panels show the same quantities as in Figure 1. The top right panel shows a measure of the strength of the non-axisymmetry,


where , etc.666It is straightforward to show that is the difference between the larger and smaller of the two principal moments of inertia of the disk when the disk has unit semi-major axis and unit mass. Thus for axisymmetric disks and for a disk in which all the eccentricity vectors have and are aligned or anti-aligned. The bottom left panel shows the pattern speed (B12) for the non-axisymmetric disks, defined to be those with .

Figure 2: The properties of axisymmetric and non-axisymmetric disks with a logarithmic potential. Each panel shows systems with dimensionless angular momentum (magenta), 0.5 (blue), 0.8 (green), and 0.95 (red). The panels show mean eccentricity (top left), difference (eq. 9) between the major and minor axes of the inertia ellipse (top right), pattern speed (bottom left) and entropy when (bottom right). The triangles show the maximum-entropy states computed via nonlinear optimization; closed triangles are axisymmetric () and open triangles are non-axisymmetric. Crosses denote the termination of the axisymmetric sequences and solid circles denote bifurcation points, both taken from Fig. 1.

First consider the magenta triangles, which outline the locus of models with zero angular momentum and equal fractions of prograde and retrograde particles. The models begin near mean eccentricity and energy , corresponding to an axisymmetric disk composed of radial orbits. As the energy is reduced the models initially follow the axisymmetric sequence shown in Figure 1. At the bifurcation point, marked by a solid circle at , , the models leave the axisymmetric sequence, which is no longer an entropy maximum beyond this point, to follow a non-axisymmetric sequence with growing mean eccentricity. This sequence terminates in a disk composed of particles on radial orbits with aligned eccentricity vectors, and . The numerical models terminate at due to the limited resolution of our grid but we believe that the non-axisymmetric sequence should extend to .

The behavior of models with angular momentum (blue triangles) is qualitatively similar, in that the axisymmetric sequence bifurcates to a non-axisymmetric sequence as the energy is decreased (at energy and mean eccentricity ). However, neither the axisymmetric nor the non-axisymmetric sequence can achieve since such a system would be composed entirely of radial orbits, which have zero angular momentum. Instead, as the energy becomes more negative, the orbits cluster more and more tightly around a single value of the eccentricity vector, with magnitude given by . The numerical models terminate at but we believe this is because of the limited resolution of our grid, and the sequence should asymptote to a horizontal line at that extends to .

For (green triangles), our earlier analysis of the thermal stability of axisymmetric systems implies that there is a bifurcation to a non-axisymmetric sequence at , , but the numerical models remain on the axisymmetric sequence for all energies. This is presumably an artifact of our limited resolution, since (i) a bifurcation point exists only for , which is close to ; (ii) the entropy curves in the bottom right panel of Figure 3 are very close together once so it is difficult for the optimization code to settle onto the non-axisymmetric sequence. For no bifurcation is expected or observed.

Figure 3: Top left: the evolution of the mean eccentricity (solid) and the norm of the mean eccentricity vector (dashed) in a 256-wire dynamical simulation of the unstable axisymmetric equilibrium at and . Top right: a sample of the non-axisymmetric maximum-entropy equilibrium with this energy and angular momentum (magenta crosses) superposed with the dynamical simulation (black circles) around . The directions of the mean eccentricity vectors of the two states were reoriented to coincide. Bottom left: the evolution of the mean eccentricity of a 128-wire sample of an unstable axisymmetric equilibrium with and (the bifurcation energy at is ). Bottom right: the final mean eccentricity of simulated ensembles (with error bars), superimposed on a zoom-in of the top left panel of Fig. 2. All initial states that are expected to be dynamically stable stayed close to their initial axisymmetric state for the duration of the simulation. All initial states that were expected to be unstable (except for the one at , ) became lopsided with a final mean eccentricity close to that of the maximum-entropy states in Fig. 2.

4 Thermal and dynamical stability

We come now to the relation between thermal and dynamical instability. In the orbit-averaged dynamics described here, dynamical instability typically proceeds on the secular time-scale, which is longer than the orbital period by a factor (i.e., in the notation of eq. 2). We do not consider possible dynamical instabilities on the time-scale of the orbital period; on this time-scale the disks should be stable since their mass is much smaller than the central mass. Thermal instability proceeds on the resonant relaxation time-scale, which is expected to be longer than the secular time by a factor (Rauch & Tremaine, 1996). Thermal stability implies dynamical stability, but thermal instability need not imply dynamical instability (Bartholomew, 1971; Ipser & Horwitz, 1979) since the collisionless Boltzmann equation conserves phase-space density and the thermal instability may not. For similar reasons, a dynamically unstable initial state does not normally evolve towards a maximum-entropy final state on the secular time. Thus, we expect that dynamical instability leads in a timescale to an “intermediate” state that is a time-independent solution of the collisionless Boltzmann equation, and that the intermediate state then evolves on a timescale to the maximum-entropy state.

Analyses of the dynamical stability of collisionless near-Keplerian stellar disks, with or without a range of semi-major axes (Touma, 2002; Sridhar & Saini, 2010; Kazandjian & Touma, 2013), generally find that if there is a sufficient number of counter-rotating particles (sufficiently small total angular momentum) the disks are dynamically unstable and settle into lopsided states on a secular time-scale. To determine whether these conclusions apply to the disks studied in this paper, we have solved the linearized collisionless Boltzmann equation for the axisymmetric models shown in Figure 1 (see Appendix B.3). We find that dynamical instabilities are present in some models, and the onset of dynamical instability occurs at the same bifurcation points at which the disk becomes thermally unstable and the sequence of non-axisymmetric maximum-entropy models begins (to within 0.3% in energy ). In other words, it appears that the axisymmetric models are dynamically unstable if and only if they are thermally unstable.

To explore further the relation between dynamical and thermal instability in these systems, we simulated the dynamical evolution of ensembles of Gaussian wires selected from the distribution functions of axisymmetric thermal equilibria. We call these N-wire simulations in analogy to N-body simulations (Touma et al., 2009). Ensembles with 128, 256 and 512 wires were simulated at and , at energies above and below the bifurcation points identified in Figure 1. In the top left panel of Figure 3, we display the mean eccentricity and norm of the mean eccentricity vector for an initially axisymmetric system with and , which is thermally unstable according to Figures 1 and 2. The mean eccentricity shows a rapid departure from its equilibrium value in the axisymmetric system () through an (overstable) cycle, which saturates after a sequence of oscillations of gradually decreasing amplitude at a mean eccentricity . The mean eccentricity vector follows suite, departing from and saturating in a lopsided configuration with . We identify these configurations with the “intermediate” equilibria described above. On timescales we expect that the intermediate equilibria should evolve toward maximum-entropy equilibria. We have not been able to detect this evolution, simply because the macroscopic properties of the intermediate equilibria are already close to those of the maximum-entropy equilibria when they first appear. For example, the mean eccentricity and mean eccentricity vectors in the intermediate state at –40 ( and ) are within a few percent of the corresponding quantities in the maximum-entropy state with the same energy and angular momentum ( and ). Similarly, the non-axisymmetry parameter (eq. 9) fluctuates around 1.0 in the simulation, close to but 10% larger than its value of 0.91 in the maximum-entropy state. In the top right panel of Figure 3, we superpose the eccentricity vectors of the 256 wires in the N-wire simulation (circles) at onto a 256-point sample of the eccentricity vectors in the maximum-entropy state with the same energy and angular momentum (magenta crosses). The distributions are similar, but the mean eccentricity vector of the maximum-entropy state is smaller and its spread around the mean is broader.

We have carried out N-wire simulations of zero angular momentum () axisymmetric equilibria at other energies, both below and above the bifurcation value , and these were equally robust in converging in the mean to states close to the expected maximum-entropy states of Figure 2. The case is more complex. As expected, the N-wire simulations showed stability and instability for values of the energy larger and smaller, respectively, than the bifurcation energy . However the dynamical evolution was far more tortuous. In the bottom left panel of Figure 3, we follow the mean eccentricity of an ensemble of 128 wires sampling an initially axisymmetric (and thermally unstable) equilibrium with and . The cluster transitions rather fast to a lopsided state with mean eccentricity then undergoes a further transition around to a more lopsided state with , almost exactly the value in the maximum-entropy state (). In both states the mean eccentricity exhibits fluctuations with an amplitude of about 0.08. By , the lopsided system is precessing with a mean pattern speed , close to the value expected in the maximum-entropy state with the same energy and angular momentum. The evolution over nearly 200 mode precession periods () shows a number of intermittent transitions to states with lower mean eccentricity and few signs of settling down to a maximum-entropy configuration; in general states with lower mean eccentricity have higher pattern speeds and vice versa. A larger N-wire simulation () showed similar transitions over the same time-scale so these are unlikely to be an artifact of small . A simulation with and , further from the bifurcation energy, lingers around the nearly axisymmetric initial state until about , before it undergoes a series of transitions to larger mean eccentricity, eventually (by ) attaining (with fluctuations of about 0.15), close to the mean eccentricity of the maximum-entropy state (). The pattern speed settled after a series of ups and downs to a mean value , close to the value expected in the maximum-entropy state with the same energy and angular momentum. Models with revealed in dynamical simulations some of the same pathologies displayed by their counter-parts in the search for maximum entropy non-axisymmetric states: in particular, states that are predicted to go unstable seemed stuck in the neighborhood of their initial near-equilibrium configuration, even in relatively lengthy simulations with wires. The final states of all of these simulations are displayed as solid squares with error bars in the bottom right panel of Figure 3), along with the maximum-entropy equilibria shown in the top left panel of Figure 1.

We conclude that in some of our models dynamical instability leads to “intermediate” states that are close to maximum-entropy states; other models, particuarly those with significant angular momentum, often seem to linger in, or oscillate between, metastable states. Possibly this behavior is associated with the small difference in entropy between the axisymmetric and non-axisymmetric entropy extrema (compare the lower-right panels of Figures 1 and 2).

5 Discussion

We have examined the maximum-entropy states of a razor-thin disk of collisionless masses orbiting a massive central body. The disks may contain particles on both prograde and retrograde orbits and particles are allowed to flip between prograde and retrograde orbits, but the total energy and angular momentum of the disk are conserved. The disk mass is assumed to be much smaller than the mass of the central body, so the interaction potential between two particles can be approximated by its orbit-averaged value. This approximation is appropriate if the disk age is shorter than the time-scale for two-body relaxation due to close encounters. The orbit-averaged interaction between particles leads to resonant relaxation, in which the angular momenta and eccentricities of the particles relax, but the semi-major axes remain fixed. For simplicity, we focus in this paper on the somewhat artificial case in which all the particles have the same semi-major axis (“Keplerian rings”), although our methods are easily adapted to more general disk models.

Although the Keplerian rings described here are artificial systems intended mainly as aids in exploring the dynamics and statistical mechanics of self-gravitating stellar systems, it is useful to relate them to the properties of a real astrophysical system to which they may offer insight. The center of the Milky Way galaxy contains a black hole surrounded by a near-Keplerian stellar system, with the following properties (taken from Kocsis & Tremaine 2011): black-hole mass ; number of stars within 0.1 pc ; orbital period at 0.1 pc ; age ; resonant relaxation time .

We construct the maximum-entropy equilibria that should be the end-state of resonant relaxation. The natural expectation is that such disks should be axisymmetric, with an eccentricity distribution given approximately by the analytic solution in (Stodólkiewicz, 1963; Ostriker, 1964), at least so long as the mean eccentricity is not too large. This expectation is not correct: for a given angular momentum we find that the maximum-entropy state has a minimum mean eccentricity (top left panel of Figure 2) which is achieved at a critical value of the energy, . The maximum-entropy state is axisymmetric for energy and lopsided for . For the axisymmetric equilibrium is an entropy extremum but not a maximum. Both the pattern speed and the temperature of the lopsided disks are generally positive. Essentially, as the disk is cooled to lower and lower energies the stellar orbits concentrate around a single eccentricity vector whose magnitude is determined by the angular momentum, .

Figure 4: As in Figure 2, except the gravitational potential is computed using the exact expression (A3) rather than the logarithmic approximation (A5).

The results presented in this paper are based on a logarithmic approximation to the orbit-averaged potential energy between two particles, an approximation that is valid only in the limit of small eccentricities. The mean eccentricities of the non-axisymmetric equilibria are large enough to cast doubt on the validity of this approximation. However, we have repeated our calculations using the exact orbit-averaged potential (computed on a three-dimensional grid in ) and we found that the maximum-entropy states produced with the logarithmic potential and the exact potential have all the same qualitative features (bifurcation points, minimum mean eccentricity, lopsided equilibria, etc.). Maximum-entropy models computed in this way are shown in Figure 4, which should be compared with Figure 2.

The numerical methods we have used need to be improved. At present we find the entropy maxima using sequential quadratic programming, defining the distribution function on a 4096-point grid in eccentricity space. In a few cases we find suspicious numerical artifacts (e.g., the small discontinuity near , in the top left panel of Figure 2), and in most cases convergence is quite slow. This said, we have confirmed our main results with Markov-chain Monte Carlo simulations, basis-function expansions of the central integral equation (see Appendix C.2), numerical solutions of the analogous differential equation (7) for axisymmetric states, and nonlinear optimization using smaller grids.

We have also examined possible dynamical instabilities in axisymmetric Keplerian rings, which are expected to occur on the secular times-scale, that is, a times-scale longer than the orbital period by the ratio of the central mass to the disk mass. We find that the rings are dynamically unstable if and only if they are thermally unstable. We showed via N-wire simulations that dynamical instability in these disks produces lopsided states. We observe that in some but not all of our experiments these are close to the maximum-entropy solutions in mean eccentricity, distribution of eccentricity vectors, and precession rates or pattern speeds; we do not have an explanation for this similarity nor do we know whether it has an illuminating physical explanation. The presence of this instability is in line with earlier findings of generic dynamical instabilities in disks containing a retrograde stellar population (Touma, 2002; Touma et al., 2009; Gulati et al., 2012)777In both N-wire and N-body simulations of unstable counter-rotating disks (Touma et al., 2009; Kazandjian & Touma, 2013), stellar orbits experience large-amplitude oscillations in inclination when their eccentricity increases beyond a critical value. Such eccentricity-inclination instabilities may operate in the disks we consider here if given the freedom to do so, and make it imperative to generalize our results to three-dimensional maximum-entropy equilibria. In addition to endowing our models with greater physical realism, the extra degree of freedom provides a natural way to resolve the otherwise singular transition from the prograde to the retrograde sector of phase space.. However, we do not know why the final state of the dynamical instability is so similar to the maximum-entropy state resulting from thermal instability, since this is not generally true in self-gravitating systems (e.g., in collapse of spherical systems that are not initially in virial equilibrium, where there is no maximum-entropy state)—perhaps part of the answer is that the phase space of the systems examined here is compact.

Our results on the equilibria of self-gravitating systems with logarithmic two-body potentials in eccentricity space have a strong kinship with the far more extensive body of work on the statistical mechanics of point vortices in compact domains (see Appendix A.3). The interaction potential for vortices is logarithmic in physical space, so physical space for vortices maps into eccentricity space for wires, and conserved circulation in point vortices to conserved semi-major axes in the secular dynamics of wires. Of course there are obvious and important differences: in self-gravitating wires the potential energy does not depend on the direction of motion—prograde or retrograde—of the particles, whereas it does depend on the sign of circulation of vortices; wires can evolve between prograde and retrograde, while vortices cannot change their circulation; negative-temperature states in vortices are prone to phase transitions, whereas negative-temperature rings appear perfectly stable in axisymmetric configurations. This said, much of the analytic machinery developed to study the existence and stability of solutions in the point vortex case should extend quite naturally to our problem.

We have found remarkable and unexpected complexity in the thermodynamics of near-Keplerian stellar disks. These results are of interest both for exploring the thermodynamics of systems with long-range forces and because they suggest that many near-Keplerian, nearly collisionless, astrophysical disks (disks near supermassive black holes, debris disks around young stars, etc.) may naturally develop a lopsided configuration.

This research was supported in part by NASA grant NNX11AF29G. JT acknowledges the support of an Arab Fund Research Fellowship for the year 2013–2014, which allowed him an extended stay at the IAS and a briefer one at the IHP (Paris), the hospitality of both institutes being greatly appreciated.

Appendix A The Keplerian ring

We assemble here the mathematical machinery, remarks, and results which underlie and amplify the results and assertions in the body of the text.

a.1 The interaction potential

The study of secular dynamics requires the time-averaged gravitational interaction energy between two particles:


where and are the semi-major axes and eccentricity vectors of the particles, and denotes a time average over both orbits. In this paper we examine the special case where all particles share the same semi-major axis . Then




here and are the true anomaly and longitude of periapsis of particle (so ), and


The function in equation (A3) is symmetric in its arguments, , and is rotationally invariant, that is, it depends on and only through . When the eccentricities are small, the integral can be evaluated analytically,


When , the potential can be expanded in powers and logarithms of ,


The integral for diverges logarithmically as , suggesting that the potential can be written in the form


where and are smooth functions. The following functional forms fit the potential with an rms fractional error of 2%:


The numerical experiments in this paper use the logarithmic potential (eq. A5) even for eccentricities of order unity, where it is not strictly valid. The reason for this is that we have also conducted experiments with the accurate averaged potential888These experiments do not use equation (A8) or other fitting formulae; instead they rely on numerical evaluations of the double integral (A3) on a uniform or grid in the space where is the Poincaré eccentricity., and found that these are in essential qualitative agreement with those obtained via the logarithmic potential (see Figure 4).

a.2 Hamilton’s equations

The Hamiltonian of a particle with eccentricity is


The Poincaré variables are a canonical coordinate-momentum pair, that is, they evolve at a rate


For prograde particles is the coordinate and is the momentum, while their roles are reversed for retrograde particles. We rescale variables, time and Hamiltonian:


In these variables, Hamilton’s equations read


The analogous equations for are


There are two conserved quantities: the energy and the angular momentum


The phase-space area element is


It is sometimes useful to think of the phase space as the surface of a sphere of unit radius, in which the azimuthal angle is and the polar angle is where


We call this the eccentricity sphere. The northern (southern) hemisphere represents prograde (retrograde) orbits, and the equator represents radial orbits (). The area element on the eccentricity sphere is proportional to the area element in phase space, .

a.3 Relation to point vortices

Equations (A13) are closely related to Kirchhoff’s equations for the motion of Helmholtz point vortices on the sphere (Kiessling & Wang, 2012). Let the unit vector denote the location of one of our particles on the eccentricity sphere. The rate of change of and is given by equations (A13). Taking the derivative of with respect to we get


The equations of motion can be rewritten in the compact vectorial form


While we believe that this formulation offers greater insight and elegance, and may well simplify the mathematical analysis of -wire systems, all the calculations in this paper are based on the (mathematically equivalent) geometry of a double copy of a disk of unit radius, corresponding to populations of prograde and retrograde wires.

Equation (A18) also governs the dynamics of a collection of point vortices all having the same circulation, evolving according to a Helmholtz-like Hamiltonian that is proportional to . The only difference between the Helmholtz Hamiltonian and ours (apart from a constant of proportionality) is that Helmholtz’s is proportional to whereas ours is proportional to (eq. A5).

More generally, our equations, methods and to some extent solutions have strong kinship with the vast body of literature dedicated to the statistical mechanics of point vortices (Onsager, 1949; Montgomery & Joyce, 1974; Kida, 1975; Miller, 1990; Robert & Sommeria, 1991), and of guiding-center plasmas (Joyce & Montgomery, 1973; Smith & O’Neil, 1990) in compact planar domains. The disks discussed in this paper provide, in their dynamics and thermodynamics, a direct bridge between collisionless stellar systems and systems of two dimensional point vortices, the formal analogies between them having been discussed at length in the literature (Chavanis et al., 1996; Chavanis, 2002).

a.4 Boundary conditions in the mean field limit

We consider boundary conditions on the density near the boundary at . To investigate these, we replace and in the equations of motion (A13) by the azimuthal and polar angles on the eccentricity sphere, and (eq. A16). Then


In general is smooth near so the quantities in brackets are also smooth. We conclude that and near . Thus the trajectories intersect the equator () along lines of constant longitude, at constant latitudinal speed.

Since the motion of particles near the equator is smooth, in a steady state the density of particles on the sphere should be smooth near the equator. Since the area element on the eccentricity sphere is proportional to the area element in a canonical phase space, the density of particles on the eccentricity sphere is proportional to the phase-space density in the northern hemisphere and in the south. Thus must join smoothly onto ; in particular at the equator. This in turn requires that as , for some function .

Appendix B Axisymmetric Keplerian rings and their perturbations

b.1 Equilibria

For axisymmetric disks, the nonlinear Poisson equation governing the mean-field potential of the maximum-entropy disk (eq. 7) simplifies to the ordinary differential equation


To solve this, we write , with , and define . The differential equation (B1) becomes


for given values of and , this can be solved by integrating outwards from with the initial conditions . The potential at the center is then given by equation (5),


Knowing and we can compute the parameter from as well as the dimensionless energy and angular momentum .

The differential equation (B2) does not appear to have a general analytic solution. However, some aspects of the behavior of these disks can be deduced analytically:

  • There is an upper limit to the dimensionless energy : the potential increases monotonically with the distance but particles are restricted to the circular area . Thus the energy of a distribution of particles with given mass is maximized if they are uniformly distributed on the circle , in which case it is straightforward to show that .

  • If the particles have small eccentricities we can replace by unity in equation (B2), to obtain


    which has the solution (Stodólkiewicz, 1963; Ostriker, 1964)


    The validity of this approximate solution requires and . The dimensionless angular momentum and energy are


    The parameter and the specific entropy is

  • In disks with (zero inverse temperature) the dimensionless angular momentum is related to the parameter by


    the specific entropy is


    and the mean eccentricity and prograde fraction are


    where is a modified Bessel function.

Numerical solutions of the differential equation (B2) are discussed in the main text.

b.2 Bifurcation to non-axisymmetric disks

We may use the differential equation (7) to investigate whether the equilibrium axisymmetric disks can remain in equilibrium under small non-axisymmetric perturbations. If the potential , where and define the potential of the unperturbed axisymmetric system defined following equation (B1), is an integer, and is small, equation (7) can be linearized to yield


The existence of a solution satisfying the boundary conditions as and as implies a bifurcation to a sequence of non-axisymmetric disks that initially have -fold symmetry.

Non-axisymmetric equilibria can be stationary in an inertial frame or in a frame precessing with some pattern speed, which we denote (relative to the physical time) or relative to the dimensionless time defined in equation (2). A rotating, non-axisymmetric equilibrium is a solution of the collisionless Boltzmann equation—which it must be, since the relaxation time is much longer than the orbital or precession time—if and only if the distribution function depends only on the Jacobi integral where and are the non-Keplerian energy and the angular momentum of a single particle. In dimensionless variables . Comparison to equation (4) implies that the dimensionless pattern speed is


b.3 Dynamical stability

In terms of the Poincaré eccentricity and the argument of periapsis , the equations of motion (A12) read


The distribution function must satisfy the collisionless Boltzmann equation


At this point we switch from the Poincaré eccentricity to the ordinary eccentricity :


We now write , , where , , and define the potential and distribution function of the unperturbed axisymmetric system, is an integer, and and are small. We then linearize equation (B15) to obtain


where , and the last equality follows from the definition (4) of the equilibrium distribution function and .

The perturbed potential and the perturbed distribution function are related by Poisson’s equation, which reads


Thus we arrive at an eigenvalue equation for the frequency ,