The statistical mechanics of selfgravitating Keplerian disks
Abstract
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 selfgravity of the disk dominates the nonKeplerian force, and the spread in semimajor axes is small. We show that with plausible approximations such disks have logarithmic twobody 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 wellbehaved. We solve for the microcanonical axisymmetric thermal equilibria and demonstrate the existence of a symmetrybreaking 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 longrange interactions, such as point vortices in twodimensional fluids.
1 Introduction
The thermodynamics of isolated, bound, selfgravitating stellar systems (systems of point particles interacting only through gravitational forces) is notoriously pathological. The infinite volume of physical space forbids maximumentropy states and the shortrange 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 Nbody system in a spherical box, and/or regularizing or “softening” the shortrange singularity. With these artifices, canonical and microcanonical equilibria of the selfgravitating gas can be constructed (LyndenBell & Wood, 1968; Thirring, 1970; LyndenBell & LyndenBell, 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; LyndenBell & LyndenBell 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 SelfGravitating 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 selfgravitating systems in the real world. What remains at the end of the day are robust results on the thermodynamics of artificially imprisoned and mutilated selfgravitating 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 selfgravitating 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 longrange 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 )^{1}^{1}1It might seem more natural to model a disk containing only prograde orbits. However, the orbitaveraged 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 timescales of order times the orbital period. Relaxation can be studied by averaging over the fast orbital timescale (the orbital period), that is, replacing each particle by a socalled 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 semimajor axes do not. Since each wire has a constant semimajor 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 periapsis^{2}^{2}2For retrograde particles, our definition of differs from the usual convention (because is always measured counterclockwise 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 orbitaveraged 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 (semimajor 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 coordinatemomentum 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) semimajor axis . Such a limit is reasonable in disks where the fractional spread in semimajor 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 semimajor axes.
Last but not least, we require the orbitaveraged 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 potential^{3}^{3}3An 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 phasespace area elements, , to write . We define a dimensionless meanfield potential of the disk by
(1) 
with . This potential is the meanfield Hamiltonian, in the sense that (see eq. A12)
(2) 
The disk’s entropy is defined by
(3) 
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
(4) 
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
(5) 
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 :
(6) 
Throughout this paper we shall approximate the potential by the logarithmic potential , and since the integral equation can be replaced by a differential one^{4}^{4}4Apart from the factor , when this is the equation for the selfgravitating isothermal cylinder (Stodólkiewicz, 1963; Ostriker, 1964; Katz & LyndenBell, 1978; Aly, 1994).,
(7) 
Rather than total angular momentum or energy, we shall work with the dimensionless quantities
(8) 
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 .
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 angularmomentum sequence terminates at a point marked by a cross. Sequences of models with nonzero 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 loweccentricity analytic limit (eqs. B5–B7), 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 nonzero angular momentum the sequence terminates at finite entropy.
We now investigate the response of these equilibria to small nonaxisymmetric 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 nonaxisymmetric 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 line^{5}^{5}5The 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 firstorder phase transition at zero angular momentum, which transitions into a secondorder 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 & LyndenBell 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.
We next construct nonaxisymmetric 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 nonaxisymmetry,
(9) 
where , etc.^{6}^{6}6It 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 semimajor axis and unit mass. Thus for axisymmetric disks and for a disk in which all the eccentricity vectors have and are aligned or antialigned. The bottom left panel shows the pattern speed (B12) for the nonaxisymmetric disks, defined to be those with .
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 nonaxisymmetric 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 nonaxisymmetric sequence should extend to .
The behavior of models with angular momentum (blue triangles) is qualitatively similar, in that the axisymmetric sequence bifurcates to a nonaxisymmetric sequence as the energy is decreased (at energy and mean eccentricity ). However, neither the axisymmetric nor the nonaxisymmetric 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 nonaxisymmetric 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 nonaxisymmetric sequence. For no bifurcation is expected or observed.
4 Thermal and dynamical stability
We come now to the relation between thermal and dynamical instability. In the orbitaveraged dynamics described here, dynamical instability typically proceeds on the secular timescale, 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 timescale of the orbital period; on this timescale the disks should be stable since their mass is much smaller than the central mass. Thermal instability proceeds on the resonant relaxation timescale, 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 phasespace density and the thermal instability may not. For similar reasons, a dynamically unstable initial state does not normally evolve towards a maximumentropy final state on the secular time. Thus, we expect that dynamical instability leads in a timescale to an “intermediate” state that is a timeindependent solution of the collisionless Boltzmann equation, and that the intermediate state then evolves on a timescale to the maximumentropy state.
Analyses of the dynamical stability of collisionless nearKeplerian stellar disks, with or without a range of semimajor axes (Touma, 2002; Sridhar & Saini, 2010; Kazandjian & Touma, 2013), generally find that if there is a sufficient number of counterrotating particles (sufficiently small total angular momentum) the disks are dynamically unstable and settle into lopsided states on a secular timescale. 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 nonaxisymmetric maximumentropy 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 Nwire simulations in analogy to Nbody 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 maximumentropy 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 maximumentropy 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 maximumentropy state with the same energy and angular momentum ( and ). Similarly, the nonaxisymmetry parameter (eq. 9) fluctuates around 1.0 in the simulation, close to but 10% larger than its value of 0.91 in the maximumentropy state. In the top right panel of Figure 3, we superpose the eccentricity vectors of the 256 wires in the Nwire simulation (circles) at onto a 256point sample of the eccentricity vectors in the maximumentropy state with the same energy and angular momentum (magenta crosses). The distributions are similar, but the mean eccentricity vector of the maximumentropy state is smaller and its spread around the mean is broader.
We have carried out Nwire 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 maximumentropy states of Figure 2. The case is more complex. As expected, the Nwire 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 maximumentropy 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 maximumentropy 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 maximumentropy configuration; in general states with lower mean eccentricity have higher pattern speeds and vice versa. A larger Nwire simulation () showed similar transitions over the same timescale 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 maximumentropy state (). The pattern speed settled after a series of ups and downs to a mean value , close to the value expected in the maximumentropy state with the same energy and angular momentum. Models with revealed in dynamical simulations some of the same pathologies displayed by their counterparts in the search for maximum entropy nonaxisymmetric states: in particular, states that are predicted to go unstable seemed stuck in the neighborhood of their initial nearequilibrium 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 maximumentropy 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 maximumentropy 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 nonaxisymmetric entropy extrema (compare the lowerright panels of Figures 1 and 2).
5 Discussion
We have examined the maximumentropy states of a razorthin 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 orbitaveraged value. This approximation is appropriate if the disk age is shorter than the timescale for twobody relaxation due to close encounters. The orbitaveraged interaction between particles leads to resonant relaxation, in which the angular momenta and eccentricities of the particles relax, but the semimajor axes remain fixed. For simplicity, we focus in this paper on the somewhat artificial case in which all the particles have the same semimajor 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 selfgravitating 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 nearKeplerian stellar system, with the following properties (taken from Kocsis & Tremaine 2011): blackhole mass ; number of stars within 0.1 pc ; orbital period at 0.1 pc ; age ; resonant relaxation time .
We construct the maximumentropy equilibria that should be the endstate 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 maximumentropy state has a minimum mean eccentricity (top left panel of Figure 2) which is achieved at a critical value of the energy, . The maximumentropy 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, .
The results presented in this paper are based on a logarithmic approximation to the orbitaveraged potential energy between two particles, an approximation that is valid only in the limit of small eccentricities. The mean eccentricities of the nonaxisymmetric equilibria are large enough to cast doubt on the validity of this approximation. However, we have repeated our calculations using the exact orbitaveraged potential (computed on a threedimensional grid in –) and we found that the maximumentropy states produced with the logarithmic potential and the exact potential have all the same qualitative features (bifurcation points, minimum mean eccentricity, lopsided equilibria, etc.). Maximumentropy 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 4096point 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 Markovchain Monte Carlo simulations, basisfunction 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 timesscale, that is, a timesscale 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 Nwire 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 maximumentropy 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)^{7}^{7}7In both Nwire and Nbody simulations of unstable counterrotating disks (Touma et al., 2009; Kazandjian & Touma, 2013), stellar orbits experience largeamplitude oscillations in inclination when their eccentricity increases beyond a critical value. Such eccentricityinclination 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 threedimensional maximumentropy 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 maximumentropy state resulting from thermal instability, since this is not generally true in selfgravitating systems (e.g., in collapse of spherical systems that are not initially in virial equilibrium, where there is no maximumentropy state)—perhaps part of the answer is that the phase space of the systems examined here is compact.
Our results on the equilibria of selfgravitating systems with logarithmic twobody 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 semimajor axes in the secular dynamics of wires. Of course there are obvious and important differences: in selfgravitating 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; negativetemperature states in vortices are prone to phase transitions, whereas negativetemperature 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 nearKeplerian stellar disks. These results are of interest both for exploring the thermodynamics of systems with longrange forces and because they suggest that many nearKeplerian, nearly collisionless, astrophysical disks (disks near supermassive black holes, debris disks around young stars, etc.) may naturally develop a lopsided configuration.
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 timeaveraged gravitational interaction energy between two particles:
(A1) 
where and are the semimajor 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 semimajor axis . Then
(A2) 
where
(A3)  
here and are the true anomaly and longitude of periapsis of particle (so ), and
(A4) 
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,
(A5) 
When , the potential can be expanded in powers and logarithms of ,
(A6)  
The integral for diverges logarithmically as , suggesting that the potential can be written in the form
(A7) 
where and are smooth functions. The following functional forms fit the potential with an rms fractional error of 2%:
(A8) 
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 potential^{8}^{8}8These 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
(A9) 
The Poincaré variables are a canonical coordinatemomentum pair, that is, they evolve at a rate
(A10) 
For prograde particles is the coordinate and is the momentum, while their roles are reversed for retrograde particles. We rescale variables, time and Hamiltonian:
(A11) 
In these variables, Hamilton’s equations read
(A12) 
The analogous equations for are
(A13) 
There are two conserved quantities: the energy and the angular momentum
(A14) 
The phasespace area element is
(A15) 
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
(A16) 
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
(A17) 
The equations of motion can be rewritten in the compact vectorial form
(A18) 
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 Helmholtzlike 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 guidingcenter 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
(A19) 
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 phasespace 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 meanfield potential of the maximumentropy disk (eq. 7) simplifies to the ordinary differential equation
(B1) 
To solve this, we write , with , and define . The differential equation (B1) becomes
(B2) 
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),
(B3) 
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
(B4) which has the solution (Stodólkiewicz, 1963; Ostriker, 1964)
(B5) The validity of this approximate solution requires and . The dimensionless angular momentum and energy are
(B6) The parameter and the specific entropy is
(B7) 
In disks with (zero inverse temperature) the dimensionless angular momentum is related to the parameter by
(B8) the specific entropy is
(B9) and the mean eccentricity and prograde fraction are
(B10) where is a modified Bessel function.
Numerical solutions of the differential equation (B2) are discussed in the main text.
b.2 Bifurcation to nonaxisymmetric disks
We may use the differential equation (7) to investigate whether the equilibrium axisymmetric disks can remain in equilibrium under small nonaxisymmetric 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
(B11) 
The existence of a solution satisfying the boundary conditions as and as implies a bifurcation to a sequence of nonaxisymmetric disks that initially have fold symmetry.
Nonaxisymmetric 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, nonaxisymmetric 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 nonKeplerian energy and the angular momentum of a single particle. In dimensionless variables . Comparison to equation (4) implies that the dimensionless pattern speed is
(B12) 
b.3 Dynamical stability
In terms of the Poincaré eccentricity and the argument of periapsis , the equations of motion (A12) read
(B13) 
The distribution function must satisfy the collisionless Boltzmann equation
(B14) 
At this point we switch from the Poincaré eccentricity to the ordinary eccentricity :
(B15) 
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
(B16) 
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
(B17) 
Thus we arrive at an eigenvalue equation for the frequency ,