Isotropic-Nematic Phase Transitions in Gravitational Systems

TL;DR

This paper models dense stellar systems with a central black hole as liquid crystal-like phases, revealing phase transitions between ordered disk-like and disordered isotropic states based on angular momentum.

Contribution

It introduces a novel analogy between gravitational systems and liquid crystal phases, demonstrating phase transitions driven by angular momentum in such systems.

Findings

Existence of ordered and disordered phases in stellar systems

Identification of a first-order phase transition at a critical angular momentum

Discovery of metastable states with mutually inclined disks

Abstract

We examine dense self-gravitating stellar systems dominated by a central potential, such as nuclear star clusters hosting a central supermassive black hole. Different dynamical properties of these systems evolve on vastly different timescales. In particular, the orbital-plane orientations are typically driven into internal thermodynamic equilibrium by vector resonant relaxation before the orbital eccentricities or semimajor axes relax. We show that the statistical mechanics of such systems exhibit a striking resemblance to liquid crystals, with analogous ordered-nematic and disordered-isotropic phases. The ordered phase consists of bodies orbiting in a disk in both directions, with the disk thickness depending on temperature, while the disordered phase corresponds to a nearly isotropic distribution of the orbit normals. We show that below a critical value of the total angular momentum, the system undergoes a first-order phase transition between the ordered and disordered phases. At the critical point the phase transition becomes second-order while for higher angular momenta there is a smooth crossover. We also find metastable equilibria containing two identical disks with mutual inclinations between $90^{\circ}$ and $180^\circ$.

Equations

$$t_{\mathrm{orb}}\ll t_{\mathrm{in}-\mathrm{plane}}\ll t_{\mathrm{vrr}}\ll t_{\mathrm{srr}}\ll t_{2-\mathrm{body}}.$$

$$l=m\sqrt{G{M}_{\bullet }a(1-e^2)};$$

$$N_i=\int f_i(\mathbf{n})d\Omega ,$$

$$N=\sum _{i=1}^K{N}_i.$$

$$\mathbf{l}=l_i\mathbf{n}\text{mbox}where{l}_i=m_i\sqrt{G{M}_{\bullet }{a}_i(1-e_i^2)},$$

$${\mathbf{L}}_i=l_i\int f_i(\mathbf{n})\mathbf{n}d\Omega$$

$$\mathbf{L}=\sum _{i=1}^K{\mathbf{L}}_i.$$

$$E=-\frac{1}{2}\sum _{i,j=1}^K{J}_{ij}\iint f_i(\mathbf{n})f_j(\tilde{\mathbf{n}})g(\mathbf{n},\tilde{\mathbf{n}})d\Omega d\tilde{\Omega },$$

$$g(\mathbf{n},\tilde{\mathbf{n}})\equiv {\left(\mathbf{n}\cdot \tilde{\mathbf{n}}\right)}^2-\frac{1}{3}.$$

$$J_{ij}=\frac{3G}{8}\frac{m_i{m}_j\min (a_i,a_j{)}^2}{\max (a_i,a_j{)}^3};$$

$$g(\mathbf{n},\tilde{\mathbf{n}})=q_{\mu \nu }(\mathbf{n})q_{\mu \nu }(\tilde{\mathbf{n}})$$

$$q_{\mu \nu }(\mathbf{n})\equiv n_{\mu }{n}_{\nu }-\frac{1}{3}{\delta }_{\mu \nu }$$

$$Q_{\mu \nu }^{(i)}={⟨q_{\mu \nu }⟩}^{(i)}\equiv \frac{1}{N_i}\int f_i(\mathbf{n})q_{\mu \nu }(\mathbf{n})d\Omega$$

$${⟨\mathbf{X}⟩}^{(i)}=\frac{1}{N_i}\int f_i(\mathbf{n})\mathbf{X}d\Omega .$$

$${\epsilon }_i(\mathbf{n})=-q_{\mu \nu }(\mathbf{n})\sum _{j=1}^K{J}_{ij}{N}_j{Q}_{\mu \nu }^{(j)}=-n_{\mu }{V}_{\mu \nu }^{(i)}{n}_{\nu },$$

$$V_{\mu \nu }^{(i)}=\sum _{j=1}^K{J}_{ij}{N}_j{Q}_{\mu \nu }^{(j)}$$

$$E=\frac{1}{2}\sum _{i=1}^K{N}_i⟨{\epsilon }_i⟩=-\frac{1}{2}\sum _{i,j}^K{J}_{ij}{N}_i{N}_j{Q}_{\mu \nu }^{(i)}{Q}_{\mu \nu }^{(j)}.$$

$$S=-k\sum _{i=1}^K\int f_i(\mathbf{n})\ln f_i(\mathbf{n})d\Omega ,$$

$$\delta S/k+\sum _{i=1}^K{\alpha }_i\delta N_i-\beta \delta E+{\gamma }_{\mu }\delta L_{\mu }=0,$$

$$\mathbf{\gamma }=\beta \mathbf{\omega }$$

$$\delta S/k$$

$$\delta N_i$$

$$\delta E$$

$$\delta \mathbf{L}$$

$$\sum _{i=1}^K\int \delta f_i(1+\ln f_{\mathrm{eq},i}-{\alpha }_i+\beta {\epsilon }_i-l_i{\gamma }_{\mu }{n}_{\mu })d\Omega =0$$

$$f_{\mathrm{eq},i}=e^{-1+{\alpha }_i-\beta {\epsilon }_i+l_i{\gamma }_{\mu }{n}_{\mu }}.$$

$$f_{{\rm eq},i}(\bm{n})=N_{i}\frac{e^{-\beta\varepsilon_{i}(\bm{n})+{l}_{i}\gamma_{\mu}n_{\mu}}}{\phantom{\Big{|}}\int e^{-\beta\varepsilon_{i}(\bm{n})+{l}_{i}\gamma_{\mu}n_{\mu}}d\Omega}.$$

$$Q_{\mu\nu}^{(i)}=\frac{\int(n_{\mu}n_{\nu}-\frac{1}{3}\delta_{\mu\nu})e^{\beta n_{\rho}V^{(i)}_{\rho\sigma}n_{\sigma}+{l}_{i}\gamma_{\sigma}n_{\sigma}}d\Omega}{\phantom{\Big{|}}\int e^{\beta n_{\rho}V^{(i)}_{\rho\sigma}n_{\sigma}+{l}_{i}\gamma_{\sigma}n_{\sigma}}d\Omega},$$

$$Q_{\mu \nu }$$

$$\displaystyle=\left(\begin{array}[]{ccc}-\frac{1}{2}Q+\frac{1}{2}W&0&0\\ 0&-\frac{1}{2}Q-\frac{1}{2}W&0\\ 0&0&Q\end{array}\right)$$

$$Q$$

$$W$$

$$q$$

Full text

1. Introduction

The thermodynamics and statistical mechanics of a self-gravitating gaseous or stellar system is an intriguing subject with a long history (Tolman, 1934; Lynden-Bell, 1967; Lynden-Bell & Wood, 1968; Shu, 1978; Tremaine et al., 1986; Padmanabhan, 1990) and many recent developments (Gurzadian & Savvidy, 1986; Nakamura, 2000; Votyakov et al., 2002b; de Vega & Sánchez, 2002a, b; Arad & Lynden-Bell, 2005; Chavanis, 2006; Axenides et al., 2012; Touma & Tremaine, 2014; Bar-Or & Alexander, 2014; Roupas, 2015; Tremaine, 2015). Astrophysical applications include star clusters, galaxies, dark matter halos, and galaxy clusters. Statistical mechanics offers the hope of describing the macroscopic equilibrium structure of astrophysical systems without tracking the microscopic, particle-by-particle, evolution. The traditional approach, in which thermodynamic equilibrium is found by maximizing the entropy, does not work straightforwardly in self-gravitating systems for several reasons. First, there is no global entropy maximum for an isolated self-gravitating gaseous or stellar system (Antonov, 1962), although long-lived metastable thermal equilibrium states (i.e., local entropy maxima) can exist if the system is confined to a finite volume (Lynden-Bell & Wood, 1968; Chavanis, 2005). In addition, depending on the number of objects and the size of the system, the stochastic process by which the system evolves toward statistical equilibrium is often too slow for the system to approach the equilibrium state during its lifetime and therefore relaxation remains incomplete (Binney & Tremaine, 2008).

In the traditional thermodynamic approach, employed in most of the above references, the fundamental microscopic entities are stars, molecules, or elementary particles, considered to be point masses, which we label “bodies”. However, in contrast to systems with only short-range forces, the phase-space trajectories of bodies in self-gravitating systems are in many cases endowed with a special structure where the location of each body is restricted over long intervals to a bounded region of phase space, usually of lower dimensionality (Sridhar & Touma, 1999). This occurs if the dynamical system admits isolating integrals of motion (Kandrup, 1998; Binney & Tremaine, 2008) and allows the introduction of the familiar concept of the orbit. Rauch & Tremaine (1996) showed that certain degrees of freedom that describe the orbits may relax much faster than others, in a process they called resonant relaxation (see also Chavanis 2012; Fouvry et al. 2015). For example, bodies bound to a dominant central point mass execute eccentric orbits described by an angular-momentum vector and an eccentricity vector, which relax much faster than the semimajor axis. The phase-space distribution may then attain partial thermodynamic equilibrium, in which the rapidly evolving degrees of freedom are driven to internal thermodynamic equilibrium while others are frozen-in at their initial values.

Another example is a system dominated by a spherically symmetric potential, in which the directions of the angular-momentum vectors relax much faster than either the magnitudes of the angular-momentum vectors or the energies of the orbits (equivalently, much faster than either the eccentricities or semimajor axes). This process is called vector resonant relaxation (VRR) (Rauch & Tremaine, 1996; Hopman & Alexander, 2006; Gürkan & Hopman, 2007; Eilon et al., 2009; Kocsis & Tremaine, 2011, 2015), and is the subject of this work.

Gravitational systems that relax through VRR are conceptually similar to liquid crystals, as we will describe in more detail later in this paper. In normal liquids the microscopic degrees of freedom are the positions and momenta of the molecules, while their size may be ignored. However, in certain liquids, the molecules’ size and shape gives rise to new types of macroscopic structures. The liquid may acquire properties similar to a crystal: a macroscopic alignment of the molecules’ orientation with discrete rotational symmetries. For example in the nematic phase of a liquid crystal, axisymmetric molecules align in parallel or antiparallel configurations (Singh, 2002). Similarly, in gravitational systems subject to VRR, the basic microscopic entities, the orbits, behave as solid bodies. Their orientational degrees of freedom relax into an internal thermodynamic equilibrium. The correspondence with liquid crystals is due to the similarity between the Coulomb-type electro- and magnetostatic interaction acting between axisymmetric molecules and the Newtonian gravitational interaction acting between orbits. We shall show that there is a gravitational phase transition between disk+halo and isotropic phases111Note that this is different from the gravitational phase transitions commonly discussed in the literature, which usually refer to a transition between collapsed (core-halo structure) and diffuse phases in spherical systems (Aronson & Hansen, 1972; Stahl & Kiessling, 1994; Chavanis, 2006)., which strongly resembles the nematic-isotropic phase transition of liquid crystals.

Furthermore, we identify one more connection with condensed-matter physics, the existence of negative-temperature equilibria for which entropy decreases with increasing energy (Kocsis & Tremaine, 2011). Generally, negative-temperature equilibria may arise if the energy has an upper bound. In nuclear spin systems, the relaxation time of spin-spin interactions is typically much smaller than that of spin-lattice interactions (Ramsey, 1956), leading to a thermodynamic equilibrium for the former while the latter remain frozen at their initial values (Landau & Lifshitz, 1980). These systems admit negative-temperature configurations, in which most spins are anti-aligned with respect to the magnetic field. In VRR the role of spin is played by the orbital angular momentum and that of the magnetic field by the thermodynamic variable conjugate to the total angular momentum, the “rotation” (Eq. 20).

For a specific example of VRR, let us consider nuclear star clusters (NSCs), which consist of stars orbiting a supermassive black hole (SMBH). NSCs are found at the centers of galaxies and include the densest stellar systems in the Universe (Norris et al., 2014). The number density of various stellar types in NSCs increases steeply inwards with radius as $r^{-1.5}$ to $r^{-2.9}$ within the radius of influence of the SMBH, which is $\sim 1\,{\rm pc}$ for the NSC at the center of the Milky Way (Schödel et al., 2009, 2014; Bartko et al., 2009; Yelda et al., 2014; Chatzopoulos et al., 2015). In the Milky Way’s NSC, the low-mass old stars are spherically distributed. There is also a prominent population of young stars at radii between about 0.03 and $0.5\,{\rm pc}$.222Based on their colors and spectra these are O-type main-sequence stars, which have an estimated age of 4-6 Myr. The distribution of orbital planes of the young stars exhibits several structures with debated physical origin, including a disk in the inner region $0.03$–$0.12\,{\rm pc}$, a clumpy, disordered, roughly spherical structure in the intermediate region $0.12$–$0.2\,{\rm pc}$ and a warped, twisted disk in the outer region $0.2$–$0.5\,{\rm pc}$ (Bartko et al., 2009; Genzel et al., 2010). A fraction of stars in the intermediate region reside in a disk that rotates clockwise as seen from the Sun; a smaller fraction appear to reside in a disk that rotates counterclockwise, with orbit normals separated by $\sim 100^{\circ}$ from the normal to the clockwise disk; and the rest constitute a roughly spherical distribution (Levin & Beloborodov, 2003; Bartko et al., 2009). Yelda et al. (2014) estimate that $25\%$ of young stars are in the clockwise disk and $75\%$ are spherically distributed. The nearly identical ages of the young stars in both disks suggest a common origin, but why is the distribution of orbital planes so complex? Is it possible that this complicated structure can exist in statistical equilibrium, or is this a transient feature?

A further motivation to understand the relaxation processes in dense stellar populations comes from the emerging field of gravitational-wave astronomy. The recent LIGO discovery of gravitational waves from the merger of two black holes has opened a new window on the Universe (Abbott et al., 2016). Most of the massive stars formed in the past history of the Galaxy would by now have turned into neutron stars and stellar-mass black holes. The number density of stellar-mass black holes is expected to be up to a billion times higher in NSCs than in the Galactic field so NSCs may dominate the rate of black hole-black hole mergers detectable by LIGO (O’Leary et al., 2009). An understanding of the expected distribution of black holes in galactic nuclei may influence gravitational-wave search strategies and eventually help to interpret detections.

The main stellar-dynamical processes in NSCs may be ordered as follows. The spherical potential due to the central SMBH of mass $M_{\bullet}$ ($4\times 10^{6}M_{\odot}$ for the Milky Way) supports eccentric Keplerian orbits with a characteristic orbital time $t_{\rm orb}\sim(GM_{\bullet}/r^{3})^{-1/2}$, between $10^{2}$ and $10^{3}$ yr for distances (more precisely, semimajor axes) $r$ between $\sim 0.03$ and 0.5 pc. The gravitational field of the spherical distribution of old, low-mass stars causes in-plane apsidal precession for each eccentric orbit with a characteristic precession time $t_{\rm in-plane}\sim t_{\rm orb}M_{\bullet}/(Nm)$ or $10^{4}$–$10^{5}$ yr in this region, where $N$ is the enclosed number of stars and $m$ is their typical mass. The stellar orbits conserve their angular-momentum vector and orbital energy over timescales of this order. However on longer timescales the time-varying higher multipole moments of the gravitational field of the stars drive chaotic mixing in phase space, through three distinct processes that operate on different timescales (see Figure 1 in Kocsis & Tremaine 2011):

• •

$t_{\rm vrr}\sim t_{\rm orb}M_{\bullet}/(\sqrt{N}m)$ ($\sim 10^{6}$–$10^{7}$ yr), the diffusion time of angular-momentum vector directions or orbital planes, called the vector resonant relaxation time,

• •

$t_{\rm srr}\sim t_{\rm orb}M_{\bullet}/m$ ($\sim 10^{8}$–$10^{10}$ yr), the diffusion time of the magnitude of the angular momenta or eccentricity, called the scalar resonant relaxation time,

• •

$t_{\rm 2-body}\sim t_{\rm orb}M_{\bullet}^{2}/(Nm^{2})$ ($\sim 10^{9}$–$10^{10}$ yr), the diffusion time of energy or semimajor axis, called the two-body relaxation time.

The age of the young stars observed in the Galactic center is 4–6 Myr, long enough that we expect that VRR has strongly affected the distribution of orbital planes, but short enough that the eccentricities and semimajor axes remain frozen at their initial values.

More generally, in most of the volume of NSCs inside the sphere of influence of the central SMBH $M_{\bullet}\gg Nm$ and $N\gg 1$, and if these inequalities are satisfied we have the timescale hierarchy (Kocsis & Tremaine, 2011):

[TABLE]

This hierarchy implies that different dynamical properties of the system relax on different timescales and they may reach internal statistical equilibrium independently from one another (see Sridhar & Touma 2016a, b for a mathematically rigorous derivation of some of these results).

In this paper, we examine the equilibrium distribution of orbital planes after the VRR process is completed ($t\gtrsim t_{\rm vrr}$) but before scalar resonant relaxation begins ($t\lesssim t_{\rm srr}$). VRR is driven by the gravitational interaction averaged over two much faster motions, the orbital motion around the SMBH and the in-plane precession induced by the spherical distribution of old stars (Kocsis & Tremaine, 2015). Since precessing eccentric stellar orbits trace out axisymmetric punctured disks, the gravitational interaction is to be calculated between two such disks, in which the local surface density is proportional to the residence time of the star at the given location. Since $t_{\rm srr}$ and $t_{\rm 2-body}$ are much longer than the timescale of interest, the eccentricity and semimajor axis are nearly conserved during this process and so the surface density, as a function of radius, of the punctured disk is fixed, though its orientation is not333During VRR, the relative change in angular momentum is ${\Delta|L|}/{|L|}=(t_{\rm vrr}/t_{\rm srr})^{1/2}\sim N^{-1/4}$, and the change in energy is ${\Delta E}/{E}=(t_{\rm vrr}/t_{\rm 2-body})^{1/2}\sim(Nm/M)^{1/2}N^{-1/4}$. The former sets the change of $\sqrt{1-e^{2}}$ where $e$ is eccentricity and the latter sets the change in semimajor axis since $E=-GM_{\bullet}/(2a)$.. This leads to an averaged or effective Hamiltonian that determines the time evolution of the orientations of the orbit normals or angular-momentum vector directions (see Appendix A). The VRR Hamiltonian is given by $\frac{1}{2}N(N-1)$ terms corresponding to pairwise interactions between the punctured disks. Each pairwise interaction is further decomposed into a sum over multipoles. In this decomposition, the coupling constants depend on the conserved quantities of the two disks, mass, semimajor axis, and eccentricity, which are set by their initial distribution. In the terminology of statistical physics mass, semimajor axis, and eccentricity are quenched random variables, so the coupling constants in the Hamiltonian are random matrices as in spin glasses. If the external perturbations of the galaxy on the NSC may be neglected (as we shall assume in this paper), the energy corresponding to the VRR Hamiltonian and the total angular-momentum vector are also conserved.

Our main goal in this paper is to describe the mean-field theory of VRR for the dominant quadrupolar interaction, arguing that the qualitative features of the theory would be similar if higher order multipoles were included. The statistical equilibrium is described by a distribution function in orbit-normal space that extremizes the Boltzmann entropy. As we already noted the quadrupolar mean-field model of liquid crystals, the Maier-Saupe model (Maier & Saupe, 1958; Plischke & Bergersen, 2006), is analogous to the mean-field model of VRR for a one-component cluster in which the bodies have the same masses, semimajor axes, and eccentricities.

The model we investigate in this paper is also reminiscent of the Hamiltonian mean-field model (Campa et al., 2014), which describes the dynamics of $N$ identical rigid rotors interacting via a cosine potential. Both are exactly soluble and exhibit phase transitions and other interesting behavior. In both cases the simplicity arises because the interactions of all pairs of bodies are identical, and the dynamical behavior is entirely determined by a few moments of the distribution function.

The remainder of the paper is organized as follows. In Section 2 we formulate the mean-field theory of vector resonant relaxation. In Section 3, we restrict our attention to the special case of a cluster composed of bodies of a single mass, eccentricity and semimajor axis; this restriction simplifies the physics so that we can explore the thermodynamics and statistical mechanics analytically. We find all of the equilibrium solutions in the microcanonical and canonical ensembles; we present the axisymmetric equilibria in Section 5.1, and the non-axisymmetric ones in Section 5.2. In Section 6 we discuss a different ensemble, which we call the $\omega TN$-ensemble, in which the system is embedded in a bath with which it can exchange not only energy but also angular momentum. Section 7 contains a brief discussion of how these results can be applied to separable multi-component systems in which the interactions between components are weak. We discuss our conclusions in Section 8.

2. Mean-Field Theory of Vector Resonant Relaxation

We wish to describe the angular distribution of orbits of gravitating bodies bound to a spherically symmetric potential, usually dominated by a central point mass. To this end we adopt the VRR Hamiltonian obtained by Kocsis & Tremaine (2015) and reproduced in Appendix A. This Hamiltonian is derived by introducing a canonical transformation from Cartesian positions and momenta to action-angle variables, called Delaunay variables, and averaging over the rapidly varying angles, the mean anomaly and argument of periapsis (see also Sridhar & Touma 2016a, b, c). To leading order in the ratio of the stellar mass to the central mass, the VRR Hamiltonian is the average potential energy between the annuli or punctured axisymmetric disks that the precessing elliptical trajectories cover over times long compared to the apsidal precession time ($t_{\rm in-plane}$ in Section 1). The surface density of the annulus at any point is inversely proportional to the time the body resides at that position during its orbit. Since the semimajor axis and eccentricity are approximately conserved during VRR, these annuli conserve their intrinsic properties (i.e., surface-density distribution, inner and outer radii), but their orientations can vary due to VRR. Since we have averaged over 2 phase-space coordinates and an additional 2 phase-space coordinates are frozen-in (conserved), the relaxation is restricted to a 2-dimensional space, which can be taken to be determined by the direction of the unit vector parallel to the angular momentum of the orbit (i.e., the normal to the annulus), which we denote by $\bm{n}$. Alternative but less convenient coordinates are the inclination and the longitude of the node of the orbit. In summary, the basic structures of VRR are concentric orbit annuli with distinct masses, inner, and outer radii, which behave as rigid bodies pinned to the dominant central point mass, and precess due to their mutual gravitational torques.

An important property of the annulus, which is conserved during VRR, is its scalar angular momentum

[TABLE]

here we have denoted the mass of each body by $m$, its semimajor axis by $a$, its eccentricity by $e$, and the mass of the central object (e.g., the SMBH in a NSC) by $M_{\bullet}$.

We assume that the system is comprised of $K$ distinct groups of bodies. The $i^{\rm th}$ group contains $N_{i}$ bodies of identical mass $m_{i}$, semimajor axis $a_{i}$, and eccentricity $e_{i}$. We assume the ordering $a_{1}\leq a_{2}\leq\cdots\leq a_{n}$. This model can be thought of as representing an approximate coarse-grained distribution, assuming that the bodies are grouped into bins with similar mass, semimajor axis, and eccentricity. We shall refer to these groups as “components”. Since every orbit of the same component has the same scalar angular momentum $l_{i}$, we can visualize our construction in angular-momentum space as follows: the tips of the angular-momentum vectors of different bodies in the same component lie on a thin spherical shell, and since the scalar angular momentum of each body is conserved during VRR, the angular momenta can only move along this spherical shell, so the bodies in a given shell can interact with bodies in other shells but can only relax within their own shell.

We denote the distribution function or number per unit solid angle within shell $i$ by $f_{i}(\bm{n})$. We neglect two-body or higher-order correlation functions (the mean-field approximation).

Let us summarize the simplifying assumptions used in this paper.

The long-term evolution is driven by the mutual gravitational torques between concentric, axially symmetric, time-independent structures, the time-averaged orbits, which change in orientation in response to these torques. This assumption is generally valid for masses which execute circular orbits in an arbitrary spherical potential or general orbits in a smooth spherical potential that is not exactly quadratic or Keplerian. 2. 2.

The multipole expansion of the interaction potential between orbits is truncated at the quadrupole order. The quadrupole is generally the strongest multipole interaction for near-Keplerian orbits, whether radially overlapping or not, but the cumulative effects of higher order multipoles may be significant for radially overlapping orbits (Kocsis & Tremaine, 2015). 3. 3.

There are $K$ components, each comprised of $N_{i}\gg 1$ bodies with the same scalar angular momentum $l_{i}$, and the quadrupole coupling ($J_{ij}$ in Eq. 8) between any body in component $i$ and any body in component $j$ is the same. 4. 4.

Multi-body correlations are negligible, so the evolution can be described by a mean-field model. This assumption is expected to be valid if $N_{i}$ is sufficiently large. 5. 5.

The total number of bodies, the mass of each body, and their total vector angular momentum are conserved, except in Section 6, where the bodies are allowed to exchange angular momentum with a surrounding reservoir. 6. 6.

Either the total energy of the system (i.e., the total gravitational potential energy arising from pairwise interactions of orbits) is conserved or the VRR temperature is fixed (see Section 2.3). These assumptions correspond to the microcanonical and canonical ensembles, respectively.

After presenting the framework for a general multi-component model in this section, we will focus most of the paper on the simple case of a one-component system444In Section 7 we will discuss briefly how these results can be extended to multi-component systems in which the interactions between different components are much weaker than the interactions within a component., which exhibits a rich phase diagram and provides an analytically tractable starting point for comparisons with future numerical studies of multi-component systems.

2.1. Basic definitions

In the following, Greek indices label coordinates and Latin indices label stellar components. We use bold characters to denote 3-dimensional vectors, and normal (non-bold) symbols for their norm ($X=\|\bm{X}\|$). If we write out the Cartesian coordinates of some vector or tensor which has a stellar component index, we put the component index in parentheses in superscript to reduce clutter, e.g., $X^{(i)}_{\mu}$ denotes a Cartesian coordinate of the vector $\bm{X}_{i}$.

The system is characterized by the total number of bodies $N$, the total angular momentum $\bm{L}$, and the total orbit-averaged interaction energy $E$. In spherical coordinates, the number of bodies in the $i^{\rm th}$ component with $\bm{n}$ between $(\theta,\phi)$ and $(\theta+d\theta,\phi+d\phi)$ is $f_{i}(\bm{n})d\Omega$, where $d\Omega=\sin\theta\,d\theta d\phi$ is the solid angle, and therefore the number of bodies of the $i$th component is

[TABLE]

where unless otherwise noted the integral is over the full unit $2$-sphere. The total number of bodies is

[TABLE]

The angular momentum of one orbit in the $i^{\rm th}$ component is

[TABLE]

and the total angular-momentum vector of all bodies in the $i^{\rm th}$ component is

[TABLE]

The total angular momentum is

[TABLE]

The total interaction energy due to the VRR Hamiltonian up to the quadrupole order is given in Appendix A:

[TABLE]

with $d\tilde{\Omega}=\sin\tilde{\theta}d\tilde{\theta}d\tilde{\phi}$ and

[TABLE]

The Hamiltonian (8) resembles that of the Maier-Saupe model for liquid crystals (Maier & Saupe, 1958; Plischke & Bergersen, 2006). We assume that $J_{ij}$ are known constant parameters of the model, which are determined by the quadrupole potential energies of interacting time-averaged orbits. For example, for circular orbits

[TABLE]

see Appendix A for the more general case of eccentric orbits in a near-Keplerian potential. The array $J_{ij}$ is symmetric.

Eq. (9) may be rewritten as

[TABLE]

where here and in the following we suppress the summation over the Greek indices (Cartesian coordinates), and we have introduced the traceless quadrupole moment tensor of angular momenta

[TABLE]

where $\delta_{\mu\nu}$ is the Kronecker-$\delta$. Let us introduce the mean-field sum

[TABLE]

where we define the ensemble average of $\bm{X}$ over component $i$ as

[TABLE]

The orbit-averaged energy of a single body in the $i^{\rm th}$ component with unit normal $\bm{n}$ is, from Eq. (8),

[TABLE]

where we defined

[TABLE]

which is traceless. The total energy is

[TABLE]

The one-body energy (15) has the following properties.

It has inversion symmetry, $\varepsilon_{i}(\bm{n}_{i})=\varepsilon_{i}(-\bm{n}_{i})$. 2. 2.

Since $n_{\mu}n_{\nu}$ has eigenvalues $\{0,0,1\}$, the average of such matrices has eigenvalues between 0 and 1. Thus $Q^{(j)}_{\mu\nu}=\left\langle n_{\mu}n_{\nu}\right\rangle^{(j)}-\frac{1}{3}\delta_{\mu\nu}$ has eigenvalues between $\frac{2}{3}$ and $-\frac{1}{3}$; similarly, the superposition $-V^{(i)}_{\mu\nu}$ has eigenvalues between $-\frac{2}{3}\sum_{j}J_{ij}N_{j}$ and $\frac{1}{3}\sum_{j}J_{ij}N_{j}$. The energy of each body is confined between these bounds with local minima (maxima) along the eigenvector of $-V^{(i)}_{\mu\nu}$ associated with its smallest (largest) eigenvalue and along the opposite directions, and saddle points along the eigenvector associated with the intermediate eigenvalue and the opposite direction. 3. 3.

The equipotential curves with fixed $\varepsilon_{i}$ are generally ellipses on the unit sphere (i.e., the intersection of the unit sphere with an ellipsoid having the same center), which enclose the local minima and maxima. 4. 4.

The eigenvectors of $V_{\mu\nu}^{(i)}$ are fixed points of the time evolution of $\bm{n}$ for any body in component $i$, provided that $V_{\mu\nu}^{(i)}$ is constant. The eigenvectors corresponding to the smallest and largest eigenvalue are stable fixed points, and the intermediate eigenvector is unstable. 5. 5.

The time evolution of $\bm{n}$ is integrable if $V_{\mu\nu}^{(i)}$ is constant in time.

However the precession of the angular momenta typically changes $Q_{\mu\nu}^{(j)}$ and $V_{\mu\nu}^{(i)}$ in time which usually leads to chaotic evolution. We expect, and shall assume, that this chaotic evolution causes the distribution function to relax toward a state of maximum entropy.

2.2. Statistical equilibrium

We calculate the equilibrium distributions $f_{{\rm eq},i}$ that extremize the Boltzmann entropy (Jaynes, 1965)

[TABLE]

for fixed energy ${E}$, fixed total angular momentum $\bm{L}$, and fixed number of bodies in each component $N_{i}$. The use of this formula for the entropy in the mean-field approximation can be justified in the general case of the self-gravitating gas (Padmanabhan, 1990; Katz, 2003) (see also Miller 1973; Sormani & Bertin 2013 for critical analysis), and also for VRR equilibria as we show in Section 3.1.

We use the method of Lagrange multipliers. For perturbations $\delta f_{i}$ about the equilibrium distributions, the extremum must satisfy

[TABLE]

where $\alpha_{i}$, $\beta$, and $\gamma_{\mu}$ are Lagrange multipliers corresponding to the constraints and $\delta$ denotes the first-order variation with respect to $\{f_{1}(\bm{n}),f_{2}(\bm{n}),\dots,f_{K}(\bm{n})\}$. We may identify the Lagrange multiplier $\beta$ with the inverse temperature $\beta=1/(kT)$, while $\bm{\gamma}$ is the thermodynamic variable conjugate to $\bm{L}$ with respect to entropy. We emphasize that in this paper “temperature” $T$ refers to the VRR temperature, which is the inverse Lagrange multiplier enforcing conservation of VRR energy in Eq. (19) when maximizing the Boltzmann entropy; thus $\beta^{-1}=kT$. It is also the inverse derivative of the Boltzmann entropy with respect to VRR energy for a series of equilibria.

For the quantity enforcing the conservation of total angular momentum, it will later prove to be convenient to replace $\bm{\gamma}$ by $\bm{\omega}$, which we define as

[TABLE]

where $\bm{\omega}$ has units of angular velocity, and we refer to it as “rotation”555For an ideal gas, $\bm{\omega}$ describes rotation in the sense that the mean velocity of the gas at any position $\bm{r}$ is $\bm{\omega}\times\bm{r}$. However, in VRR $\omega$ is not simply related to the mean angular velocity. (see Section 6, and Votyakov et al. 2002b, a; de Martino et al. 2003). We will use this quantity to illustrate the analogy with paramagnetism and spin systems.

From the definition of $S$, $E$, $N_{i}$ and $\bm{L}$, we get to first order (for the second-order variation see Appendix C)

[TABLE]

where $f_{{\rm eq}}(\bm{n})$ is the equilibrium distribution function. Eq. (19) becomes

[TABLE]

Since the variations $\delta f_{i}$ are independent, the quantities in parentheses must vanish, which implies

[TABLE]

Using the constraint on $N_{i}$, Eq. (3), we may eliminate $\alpha_{i}$ to get finally

[TABLE]

where $\varepsilon_{i}$ depends on the orbit normal $\bm{n}$ through the tensor $q_{\mu\nu}$ as given by Eqs. (12), (13), and (15). In the following, we drop the subscript “eq” for brevity.

Eq. (13) gives the self-consistency equations

[TABLE]

where $V_{\mu\nu}^{(i)}$ is defined in Eq. (16). These self-consistency equations are subject to the constraints (7) and (17) for given $\bm{L}$ and $E$, which determine $\beta$ and $\bm{\gamma}$. These equations define the equilibria of the system and may be solved numerically, or in many cases analytically (see Appendix E).

2.3. On inequivalence of statistical ensembles

In statistical physics the term canonical ensemble is used to describe a system in equilibrium with a large heat reservoir having a fixed temperature. A microcanonical ensemble, on the other hand, describes an isolated system with constant total energy. The two ensembles are equivalent in the thermodynamic limit $N\rightarrow\infty$ for systems that are governed by so-called “short-range” interactions (Campa et al., 2014). The essential property of short-range systems is that they can be divided in an arbitrary manner into subsystems whose mutual interaction energy can be neglected with respect to the total energy of the system. Thus the subsystem energies are additive (Padmanabhan, 1990).

For systems with long-range interactions, such as the self-gravitating gas, the two ensembles may be inequivalent. The inequivalence leads to different stability properties for the two ensembles (see also Appendix C). A second condition must also be met for inequivalence to appear666More specifically, inequivalence may appear in regions where the entropy is not a concave function of energy for a series of equilibria.: the existence of a first-order phase transition in the canonical ensemble (Ellis et al., 2000; Barré et al., 2001; Ellis et al., 2004; Bouchet & Barré, 2005; Touchette & Ellis, 2005). The phase transition in the canonical ensemble is replaced by a stable region with negative specific heat in the microcanonical ensemble (Lynden-Bell & Wood, 1968; Padmanabhan, 1990; Campa et al., 2014). In the canonical ensemble the phase transition takes place in an out-of-equilibrium process, where the system absorbs or emits the whole latent heat from the heat bath needed for the transition. The two phases cannot coexist in equilibrium because the mixed phase has a higher free energy than a single pure phase in the canonical ensemble due to the interaction energy (Chavanis, 2006; Campa et al., 2014), and there cannot be a phase separation. In other words, there are negative specific heat equilibria in the phase transition region which are therefore unstable in the canonical ensemble777A system with negative specific heat cannot be in equilibrium with a heat bath. Any energy loss of the system to a colder heat bath will increase its temperature, making it lose even more energy, while any energy absorption of the system from a hotter heat bath will make the system colder causing a further energy absorption.. However, these equilibria are stable in the microcanonical ensemble, and the system can lie in any of those configurations between the two phases.

In the case of inequivalence, the equilibrium state depends on whether the system is under conditions of constant temperature (external environment has a strong influence which resembles a heat bath) in this state, in which case we refer to the canonical ensemble, or under conditions of constant energy (isolated system), in which case we refer to the microcanonical ensemble.

The self-gravitating gas is generally non-additive and presents inequivalence of ensembles (Padmanabhan, 1990; Campa et al., 2014). On the other hand, in VRR, the spatial and velocity distributions of bodies are not determined by the relaxation process, apart from the orientation of the orbits, leaving open the possibility for specific distributions to be approximately additive with respect to VRR energy. However, in more typical cases, VRR in a multi-component system is non-additive. For the one-component systems that are the focus of this paper, VRR is not only non-additive, but also presents a first-order phase transition and therefore inequivalence of ensembles. For this reason the microcanonical and canonical ensembles will be studied separately.

3. One-component systems

In the next three sections, we restrict ourselves to one-component systems; these are the simplest VRR systems yet they exhibit a remarkably rich phenomenology. The quadrupole mean-field VRR Hamiltonian of the one-component system is equivalent to the Maier-Saupe model of liquid crystals, which allows us to follow standard textbooks on the subject (Plischke & Bergersen, 2006)888The main difference is that we also require the conservation of total angular momentum.. First, we write down the one-component version of the relevant equations of Section 2, in particular the mean-field quadrupole moment, the total energy, the total angular momentum, the distribution function, and the self-consistency equation, and evaluate the entropy and the free energy. We solve the resulting system of equations in the following section.

We adopt a coordinate basis in which the axes are aligned with the eigenvectors of $Q_{\mu\nu}$. In the next subsection 3.1 and Appendix D we show that at equilibrium $\bm{L}$ and $\bm{\gamma}$ are parallel to one of these eigenvectors. Therefore, without loss of generality, we may choose the third coordinate axis to lie along $\bm{L}$. In this case we may write

[TABLE]

where we define

[TABLE]

Here $(\theta,\phi)$ are spherical coordinates in orthogonal axes aligned with the eigenvectors of $Q_{\mu\nu}$ and $\bm{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$. We have introduced the angular-momentum vector direction coordinates $(q,w)$, which are proportional to the $\ell=2$ real spherical harmonics999Specifically, $q=(16\pi/45)^{1/2}Y_{2,0}(\bm{n})$ and $w=(16\pi/15)^{1/2}\times Y_{2,2}(\bm{n})$, where $Y_{\ell,m}(\bm{n})$ are real spherical harmonics, which form an orthonormal basis on the sphere, $\int Y_{\ell,m}(\bm{n})Y_{\ell^{\prime},m^{\prime}}(\bm{n})d\Omega=\delta_{\ell,\ell^{\prime}}\delta_{m,m^{\prime}}$.. These coordinates are defined in the region

[TABLE]

From Eq. (15), the mean-field potential energy of a single body is simply

[TABLE]

All of our results are parameterized by $J$, so its value need not be specified for this analysis. For concreteness $J_{ij}$ is given for general near-Keplerian orbits in Appendix A, and for circular orbits Eq. (10) implies

[TABLE]

At equilibrium, the total energy from Eq. (17) is

[TABLE]

In Appendix B we show that the extrema of the entropy have the property that the Lagrange multiplier $\bm{\gamma}$ and the total angular momentum $\bm{L}$ are parallel. Without loss of generality we may label the direction of $\bm{L}$ as the positive $z$-axis. Letting

[TABLE]

we get $\bm{\gamma}\cdot\bm{n}=\gamma s$ and the total scalar angular momentum from Eq. (6) is

[TABLE]

Thus, $\left\langle s\right\rangle$ is the dimensionless total angular momentum of the system normalized to the configuration in which all angular momenta are aligned. This satisfies $0\leq L/(Nl)\leq 1$, where $L/(Nl)=1$ represents a razor-thin disk in physical space with all bodies orbiting counterclockwise as seen from the positive $z$-axis, and $L=0$ represents configurations with an equal net angular momentum in clockwise and counterclockwise orbits.

The parameters $Q$ and $W$ are bounded by the inequalities101010These follow from the bounds $0\leq\left\langle n_{\mu}^{2}\right\rangle\leq 1$ (since $n_{\mu}$ are Cartesian coordinates of a unit vector) and $0\leq\left\langle\cos\theta\right\rangle^{2}\leq\left\langle\cos^{2}\theta\right\rangle\leq 1$ (due to the Cauchy–Schwarz inequality).

[TABLE]

and the total energy is bounded by the inequalities

[TABLE]

The mean-field equilibrium distribution function (27) is written as

[TABLE]

In these equations the angular dependence of $f$ is implicit in $q$, $w$, and $s$. We now evaluate $Q$ and $W$ using their definitions in Eqs. (33) and (34). In so doing we use the identity $q=s^{2}-\frac{1}{3}$ and simplify the integrals over the azimuth angle $\phi$ using modified Bessel functions

[TABLE]

The self-consistency equations (33)–(34) become

[TABLE]

where the arguments of the Bessel functions are $\frac{1}{2}JN\beta W(1-s^{2})$. The moments $Q$ and $W$ are the order parameters of the system. Configurations with $W=0$ are axisymmetric around the $z$-axis while configurations with $W=\pm 3\,Q$ are axisymmetric around the $x$- and $y$-axes, respectively (recall that the positive $z$-axis points along the total angular-momentum vector). Therefore $W$ may be regarded as a measure of the deviation from axisymmetry around the angular-momentum axis. Configurations with $W=Q=0$ are isotropic.

The entropy of equilibrium states is calculated from Eqs. (18), (33), (34), and (45) to be

[TABLE]

From now on we drop the constant $N\ln N$ term. In this equation

[TABLE]

The argument of the Bessel function $I_{0}$ is given below Eq. (48). The total energy $E$ and angular momentum $\bm{L}=L\bm{e}_{z}$ in Eq. (49) are given in Eqs. (40) and (42), respectively.

In what we call the canonical ensemble, the system can exchange energy but not angular momentum with a reservoir. The corresponding thermodynamic potential is the Helmholtz free energy $F=E-TS$. Using Eq. (49) we get

[TABLE]

In this case, the parameter $\gamma=\gamma(T,L)$ appearing explicitly here and implicitly in $Z_{0}$ (see Eq. 50) is to be determined by Eqs. (42) and (45). In what we call the $\omega TN$-ensemble (see Eq. 20 and Section 6), the system can exchange both energy and angular momentum with a reservoir. Then the variables $\beta$ and $\gamma$ conjugate to $E$ and $L$ are conserved during the evolution of the system. The corresponding thermodynamic potential, analogous to the Gibbs free energy, is

[TABLE]

Typically we would like to solve the system of equations (47)–(48) and (42) for $Q$, $W$, and $\gamma$, given values of the angular momentum $L$, and the energy $E$ (for the microcanonical ensemble) or temperature $T$ (for the canonical ensemble). It is helpful to observe that $Q$, $W$, and $\gamma$ only appear in the combinations (see Appendix E)

[TABLE]

Thus the ratio of Eqs. (47) and (48) gives

[TABLE]

where the argument of the Bessel functions is $\kappa_{2}(1-s^{2})$. Moreover the angular momentum constraint (42) can be written as

[TABLE]

This last equation is a monotonic and therefore invertible function of $c\geq 0$ for given $\kappa_{1}$ and $\kappa_{2}$, mapping onto the interval $[0,1)$. Thus for any triple $(\kappa_{1},\kappa_{2},L)$, we can solve this numerically, and substitute the resulting value of $c(\kappa_{1},\kappa_{2},L)$ into Eq. (54) to get a relation between $\kappa_{1}$ and $\kappa_{2}$, although this relation may have multiple solutions for $\kappa_{2}$ at a given value of $\kappa_{1}$. For each $(\kappa_{1},\kappa_{2})$ pair we can evaluate the temperature, energy or other quantities using Eqs. (47), (53), and (40). Once we have determined the value(s) of $(\kappa_{1},\kappa_{2})$ corresponding to the desired temperature, the order parameters $Q$ and $W$ the entropy, and the free energies follow from Eqs. (53), (49), (51) and (52. In Appendix E we simplify the integrals analytically in Eqs. (47), (50), and (55) for axisymmetric states (see Eqs. E15–E18) and derive the asymptotic behavior of macroscopic variables.

3.1. Positive-temperature equilibria

In this Section we calculate the positive-temperature VRR mean-field free energy (51) as the steepest-descent approximation to the partition function and we describe the equilibrium states and their stability properties. The negative-temperature equilibria are discussed in Section 3.2 and Appendix C.

We work in the $\omega$TN-ensemble, that is we assume that the system is embedded in a heat bath with constant temperature $T$ and rotation $\omega$, with which it can exchange VRR-energy and angular momentum but not bodies. Therefore, in this ensemble $T$, $N$ and $\omega$ are held fixed (see also Section 6). The partition function is

[TABLE]

strictly the sum should not contain terms with $i=j$ but $\sum_{\mu\nu}q_{\mu\nu}^{(i)}q_{\mu\nu}^{(i)}=\frac{2}{3}$ for any $\Omega_{i}$ so these only contribute a constant factor to the partition function. Note that (i) the Latin indices run through different bodies in the single-component case and not through different components as in Section 2; (ii) for clarity, in this Section we explicitly write out sums over Cartesian coordinates indicated by Greek indices.

We will use the Hubbard–Stratonovich method to calculate this integral. We may write the first sum in Eq. (3.1) as

[TABLE]

Now let $Q_{\mu\nu}$ be an arbitrary $3\times 3$ matrix. For $T>0$ and $J>0$ we can write

[TABLE]

which can be shown by completing the square of $Q_{\mu\nu}$ on the right-hand side. Then, Eqs. (3.1) and (58) become

[TABLE]

Here we have used the fact that the integrals over $d\Omega_{1}\cdots d\Omega_{N}$ are separable and the integral is over the unit sphere, $n_{\mu}n_{\mu}=1$. We have also used the definition (12) and in the last line we introduced

[TABLE]

In this formula we identify $\{Q_{\mu\nu}\}$ as the order parameters.

We now specialize to the thermodynamic limit $N\rightarrow\infty$. We note that the integral in Eq. (60) is invariant as $N\to\infty$ for fixed $Q_{\mu\nu}$ as long as $\beta$ is rescaled as $\beta\propto N^{-1}$. With this rescaling $\beta g$ is independent of $N$. But since the integrand in (59) is $\exp(-N\beta g)$, this can be evaluated asymptotically for large $N$ by the method of steepest descent. We have

[TABLE]

The equilibrium free energy $G_{\rm eq}$ is given by $\ln\Xi=-\beta G_{\rm eq}$ so by dropping the fractional correction of order $(\ln N)/N$ we have

[TABLE]

Note that the extrema of $g(Q_{\mu\nu},\omega,T)$ occur when $\partial g/\partial Q_{\mu\nu}=0$, and these conditions give exactly the self-consistency equations (28), that is, the same distribution function as specified by use of the Boltzmann entropy (18).

In Appendix D we show that a coordinate system aligned with $\bm{\omega}$ can always be rotated in such a way that the off-diagonal elements of $Q_{\mu\nu}$ are zero at equilibrium and the system is stable with respect to perturbations in the off-diagonal elements in this coordinate system. Therefore we can finally write

[TABLE]

where

[TABLE]

which gives exactly Eq. (52).

More generally, a global minimum of free energy with respect to the order parameters is a stable equilibrium, while a local minimum is metastable. Saddle points and local/global maxima are unstable equilibria. A similar result holds for the free energy of the canonical ensemble $F$, given by

[TABLE]

where $G(Q,W,\omega,T)=Ng(Q,W,\omega,T)$ and now $\omega$ is determined implicitly by the requirement that the total angular momentum is $L$.

The top left panel of Figure 1 shows the free energy in the canonical ensemble as a function of the order parameters $Q$ and $W$ for $L/(Nl)=0.2$ at a low temperature, $kT/(JN)=0.05$. Points $\Sigma_{1}$–$\Sigma_{5}$ are the five distinct equilibrium configurations found at this temperature; the corresponding angular-momentum unit vector distribution functions $f(\bm{n})$ are shown below the panel. The dashed lines represent axisymmetric configurations. In particular,

• •

$\Sigma_{1}$, $\Sigma_{2}$ and $\Sigma_{3}$ are axisymmetric equilibria. $\Sigma_{4}$ and $\Sigma_{5}$ are non-axisymmetric (even though they lie close to the dashed lines).

• •

$\Sigma_{1}$ is stable since it is the global minimum of the free energy. We call this the “ordered phase” of the system or “uniaxial”, since it possesses one axis of symmetry (following the terminology of liquid crystals; Gramsbergen et al. 1986; Luckhurst & Sluckin 2015).

• •

$\Sigma_{2}$ is a saddle point, which is unstable with respect to perturbations that are not axisymmetric with respect to the $z$-axis.

• •

$\Sigma_{3}$ is a local maximum and hence an unstable equilibrium at this low temperature; however it is stable at higher temperature or at negative-temperature. It is nearly isotropic and we call it the “disordered phase”.

• •

$\Sigma_{4}$ is a saddle point and therefore an unstable equilibrium.

• •

$\Sigma_{5}$ is a local minimum and hence a metastable equilibrium. We call this state “biaxial”, since at low temperatures the bodies are organized into two distinct disks with different symmetry axes (Section 5.2). The state is long-lived for a significant range of temperatures and total angular momenta.

The top right panel of Figure 1 shows the canonical free energy at high temperature for the same mean angular momentum, $L/(Nl)=\left\langle s\right\rangle=0.2$. All equilibria merge at the point labeled $\Sigma$ which represents a stable disordered phase.

The geometry of the equilibrium states may be understood from the bottom panels in Figure 1, or from the eigenvalues of $\left\langle n_{\mu}n_{\nu}\right\rangle$ which from Eqs. (13) and (3)–(36) are $-\frac{1}{2}Q+\frac{1}{2}W+\frac{1}{3}$, $-\frac{1}{2}Q-\frac{1}{2}W+\frac{1}{3}$, and $Q+\frac{1}{3}$. For example, a razor-thin disk with all angular momenta pointing along the positive $z$-axis has distribution function $f(\bm{n})=N\delta(\bm{n}-\bm{e}_{z})$, where $\bm{e}_{z}=(0,0,1)$ is the unit vector along the $z$-axis, and so $\left\langle n_{\mu}n_{\nu}\right\rangle$ has eigenvalues $\{0,0,1\}$. In contrast, an isotropic distribution has $f(\bm{n})=N/(4\pi)$ and eigenvalues $\{\frac{1}{3},\frac{1}{3},\frac{1}{3}\}$. For the stable equilibrium $\Sigma_{1}$ in Figure 1 the eigenvalues are $\{0.03,0.03,0.93\}$. In physical space, this represents a relatively thin, axisymmetric disk containing both prograde and retrograde bodies, oriented perpendicular to the $z$-axis. The metastable equilibrium $\Sigma_{5}$ has eigenvalues $\{0.9,0.03,0.07\}$ which in physical space represents two disks with $157^{\circ}$ mutual inclination (see Eq. 71 below).

3.2. Negative-temperature equilibria

We remark that the internal energy given by Eq. (40) is bounded from above (Eq. 44). It is known that in this case, negative-temperature equilibrium configurations may exist for isolated systems (Ramsey, 1956; Pathria & Beale, 2012). Negative absolute temperature states were first introduced by Onsager (1949) in a statistical approach to turbulent flow aiming to explain long-lasting large vortices. This phenomenon may also arise in quantum systems with an upper energy bound; such a bound may arise for some degrees of freedom as long as they evolve in isolation from the motional degrees of freedom, as in the case of nuclear spins. Indeed, negative-temperature configurations were first produced in the laboratory in nuclear spin systems (Purcell & Pound, 1951; Ramsey & Pound, 1951; Ramsey, 1956). The spin-spin relaxation timescale is about six orders of magnitude smaller than that of the spin-lattice interactions in the system of Purcell & Pound (1951). A negative spin-temperature state was achieved for the spinorial degrees of freedom by cooling the system to very low temperature in the presence of a magnetic field and then rapidly reversing its orientation. The system maintained negative spin temperature for several minutes until relaxation of the motional degrees of freedom restored positive temperature (Ramsey, 1956).

Recently, a negative-temperature state was achieved in the laboratory for the motional degrees of freedom of ultra-cold quantum boson gases (Braun et al., 2013). This achievement triggered a debate (Schneider et al., 2014; Hänggi et al., 2015; Frenkel & Warren, 2015; Campisi, 2015; Cerino et al., 2015; Poulter, 2016), initiated by Dunkel & Hilbert (2014), on whether negative absolute thermodynamic temperatures are observable after all or maybe the Boltzmann’s entropy formula should be abandoned in favor of Gibbs’ entropy, which allows only positive temperatures. However, the arguments in favor of the validity of Boltzmann’s formula are more convincing at least for the case of ensemble-equivalence where it was proven recently that Boltzmann’s formula is appropriate and negative temperatures do occur (Buonsante et al., 2016).

Negative-temperature states in nuclear spin systems present a close analogy with the quadrupole VRR system discussed here. The VRR degrees of freedom, the unit normals $\bm{n}$, are similar to spinorial degrees of freedom. Orbits averaged over a precession period represent fixed density annuli which are characterized by their normals and relax similarly to spin vectors in quantum mechanics; the difference is that the interaction is a sum of terms proportional to $(\bm{n}_{i}\cdot\bm{n}_{j})^{2}$ for VRR and $\bm{s}_{i}\cdot\bm{s}_{j}$ for spin systems. The analogy goes even further since $\bm{\omega}$, defined in Eq. (20), enters the distribution function through $\bm{\omega}\cdot\langle\bm{n}\rangle$ similar to a paramagnetic term, $\bm{B}\cdot\left\langle\bm{s}\right\rangle$.

Negative-temperature equilibria also arise for an isolated thin circular disk of bodies undergoing VRR (Kocsis & Tremaine, 2011). That model is recovered in the limit $L/(Nl)\approx 1$ of the model presented here. With this perspective, negative-temperature equilibria arise naturally for VRR.

At negative temperature, we find a single solution to the self-consistency equation (47)–(48) at fixed angular momentum $\bm{L}$. In Appendix C, we show that these negative-temperature equilibria are stable in both the canonical and microcanonical ensembles. These states are always axisymmetric (Appendix E.2), and they are similar to the large positive-temperature disordered state shown in the top right panel of Figure 1 although negative-temperature states always exhibit a population inversion: the occupation number is an increasing function of one-body VRR-energy (Eq. 38).

Before closing this section, let us introduce some quantities for further use. We denote by $N_{+}$ the number of orbits with angular-momentum vectors on the $s\geq 0$ hemisphere (i.e. prograde orbits with respect to $\bm{L}$) and by $N_{-}$ the number of orbits with $s\leq 0$ (retrograde orbits with respect to $\bm{L}$)

[TABLE]

We also introduce the expectation value of a quantity $X$ over bodies on the $s\geq 0$ hemisphere as

[TABLE]

and similarly $\left\langle X\right\rangle_{-}$ for the $s<0$ hemisphere. In particular, we shall use $L_{\pm}$ to denote the mean angular momentum of the prograde and retrograde orbits with respect to the $z$-axis.

4. Zero angular momentum

It is instructive to start with the configurations having zero total angular momentum, $L=\gamma=0$ in Eqs. (42)–(45). In the one-component case, this is completely equivalent to the Maier-Saupe model of liquid crystals.

Figure 2 shows the free energy $F/(NkT)$ as a function of the order parameters $Q$ and $W$ at the same temperature $kT/(JN)=0.05$ as in the left panel of Figure 1. We may observe a similar structure as in that Figure with correspondence $\Sigma_{1}\rightarrow S_{1}$, $\Sigma_{2}\rightarrow S_{2}$, $\Sigma_{3}\rightarrow S_{3}$, $\Sigma_{4}\rightarrow S^{\prime\prime}_{2}$, $\Sigma_{5}\rightarrow S^{\prime}_{1}$. The main difference is that all equilibria are exactly axisymmetric ($W=0$ or $Q=\pm\frac{1}{3}\,W$, see Eq. 3). The primed and double primed states may be obtained by a $90^{\circ}$ rotation of the unprimed distribution.

The left panel of Figure 3 shows the order parameter $Q$ for a series of axisymmetric equilibria as a function of temperature (only positive temperatures are shown). There are three zero-temperature states $O_{1}$, $O_{2}$ and $O_{3}$ with $Q=\frac{2}{3}$, $-\frac{1}{3}$, and 0. The state $O_{3}$ corresponds to the isotropic distribution function

[TABLE]

which satisfies the self-consistency conditions (42), (48), and (47) when $Q=W=\gamma=0$ for any inverse temperature $\beta$; thus $O_{3}$ is identical to $S_{3}$. The state $O_{1}$ has $Q=\frac{2}{3}$, and since this is the maximum allowed value of $q$ by Eq. (37) all bodies in the system must have $q=\frac{2}{3}$. In other words all of the orbits have unit normals $\bm{n}=(0,0,\pm 1)$ so they form a razor-thin disk in the $x$-$y$ plane with equal numbers of prograde and retrograde orbits and energy $E=-\frac{1}{3}JN^{2}$. This is the most extreme configuration of what we call the ordered phase. The state $O_{2}$ has $Q=-\frac{1}{3}$ and since this is the minimum allowed value of $q$ all bodies must have $q=-\frac{1}{3}$, that is, their orbit normals lie in the equatorial plane $\theta=\frac{1}{2}\pi$. The azimuths of the orbit normals are uniformly distributed. These three states are respectively analogous to $S_{1}$, $S_{2}$, and $S_{3}$ in Figure 2 for zero temperature, which indicates that $O_{1}$ is stable, while $O_{2}$ and $O_{3}$ are unstable. A mathematical description of these states is given in Appendix E. The right panel of Figure 3 shows the mean angular momenta of the clockwise and counterclockwise orbits respectively, $L_{\pm}=Nl\left\langle s\right\rangle_{\pm}$ as defined in Eq. (67). As we have already shown, the stable zero-temperature state $O_{1}$ corresponds to a razor-thin disk with equal numbers of bodies orbiting clockwise ($s=-1$) and counterclockwise ($s=1$). At higher temperatures most of the bodies remain in a disk, which thickens as the temperature increases. In addition a few bodies are found at higher inclinations forming a “halo” as shown in Figure 4. As Eq. (68) implies, the disordered phase corresponding to $O_{3}$ is isotropic at all temperatures for zero total angular momentum ($Q=W=\left\langle s\right\rangle=0$).

Figure 3 marks the states $P_{\rm ord}$, $A$, $B$, and $P_{\rm dis}$, where stability properties change in the canonical ensemble, as shown in Figures 5 and 6 and in Appendix E.1111111For the case of zero total angular momentum discussed here, $B$, $P_{\rm dis}$, $D$, and $O_{3}$ are identical: they all represent the isotropic configuration, $f(\bm{n})=N/(4\pi)$. However, these states are distinct for nonzero angular momentum $L$. The isotropic disordered state is unstable at low temperature, but in the canonical ensemble it becomes metastable for temperatures in the range $T_{B}<T<T_{P}$ where $kT_{B}/(JN)=\frac{2}{15}=0.133333$ and $kT_{P}/(JN)=0.146796$. These states become stable for $T>T_{P}$. The ordered phase (branch $O_{1}P_{\rm ord}$) becomes metastable in the canonical ensemble for temperatures in the range $T_{P}<T<T_{A}$, where $kT_{A}/(JN)=0.148556$. At $T_{P}$, the system suffers a first-order phase transition from the ordered phase $P_{\rm ord}$ to the disordered phase $P_{\rm dis}$. At the transition point $P_{\rm ord}$, $Q=0.286014$, $L_{\pm}/(Nl)=\pm 0.747913$ and $E/(JN^{2})=0.0613553$, while at $P_{\rm dis}$, $Q=0$, $L_{\pm}/(Nl)=\pm\frac{1}{2}$, and $E/(JN^{2})=-\frac{1}{3}$. The distribution functions at $P_{\rm ord}$ and $P_{\rm dis}$ are shown in Figure 4.

To understand the nature of the phase transition, let us imagine slowly heating a one-component system from low temperature across the phase transition in the canonical ensemble121212This would occur, for example, if the system interacts with massive distant objects with a nearly isotropic angular-momentum distribution.. At low temperatures, the equilibria along the sequence starting with $O_{1}$ are minima of the free energy and therefore stable, while those along the sequences starting with $O_{2}$ and $O_{3}$ are saddle points and maxima, respectively, and therefore unstable. A comparison of Figures 2, 3, and 5 shows how the equilibria along these sequences change as the temperature increases: the sequences starting at $O_{1}$ and $O_{2}$ move toward the isotropic configuration, in the sense that the disk thickens and $|Q|$ decreases for both. At $kT_{B}/(JN)=2/15=0.133333$, the unstable sequences starting from $O_{2}$ and $O_{3}$ intersect. As the temperature increases across $T_{B}$, the equilibria on the isotropic sequence that started at $O_{3}$ change from maxima to minima of the free energy so the sequence becomes metastable, while the sequence starting from $O_{2}$ remains unstable. The stable ordered sequence that begins at $O_{1}$ remains stable until the temperature increases to $T_{P}=0.146796\,JN/k$. At $T_{P}$, the free energies of the ordered and disordered states, $P_{\rm ord}$ and $P_{\rm dis}$ in Figure 3, become equal. When the temperature is increased above $T_{P}$, a first-order phase transition occurs in the canonical ensemble, in which the isotropic disordered state becomes stable and the ordered state becomes metastable. As the temperature continues to increase, the sequences that began at $O_{1}$ and $O_{2}$ approach the same state, and they eventually coincide at $T_{A}=0.148556\,JN/k$, For $T>T_{A}$ there is no equilibrium other than the isotropic disordered state, which is stable.

The behavior near the phase transition is explored further in Figure 5. The left panel shows the free energy as a function of the order parameters $Q$ and $W$ at the temperature of the phase transition, $T_{P}$. Since the equilibria are all axisymmetric, we can set $W=0$ and plot the free energy as a function of one order parameter $Q$, which we do in the right panel. The free energy is shown for three temperatures. The bottom curve is for temperatures between $T_{B}$ and $T_{P}$, where the minimum at $P_{\rm ord}$ corresponds to the stable ordered sequence starting at $O_{1}$, the minimum at $P_{\rm dis}$ corresponds to the metastable disordered sequence starting at $O_{3}$, and the intervening maximum corresponds to the unstable sequence starting at $O_{2}$. The middle curve is at the temperature $T_{P}$ of the phase transition, where the free energies at $P_{\rm ord}$ and $P_{\rm dis}$ are equal, and the upper curve is for temperatures between $T_{P}$ and $T_{A}$.

Figure 6 shows $F/(NkT)$ near the phase transition. At $T_{P}$, the free energy is the same for the ordered and disordered phases ($P$ in the left panel), but the energies of these states are different ($P_{\rm ord}$ and $P_{\rm dis}$ in the right panel). This panel shows the caloric curve and Maxwell’s construction for the first-order phase transition. In a nonadditive system multiple phases cannot coexist in equilibrium (see Section 2.3) and so the dotted line $P_{\rm ord}P_{\rm dis}$ in the right panel of Figure 6 is unphysical. However, phase separation can occur in separable multi-component systems, as described in Section 7; in this case individual components are either in $P_{\rm ord}$ or in $P_{\rm dis}$ but the full system can lie anywhere along $P_{\rm ord}P_{\rm dis}$.

This analysis is different in the microcanonical ensemble for a one-component system, i.e., an isolated system under conditions of constant energy. In that case, we find that the series of equilibria $P_{\rm ord}ABP_{\rm dis}$ represent the highest entropy states at fixed energy, while the states along $O_{2}B$ and $O_{3}B$ at fixed energy have an entropy minimum. Thus, the branch $P_{\rm ord}ABP_{\rm dis}$ is stable in the microcanonical ensemble, and a continuous second-order phase transition takes place at point $B$ in the microcanonical ensemble at temperature $kT_{B}/(JN)=2/15=0.133333$ (see Figure 6). Branches $BO_{2}$ and $BO_{3}$ are unstable in both canonical and microcanonical ensembles.

5. Non-zero angular momentum

Let us turn to the more general case of a one-component system with non-zero total angular momentum, $L\neq 0$ in Eq. (42). In this case, the angular momentum constraint can only be satisfied if the Lagrange multiplier $\bm{\gamma}=\beta\bm{\omega}$ is non-zero, which introduces a factor $\exp(l\beta\bm{\omega}\cdot\bm{n})$ in the distribution function (Eq. 45). This is similar to the effect of a paramagnetic term in the Hamiltonian, and hence the distribution function, of a spin system in a magnetic field, where the role of magnetic field is played by $\bm{\omega}$, the magnetic moment is replaced by $l\bm{n}$ and the spin is replaced by $\bm{n}$.131313This is different in the Maier-Saupe model of liquid crystals, where an external magnetic field $\bm{B}$ gives rise to a term proportional to $(\bm{B}\cdot\bm{n})^{2}$ where $\bm{n}$ describes the orientation of the molecules (Wojtowicz & Sheng, 1974; Palffy-Muhoray & Dunmur, 1982)..

We work in a coordinate system in which the total angular momentum is parallel to the positive $z$-axis. The equilibrium distribution function $f(\bm{n})$ may be either axisymmetric around the $z$-axis (e.g., $\Sigma_{1}$, $\Sigma_{2}$, $\Sigma_{3}$ in Figure 1) or else non-axisymmetric, with three distinct eigenvalues (i.e., $Q\neq\pm\frac{1}{3}W$ and $W\neq 0$ as in $\Sigma_{4}$, $\Sigma_{5}$ in Figure 1). We study these cases separately.

5.1. Axially symmetric equilibria

Let us consider first axisymmetric configurations (e.g., $\Sigma_{1}$, $\Sigma_{2}$, $\Sigma_{3}$ in Figure 1). These constitute the most important cases, since global free energy minima at fixed temperature and global entropy maxima at fixed total energy are always axisymmetric as we will see below. We derive analytic expressions for the order parameters valid at any temperature and angular momenta in Appendix E.2.

Figure 7 shows the order parameter $Q$ with respect to temperature for axisymmetric equilibria, at several values of the dimensionless total angular momentum. For small values of $L$, three axisymmetric ground states exist at zero temperature, similar to the case with $L=0$ shown in Figure 3. The stable ordered state $O_{1}$ has $Q_{O_{1}}=\frac{2}{3}$ at $T\rightarrow 0^{+}$, for any $L$. This state is a razor-thin disk in physical space, with a fraction of prograde orbits equal to $\frac{1}{2}(1+L/(Nl))$. As $L$ increases, $Q_{O_{2}}$ at $T=0^{+}$ moves toward $Q_{O_{3}}=0$ and they merge at $L/(Nl)=1/\sqrt{3}$. This behavior is explained using the asymptotics of the partition function in Appendix E. Figure 8 shows the mean angular momentum $L_{\pm}$ and the ratio of prograde or retrograde orbits $N_{+}/N_{-}$ at the critical total angular momentum, $L_{\rm cr}/(Nl)=0.13714$. This is the maximum total angular momentum for which a phase transition occurs.

The first-order phase transition in the canonical ensemble presented in Figure 5 for $L=0$ is similar for any $0<L<L_{\rm cr}$. As $L$ is increased, the jumps $Q_{P_{\rm ord}}-Q_{P_{\rm dis}}$ and $E_{P_{\rm dis}}-E_{P_{\rm ord}}$ shrink, and above $L_{\rm cr}$ there is no phase transition. This is illustrated in Figure 10, which shows the region in energy-angular momentum space that corresponds to the phase transition.

On the critical curve $L=L_{\rm cr}$, the phase transition in the canonical ensemble becomes second-order. Figure 9 shows the order parameter with respect to temperature for several values of angular momentum below and above the critical value, and the free energy and the caloric curve are shown in Figure 11. The critical point is labeled with $C$ on the $L_{\rm cr}$ critical curve. In the left panel of Figure 9 we also plot the curve defined by the phase transition points (dashed line), the “coexistence curve”141414Although in the one-component model there cannot be phase coexistence in equilibrium, see Section 2.3., which is found empirically to follow a parabolic rule

[TABLE]

where $T_{P}$ is the temperature of the corresponding first-order phase transition and $T_{C}$ the temperature at the critical point. The right panel of Figure 9 shows the temperature-angular momentum phase diagram. In the vicinity of the critical point and on the critical curve, where the phase transition becomes second-order, we get the scaling

[TABLE]

These scalings are the ones predicted by mean-field theory (Stanley, 1987; Papon, 2002) and the critical point as described above is completely analogous to the critical points displayed not only by liquid crystals, but also by ferromagnetic and liquid-vapor systems (Stanley, 1987; Wojtowicz & Sheng, 1974; Palffy-Muhoray & Dunmur, 1982). The thermodynamic quantities $Q-Q_{C}$, $T$, $L$ of the one-component quadrupole VRR system discussed here correspond to $S-S_{C}$, $T$ and $H^{2}$, respectively, for liquid crystals in a magnetic field (see Palffy-Muhoray & Dunmur 1982 for the definitions), or $M$, $H$ and $T$ for ferromagnets, or $(\rho-\rho_{C})$, $P$ and $T$ for liquid-vapour systems (see Stanley 1987).

The first-order phase transition shown for $L\neq 0$ in Figure 11 is qualitatively similar to that for $L=0$ shown in Figure 6. Figure 12 (left panel, black curve) shows that the disordered phase is nearly isotropic (for all $L<L_{C}$ for which the phase transition is defined), while the ordered phase corresponds again to a disk+halo structure in physical space (cf. Figure 4 and see Section 3.1). In the ordered phase the system consists of a disk with two components, rotating in opposite directions, and a dilute halo. Due to our definition of the coordinate system, more bodies have $L_{z}>0$ than $L_{z}<0$ as shown by the asymmetry of the curve in the left panel of Figure 12. For high values of $L$ the system is a thin disk rotating in practically one direction only as depicted in Figure 13. Figure 14 shows the distribution function for states along the unstable branches.

For non-zero angular momentum, there is no phase transition in the microcanonical ensemble, but the system passes continuously from the ordered to the disordered phase as the total VRR energy is slowly increased. The caloric curve is shown in Figure 15 for the microcanonical ensemble, while the entropy with respect to energy is presented in Figure 16. The equilibrium configurations are identical in the microcanonical ensemble to those in the canonical ensemble apart from the phase-transition region of the canonical ensemble. In both ensembles, the system presents a disk+halo structure at low temperature and a nearly isotropic structure at high temperature for low total angular momentum. However, there are no stable configurations in the phase-transition region in the canonical ensemble between $P_{\rm ord}$ and $P_{\rm dis}$ and the system passes discontinuously from the one phase to the other, while in the microcanonical ensemble the corresponding configurations (branch $P_{\rm ord}ABP_{\rm dis}$ of Figure 11) are stable.

5.2. Non-axisymmetric equilibria

In Figure 1, we have seen that VRR free energy extrema exist with $W\neq 0$ and $W\not=\pm 3Q$ ($\Sigma_{4}$, $\Sigma_{5}$, and their rotations $\Sigma^{\prime}_{4}$, $\Sigma^{\prime}_{5}$). For $L\not=0$, these states represent non-axisymmetric configurations in the sense that the tensor $\left\langle n_{\mu}n_{\nu}\right\rangle$ has three distinct eigenvalues (see Eq. 3).

Figure 17 shows the order parameters of the non-axisymmetric equilibria as a function of temperature, for four values of the total angular momentum. For $L/(Nl)\leq 1/\sqrt{2}$ and sufficiently low temperature, there are two equilibria shown by solid and dashed lines, which correspond respectively to $\Sigma_{5}$ and $\Sigma_{4}$ in Figure 1. The solid lines represent local free-energy minima which are therefore metastable, and the dashed lines represent free-energy saddle points which are unstable (see Figure 1).

Figure 18 shows the two-dimensional probability density distribution of the orbit normals for the non-axisymmetric equilibria. When the angular-momentum vectors are concentrated near a single direction, as in the left panels, the distribution resembles a disk in physical space in which each object orbits in the same sense. Thus, in physical space, the metastable non-axisymmetric equilibria in the left panels are comprised of two misaligned but otherwise identical disks with axes along $(\cos\theta,\phi)=(\left\langle s\right\rangle,0)$ and $(\left\langle s\right\rangle,\pi)$; half of the bodies are in each disk. The mutual inclination of the disks is approximately

[TABLE]

We call these “biaxial” equilibria.

Let us now perform a thought experiment in which we heat the system, and visualize how the distribution of angular momenta changes. At $T=0$, the metastable state is $O_{5}$ (analogous to $\Sigma_{5}$ in Figure 1), consisting of two razor-thin disks in physical space with angular momenta pointing along directions $(\theta,\phi)=(\cos^{-1}(L/(Nl)),\pm\pi)$151515This configuration is metastable up to a global rotation of the coordinate system since the maximum entropy equilibrium state does not specify the direction of the eigenvectors of $\bm{Q}$ perpendicular to $\bm{\omega}$. Integrating Hamilton’s equations of motion for the VRR system shows that the two disks precess uniformly around the total angular-momentum vector (Kocsis & Tremaine, 2015).. As we heat the system, we increase the scatter of the angular-momentum vectors, or in physical space, the thickness of the disks. As the temperature continues to increase, the scatter of the angular-momentum vectors eventually becomes so large that the two angular-momentum direction distributions become connected, forming an inhomogeneous thick arc (or possibly a complete ring for small $\left\langle s\right\rangle$) in the $x$–$z$ plane with two density maxima. Eventually, above some maximum temperature $T_{I}$, only axisymmetric states exist. A second non-axisymmetric sequence begins at $T=0$ with a razor-thin arc along the $x$–$z$ plane, which is the unstable $O_{4}$ ground state analogous to $\Sigma_{4}$ in Figure 1. As the temperature increases the arc thickens. At point $G$ in Figure 17, the distribution function becomes bimodal and continues to thicken, until it finally connects with the $O_{5}$ sequence at $T_{I}$.

Figure 17 shows that for any fixed total angular momentum, non-axisymmetric equilibria only exist below a maximum temperature $T_{I}$ corresponding to point $I$ in the figure. For $T<T_{I}$ both metastable and unstable configurations are allowed in the canonical ensemble. $T_{I}$ decreases with increasing $L$ or $\left\langle s\right\rangle$. These equilibria are similar to161616$O_{1}$ and $O_{5}$ are equivalent in the limit $L\rightarrow 0$, and similarly for $O_{2}$ and $O_{4}$. the $O_{1}$ and $O_{2}$ branches of zero angular momentum axisymmetric equilibria for $T<T_{A}$ shown in Figure 3. At this maximum temperature an instability171717 The term “gravitational zeroth order phase transition” has been applied to analogous instabilities of the self-gravitating gas in de Vega & Sánchez (2002a) where it is used to describe the collapse of classical Newtonian self-gravitating gas and has been used in Chavanis (2002) to describe the condensation of a fermionic Newtonian self-gravitating gas.

arises in the canonical ensemble, which is also evident from the caloric curve shown in the left panel of Figure 19, and the system’s transition to the $W=0$ axisymmetric phase. Indeed, the left panel of Figure 19 shows that the stable axisymmetric equilibrium has lower free energy than the metastable non-axisymmetric equilibrium of the same temperature and angular momentum. Similarly, Figure 20 shows that axisymmetric equilibria have higher entropy than non-axisymmetric equilibria of the same energy and angular momentum.

The left panel of Figure 21 displays the region in the VRR temperature and total angular momentum space where biaxial equilibria exist. The maximum biaxial temperature $T_{I}$ is a decreasing function of $\left\langle s\right\rangle\equiv L/(Nl)$, and there are no biaxial states above $\left\langle s\right\rangle=\left\langle s\right\rangle_{\rm cr}\sim 0.7$. The value of $\left\langle s\right\rangle_{\rm cr}$ can be derived analytically by using the distribution function in angular-momentum direction space shown in Figure 18, and exploiting the fact that $T_{I}=0$ at $\left\langle s\right\rangle_{\rm cr}$. At $T=0$, we may assume that the bodies are arranged in $K$ razor-thin disks with normals $(\theta_{i},\phi_{i})$ and $s_{i}\equiv\cos\theta_{i}$:

[TABLE]

In particular, the biaxial $O_{5}$ state has $K=2$ with

[TABLE]

To determine whether this state is stable or unstable to splitting each peak of the distribution function to two peaks centered at slightly smaller and larger $s_{i}$, we examine if the energy of those configurations is lower or higher than that of $O_{5}$ using Eqs. (72–76). At fixed $N$, $L$, and $T=0$, we find that these perturbations result in an energy decrease exactly if $L/(Nl)\geq 1/\surd{2}$ corresponding to $\left\langle s\right\rangle_{\rm cr}=1/\surd{2}$. Hence the mutual inclination between the two axes in the zero-temperature biaxial state must be larger than $90^{\circ}$ for metastability. We conjecture that $O_{5}$ is metastable at higher inclinations and smaller $\left\langle s\right\rangle$ against arbitrary perturbations of the distribution function181818We find stability against a variety of perturbations including more general bifurcations into four modes with $N_{i}\rightarrow N_{i}+\delta N_{i}$, $\bm{n}_{i}=\bm{n}_{i}+\bm{\delta n}_{i}$ for $i\leq 4$ where $\sum_{i}\delta N_{i}=0$ and for the additional modes $N_{3}=N_{4}=0$ with $\delta\bm{n}_{i}$ arbitrary..

A rough estimate of the lifetime of the biaxial metastable states is obtained by assuming that the transition probability for each body in a unit VRR relaxation time $t_{\rm vrr}$ (see Section 1) is of order $\exp[\Delta F/(NkT)]$, where $\Delta F$ is the free-energy barrier between the metastable biaxial state $\Sigma_{5}$ and the stable axisymmetric state $\Sigma_{1}$. To find $\Delta F$, note that the transition from $\Sigma_{5}$ to $\Sigma_{1}$ requiring the least free energy passes through $\Sigma_{4}$, the unstable non-axisymmetric state (see Figure 1). In the most conservative estimate (Chavanis, 2005), the $\Sigma_{5}\rightarrow\Sigma_{4}\rightarrow\Sigma_{1}$ transition occurs when all $N$ objects climb the barrier coincidentally during their thermal motion, which happens in a time

[TABLE]

The right panel of Figure 21 shows $\Delta F/(NkT)$. Based on this simple estimate, the lifetime of the biaxial metastable states increases exponentially with $N$ (recall that at fixed values of the order parameters and coupling constant $F\sim N^{2}$ and $T\sim N$), and thus the metastable state is very long-lived for $N\gg 1$. In the opposite extreme, the transition to the $\Sigma_{1}$ state could occur gradually with different objects jumping over the barrier in succession. This happens over time $\exp[(\Delta F/(NkT)]t_{\rm vrr}$, which is much shorter. Numerical simulations are needed to determine the correct scaling with $N$.

6. $\omega TN$-ensemble

In Section 3.1 we introduced the $\omega TN$-ensemble. The system is assumed to be in a heat bath, a much bigger system that surrounds our small one-component system, with which it can exchange VRR-energy and angular momentum but not bodies. Thus in this ensemble the temperature $T$, the number of bodies $N$, and the rotation parameter $\bm{\omega}$ (Eq. 20) conjugate to angular momentum are held fixed. The generalized thermodynamic potential of this ensemble is given in equation (52), $G(\bm{\omega},T,N)=E-\bm{\omega}\cdot\bm{L}-TS.$ We focus our attention to the axisymmetric case here. We derive a parametric solution to the self-consistency equations in Appendix E.2. Similar to the canonical ensemble, we find a phase transition, but the values of the order parameters at the transition points are modified as shown in Figure 22 (cf. Figures 9 and 10). A first-order phase transition occurs for $\omega<\omega_{C}$; the transition becomes second-order for $\omega_{C}=0.057285JN/l$. This behavior resembles the canonical ensemble but the critical values are different with $T_{C}=0.15394JN/k$, $L_{C}=0.17173Nl$ and $Q_{C}=0.14273$ as depicted in Figures 22 and 23.

7. Separable multi-component systems

Let us define a separable multi-component system to be a system in which the energy of each individual component is dominated by its self-energy and the small interaction between components leads to thermodynamic equilibrium among different components. More specifically, we require that $V^{(i)}_{\mu\nu}\approx J_{ii}Q^{(i)}_{\mu\nu}$ for all $i$ in Eq. (15), which is typically191919This does not hold if $Q^{(i)}_{\mu\nu}\lesssim\sum_{j}N_{j}J_{ij}Q^{(j)}_{\mu\nu}/(N_{i}J_{ii})$ for some $i$. satisfied if

[TABLE]

In such a system each individual component is subject to the one-component model developed above and thus the analysis (including all of the figures) applies to each one of them. In the special case where the mutual interaction of a particular component $i$ with all other components is exactly zero then this component behaves as an isolated system and it may have a different temperature $T_{i}$ from the rest of the system in equilibrium in a microcanonical ensemble. Otherwise, if the mutual interactions between components are small but non-zero, then different components settle at a common temperature $T_{i}=T$ and $\gamma_{i}=\gamma$. However note that for each component, the equilibrium distribution function (Eq. 45) depends explicitly on the dimensionless quantities

[TABLE]

as shown in the figures above. Thus, different components may exhibit different phases (uniaxial, biaxial, or disordered) in equilibrium, depending on their respective values of $J_{ii}N_{i}$ and $c_{i}$ at the same $T$ and $\gamma$.

In the special case where $J_{ii}N_{i}=J_{jj}N_{j}$ for all $i$ and $j$ and either $L=0$ or $l_{i}=l_{j}$ for all $i$ and $j$, all components’ phase transition occurs at the same temperature. Then depending on the microscopic initial conditions, some components lie in the ordered and the rest in the disordered phase. In this sense, during the phase transition, different phases coexist.

In multi-component separable systems with radially non-overlapping components, those with larger values of $J_{ii}N_{i}\propto N_{i}m_{i}^{2}/a_{i}$ have smaller $\tau_{i}$ and therefore form thinner disks than components with lower $N_{i}m_{i}^{2}/a_{i}$ (see Figures 3 and 7). In particular, if the radial number density follows $n\propto a^{-\gamma}$ and the stellar mass distribution is independent of radius then $N_{i}\propto a_{i}^{3-\gamma}$ and $\tau_{i}\propto N_{i}m_{i}^{2}/a_{i}\propto a_{i}^{\gamma-2}$ (see Eq. 79). Thus if the radial number density is steeper than $a^{-2}$, as it is for the massive stars in the Galactic center (Bartko et al., 2009; Yelda et al., 2014), then the fractional thickness of the equilibrium disk increases outwards (i.e., the disk flares) and the distribution may be nearly isotropic (disordered) beyond a transition radius. Conversely, if the number density is shallower than $a^{-2}$ then the disk thickness increases inwards and the system is nearly isotropic inside some transition radius. Similarly, higher mass objects form thinner disks. We emphasize however that the multi-component system observed in the Galactic center is probably non-separable and higher order multipoles beyond the quadrupole are likely to play a significant role in the dynamics. We leave further study of multi-component systems to future work.

8. Conclusions

In self-gravitating systems subject to a dominant central potential, such as that of a massive central object, the bodies typically follow bounded planar orbits about the center. Due to rapid in-plane precession, each planar orbit can be represented by an axisymmetric surface density profile (over timescales $t$ given by $t_{\rm in-plane}\ll t\ll t_{\rm vrr}$, see Eq. 1). On longer timescales the angular-momentum directions, which describe the orientation of the orbital planes, relax and attain thermal equilibrium long before other degrees of freedom do (see Section 1). This process is called vector resonant relaxation (VRR). We determined the mean-field thermodynamic VRR equilibrium states using the quadrupole approximation for the gravitational interactions between orbits. For a list of our assumptions see the introduction to Section 2. The equilibria exhibit a remarkable variety of behavior including three qualitatively different VRR-phases:

1. (i)

a uniaxial “nematic” ordered phase which represents a disk in physical space containing both prograde and retrograde orbits, surrounded by a dilute halo (Figures 12 and 13); 2. (ii)

a biaxial “nematic” ordered phase consisting of two disks in physical space; the two disks have the same thickness and mass and their mutual inclination is between $90^{\circ}$ and $180^{\circ}$ (Figure 18); 3. (iii)

an axially symmetric disordered phase representing a nearly isotropic distribution of orbits in both physical space and angular-momentum space (Figure 12).

The term “nematic” highlights the analogy with liquid crystals, where the molecules specified by their symmetry axes $\bm{n}_{i}$ are concentrated towards a symmetry axis $\hat{\bm{s}}$, in the sense that $\bm{n}_{i}$ is preferentially parallel or antiparallel to $\hat{\bm{s}}$. In both liquid crystals and the quadrupole VRR model the Hamiltonian is invariant under the inversion $\bm{n}_{i}\to-\bm{n}_{i}$, but the direction of the total angular momentum $\bm{L}$ in the VRR model breaks this symmetry.

The ordered and disordered phases are stable at low and high temperatures, respectively202020See definition of temperature below Eq. (19.. The biaxial ordered phase, consisting of two highly inclined disks, requires the temperature to be limited to the range shown in Figure 21 and the mean angular momentum to be limited to $L\leq Nl/\surd 2$ where $l$ is the angular momentum of one object (Eq. 2). Biaxial states are metastable in both the canonical and microcanonical ensembles (i.e., in the canonical [microcanonical] ensemble they are a local minimum of the VRR free energy [negative of the entropy] but this minimum is larger than that of axisymmetric ordered states with the same temperature [energy]). Nevertheless, for practical purposes the lifetime of the two-disk state may be very long, especially if the disks are thin.

The canonical ensemble exhibits a first-order gravitational phase transition between the axisymmetric ordered and disordered states for a limited range of total angular momentum $L<0.137\,Nl$ (where $l$ is the magnitude of the angular-momentum vector of a single object, see Figure 9). The phase-transition temperature $T_{P}$ depends on the total angular momentum ($0.146\leq kT_{P}(L)/(JN)\leq 0.151$). The disordered state is only stable if either $T>T_{P}$ or $T<0$. Negative-temperature states are allowed due to the fact that the VRR energy has a maximum possible value, which corresponds to the isotropic distribution for $L/Nl<1/\surd 3$ (see Section 3.2 and Figures 16 and 20).

At a critical total angular momentum $L_{\rm cr}=0.137Nl$ the gravitational phase transition becomes second-order in the canonical ensemble and there is a smooth crossover between the ordered and disordered states for larger $L$. For high angular momentum, equilibria resemble an axisymmetric disk in physical space for any positive temperature.

Due to the non-additivity of the quadrupolar mean-field model and the existence of negative heat-capacity equilibrium states, the canonical and microcanonical ensembles are inequivalent. In particular, the phase-transition of the canonical ensemble is replaced by a stable sequence of equilibria without a phase transition in the microcanonical ensemble including a region of negative specific heat (between $A$ and $B$ in Figure 11). Furthermore, negative heat capacity biaxial equilibria (between $I$ and $G$ in Figure 19) are metastable in the microcanonical ensemble but unstable in the canonical ensemble. Outside these regions, all VRR equilibrium states are identical in the microcanonical and canonical ensembles.

We furthermore introduced an ensemble, the $\omega TN$-ensemble, in which the system is embedded in a much bigger bath with which it can exchange not only VRR energy but also angular momentum. This ensemble also exhibits a first-order phase transition and a critical point where it becomes second-order, but at different values of the thermodynamic variables (see Section 6).

The one-component quadrupolar VRR model discussed in this paper is reminiscent of the Hamiltonian Mean Field (HMF) model (Dauxois et al., 2002), although the VRR model has two degrees of freedom per body while the HMF model has only one, and the HMF model has a kinetic energy term in the Hamiltonian where VRR does not. Both the VRR and the HMF model exhibit ensemble inequivalence, first- and second-order phase transitions, and negative heat capacity, but the biaxial metastable equilibria and the negative-temperature equilibria do not appear in the HMF model.

There are strong similarities between the equilibrium configurations of VRR and liquid crystals. The Newtonian gravitational interaction between rapidly precessing elliptical orbits, which trace out axisymmetric punctured disks, is similar to the Coulomb interaction between axisymmetric molecules. In particular, the quadrupole mean-field Hamiltonian of VRR for a one-component system with zero total angular momentum is equivalent to the Maier–Saupe model of liquid crystals (Maier & Saupe, 1958), which exhibits a nematic-isotropic phase transition. In the nematic phase of VRR, the bodies are configured as a disk embedded in a dilute halo. The onset of the nematic-isotropic transition in VRR depends on the total angular momentum, just as the onset of the transition in liquid crystals depends on the external magnetic field, as can easily be verified by comparing our Figure 9 with the corresponding ones of Wojtowicz & Sheng (1974); Gramsbergen et al. (1986) for example.

However, there are important differences between the quadrupolar VRR system and liquid crystals. In contrast to the diamagnetic term $-{\textstyle\frac{1}{2}}\Delta\chi(\bm{B}\cdot\bm{n}_{i})^{2}$ that arises in the free energy due to an external magnetic field in liquid crystals, the term that arises in VRR due to the angular-momentum constraint is $-l\bm{\omega}\cdot\bm{n}_{i}$, which resembles a paramagnetic term of a spin system in an external magnetic field, $-\chi\bm{B}\cdot\bm{n}_{i}$.

We have restricted our attention in this paper to the mean-field approximation, in which each body is drawn independently from a distribution function and correlations are ignored. In the limit of zero temperature ($T\rightarrow 0^{+}$), we find two distinct distributions that are in stable or metastable thermodynamic equilibrium: (i) a single razor-thin disk containing bodies with orbit normals that are both aligned and anti-aligned to the total angular momentum, with the relative numbers of each population depending on the total angular momentum (state $O_{1}$); and (ii) two razor-thin disks of equal mass, with normals inclined to the total angular-momentum vector by ${\textstyle\frac{1}{2}}\iota$ where $\iota$ depends on the total angular momentum and lies between $90^{\circ}$ and $180^{\circ}$ set by $L/(Nl)$ (state $O_{5}$). We call these states “uniaxial" and “biaxial", respectively.

The zero-temperature equilibria $O_{1}$ and $O_{5}$ are discrete configurations of angular-momentum vector directions $\bm{n}$ (i.e., the distribution function has compact support on the unit sphere) that are time-invariant up to a rigid-body rotation at uniform angular speed. These are not the only such configurations. When $L=0$ other examples include (i) $\bm{n}$ oriented along the vertices of a regular polyhedron, (ii) $\bm{n}$ oriented along the vertices of a regular planar polygon in an arbitrary plane, (iii) $\bm{n}$ oriented along the vertices of a regular planar polygon and either or both of the two directions perpendicular to the polygon’s plane. In the simplest of these configurations the mass in each component is equal. These configurations do not appear in our maximum-entropy analysis, but in the limit where the number of components is large, we may recover the axisymmetric continuous zero-temperature states derived earlier in this paper: for example, configuration (ii) $\rightarrow O_{2}$.

Our analysis is based on two important simplifications in addition to the mean-field approximation (see Section 2 for a detailed inventory of our simplifying approximations). First, we approximated the gravitational torques between orbit annuli by their quadrupole component. Multipoles beyond the quadrupole may be important for systems with radially overlapping or closely spaced orbits. Second, we focused on a one-component model, in which all bodies have the same scalar angular momentum and all of the pairwise interactions have the same coupling coefficients (we briefly discussed a simple extension of this model to multi-component systems in Section 7). As we have shown in this paper, the statistical mechanics of this simplified model can be described completely yet exhibits a remarkable range of behavior. We expect that many of the features we have observed will also be present in more general models of stellar systems dominated by a central mass that do not depend on these approximations.

We are grateful to Julien Barré for generously helping us to clarify several aspects of the stability analysis and to Benjamin Beri for pointing out that the quadrupole mean-field Hamiltonian is described by the Maier-Saupe model of liquid crystals, and to the referee for thoughtful comments that substantially improved the paper. We thank Zoltan Rácz for useful discussions. This work was supported in part by the European Research Council under the European Union’s Horizon 2020 Programme, ERC Starting Grant #638435 (GalNUC), by the U.S. National Science Foundation through grant AST-1406166, and by NASA through grant NNX14AM24G.

Appendix A Hamiltonian for vector resonant relaxation

Here we describe the interaction Hamiltonian that governs VRR, that is, the interaction energy between Keplerian orbits after averaging over orbital phase and apsidal precession (equivalently, mean anomaly and longitude of pericenter). The most general form of this Hamiltonian may be written (Kocsis & Tremaine, 2015)

[TABLE]

where $P_{\ell}(x)$ are Legendre polynomials and $p$ and $q$ label the $N$ bodies, which orbit around the central point mass with angular-momentum unit vectors $\bm{n}_{i}$. In this formula $\cos\theta_{pq}=\bm{n}_{p}\cdot\bm{n}_{q}$,

[TABLE]

and

[TABLE]

Here $m_{p}$, $a_{p}$, $e_{p}$ denote the mass, semimajor axis, and eccentricity of the orbit of body $p$ around the central object, all of which are fixed during VRR, and “${\rm out}$” and “${\rm in}$” label the index $p$ or $q$ with the larger and the smaller semimajor axis, respectively, i.e., $\alpha_{pq}\leq 1$. Note that only terms with even $\ell$ contribute to the sum, and the $\ell=0$ term only contributes an unimportant constant to the Hamiltonian.

For circular orbits $s_{pq\ell}=1$ for all $\ell$, and more generally, for eccentric orbits with $a_{\rm out}(1-e_{\rm out})>a_{\rm in}(1+e_{\rm in})$ (radially non-overlapping orbits)

[TABLE]

for $\ell>0$. Here $\chi_{p}=a_{p}/b_{p}=(1-e_{p}^{2})^{-1/2}$. For radially non-overlapping orbits the sum in Eq. (A1) converges exponentially as a function of $\ell$ and the $\ell=2$ term dominates the dynamics for arbitrary mutual inclinations. In contrast, for radially overlapping orbits, $a_{\rm out}(1-e_{\rm out})\leq a_{\rm in}(1+e_{\rm in})$, the contribution of the $\ell^{\rm th}$ multipole decays as $\ell^{-2}$, and multipoles up to order $\ell_{\max}\sim(1-\bm{n}_{p}\cdot\bm{n}_{q})^{-1/2}$ must be included when calculating the Hamiltonian (see Figure B1 in Kocsis & Tremaine 2015). However, if most orbits have relatively large inclinations and $N\gg 1$, the net torque on an object is still dominated by the $\ell=2$ quadrupole interaction. In particular, for a spherical cluster the contribution of the $\ell^{\rm th}$ multipole to the net torque decays as fast as $\ell^{-3}\ln\ell$ (see Appendix D in Kocsis & Tremaine 2015). For clusters with overlapping orbits, the effects of multipoles beyond the quadrupole may be significant.

For simplicity, let us keep only the $\ell=2$ term and assume that a large number $N_{i}$ of bodies have similar $(m_{i},a_{i})$ but different $\bm{n}$ with several such components $i=1,2,\dots,K$. Then the interaction energy is a sum over the components and the bodies therein,

[TABLE]

where $J_{ij}={\textstyle\frac{3}{2}}P_{2}(0)^{2}C_{ij}s_{ij2}\alpha_{ij}^{2}={\textstyle\frac{3}{8}}C_{ij}s_{ij2}\alpha_{ij}^{2}$ and $g(\bm{n}_{p},\bm{n}_{q})={\textstyle\frac{2}{3}}P_{\ell}(\bm{n}_{p}\cdot\bm{n}_{q})=(\bm{n}_{p}\cdot\bm{n}_{q})^{2}-\frac{1}{3}$. In the second line the sums over $p$ and $q$ of components $i$ and $j$ respectively are written as integrals using the corresponding distribution functions $f_{i}$ and $f_{j}$ of the angular-momentum unit vectors. For circular orbits, $s_{ij2}=1$, and this is Eq. (8) used in the main text. For eccentric, radially overlapping orbits $s_{ij2}$ is given analytically by Eqs. (B12) and (B31) in Kocsis & Tremaine (2015).

Appendix B Alignment of total angular momentum

In this appendix we show that, at equilibrium, the angular momentum and the Lagrange multiplier $\bm{\gamma}$ lie along one of the eigenvectors of the matrix $Q_{\mu\nu}$. We work in a coordinate system aligned with these eigenvectors so $Q_{\mu\nu}$ is diagonal (Eq. 3) and we can write $Q_{\mu\nu}\equiv Q_{\mu}\delta_{\mu\nu}$ (there is no summation over $\mu$ here, that is, $Q_{\mu\nu}\equiv\mbox{diag}[\mathbf{Q}]$, where $\mathbf{Q}$ is the vector whose elements are the eigenvalues of $Q_{\mu\nu}$). The self-consistency condition (28) is then

[TABLE]

In the case $\mu=1$, $\nu=2$ this condition simplifies to

[TABLE]

in the last equation we have introduced a spherical polar coordinate system. Now divide the interval $\phi\in[0,2\pi)$ into four quadrants: in the first quadrant, $0\leq\phi<\pi/2$, use $\psi=\phi$ as the integration variable; in the second quadrant, $\pi/2\leq\phi<\pi$, use $\psi=\pi-\phi$; in the third quadrant, $\pi\leq\phi<3\pi/2$, use $\psi=\phi-\pi$; and in the fourth quadrant, $3\pi/2\leq\phi<2\pi$, use $\psi=2\pi-\phi$. Adding the sub-integrals we find

[TABLE]

Since the integration range is from 0 to $\pi/2$, $\sin\psi$ and $\cos\psi$ are positive, so the integral is positive-definite if $\gamma_{1}\gamma_{2}>0$, negative-definite if $\gamma_{1}\gamma_{2}$ is negative, and zero (as required) if and only if at least one of $\gamma_{1}$ and $\gamma_{2}$ is zero. Repeating this argument for $\mu=1$, $\nu=3$ and $\mu=2$, $\nu=3$ we conclude that at least two of the $\gamma_{\mu}$ must be zero. Therefore $\bm{\gamma}$ must be aligned with one of the coordinate axes and thus with one of the eigenvectors of $Q_{\mu\nu}$.

Without loss of generality we may assume that $\bm{\gamma}$ points along $\bm{n}_{3}$. Then the angular-momentum vector is

[TABLE]

If $\mu=1$ or 2, the integrand in the numerator is odd under the transformation $\phi\to\phi+\pi$ so the integral must vanish. Thus only $L_{3}$ is non-zero so the angular momentum $\bm{L}$ is parallel or anti-parallel to $\bm{\gamma}$.

Appendix C Second-order variation

Here we calculate the second variation of entropy and free energy under general perturbations of the mean-field distribution function and show the stability of the negative-temperature configurations and the disordered phase at infinite temperature, labeled $D$ in Figure 3.

Let $f$ be the equilibrium distribution function,

[TABLE]

and let $\delta f$ be a perturbation about this equilibrium. The entropy changes to

[TABLE]

Expand this to second order in $\delta f$. Then the change in entropy is

[TABLE]

Substituting Eq. (C1) and dropping terms higher than second order,

[TABLE]

Conservation of the number of bodies and the total angular momentum implies that

[TABLE]

Thus the equation for the entropy variation simplifies to

[TABLE]

The energy change can be written as

[TABLE]

which implies212121We denote $M_{\mu\nu}^{2}=\sum_{\mu\nu}M_{\mu\nu}M_{\mu\nu}$.

[TABLE]

where the second line follows from the relation $\epsilon=-q_{\mu\nu}\int\,f^{\prime}q_{\mu\nu}^{\prime}d\Omega^{\prime}$ (Eq. 15). Thus if the energy is fixed

[TABLE]

Substituting in Eq. (C6) to eliminate the terms depending on $\varepsilon$, we find

[TABLE]

where $\delta^{2}S$ is the second-order variation of entropy

[TABLE]

As it should, the entropy variation (C10) contains no first-order terms, since the perturbation is performed about the equilibrium distribution, at which the entropy has vanishing first-order variation by construction, if the constraints are satisfied (see Section 2.2).

In the canonical ensemble the equilibria are extrema of the free energy $F=E-TS$ at constant temperature $T$. Using Eqs. (C6) and (C9) we find the second-order variation of free energy to be

[TABLE]

The condition for stability of an equilibrium with respect to any perturbation is given in the microcanonical ensemble by the condition

[TABLE]

that is, the entropy must be a local maximum. In the canonical ensemble the stability condition is

[TABLE]

that is, the free energy must be a minimum if the temperature is positive, or a maximum if the temperature is negative.

Now, since $\delta^{2}S=-k\beta\delta^{2}F$, as is evident from (C11) and (C12), stability in the microcanonical ensemble is determined by the same inequality condition as in the canonical ensemble (Eqs. C13 and C13)). However, equilibria do not necessarily possess the same stability properties in the two ensembles, because the perturbations are subject to different constraints. In the canonical ensemble, the perturbations are subject only to constraints (C5), while in the microcanonical ensemble they are subject additionally to the fixed energy constraint (C9). Thus, modes that render the canonical ensemble unstable may not exist in the microcanonical ensemble. If this happens, the system presents inequivalence of ensembles, as described in Section 2.3.

Negative-temperature states

For negative temperature, $\beta<0$, both of the terms in the entropy second variation (C11) are negative-definite and so $\delta^{2}S<0$ for perturbations that keep $E$, $N$, and $\bm{L}$ fixed. Therefore negative-temperature equilibria are local maxima of the entropy and hence they are (meta)stable in the microcanonical ensemble.

Similarly, in the canonical ensemble both terms in (C12) are positive-definite and so $\beta\delta^{2}F>0$ for perturbations which keep $T$, $N$, and $\bm{L}$ fixed222222The temperature $T$ is defined by the angular momentum vector distribution function of the heat bath; it is constant because the heat bath has a much larger angular momentum than the system we are examining. Note that fixing $T$ does not impose any constraint on $\delta f$ in Eq. (C12).. Therefore negative-temperature equilibria are local minima of $\beta F$ and hence they are (meta)stable in the canonical ensemble.

State $D$

The state $D$ is the disordered phase at infinite temperature and therefore corresponds to axisymmetric solutions with $\beta\rightarrow 0^{+}$. Its stability is straightforwardly implied by Eq. (C11)

[TABLE]

Appendix D Stability of off-diagonal elements

We show here that for positive temperature, there can always be found a coordinate system with one axis aligned with $\bm{\omega}$, such that the off-diagonal elements of $Q_{\mu\nu}$ do not give rise to instabilities and are zero in equilibrium. We follow the notation of Section 3.1.

The dominant contributions to the integral (59) for the partition function of the $\omega TN$-ensemble come from the minima of $g(\mathbf{Q},\beta,\gamma)$, which occur at a subset of the extrema where $\partial g/\partial Q_{\mu\nu}=0$. The extrema are located at $\mathbf{Q}^{\rm eq}$ given by the implicit equation

[TABLE]

here the angle brackets denote the average

[TABLE]

We now describe some properties of the matrix $\mathbf{Q}^{\rm eq}$. (i) Since $q_{\mu\nu}=q_{\nu\mu}$ we have $Q_{\mu\nu}^{\rm eq}=Q_{\nu\mu}^{\rm eq}$, that is, $\mathbf{Q}^{\rm eq}$ is symmetric. (ii) Since $\sum_{\mu=1}^{3}q_{\mu\mu}=0$, we have $\mbox{Tr\,}\mathbf{Q}^{\rm eq}=0$, that is, $\mathbf{Q}^{\rm eq}$ is traceless. (iii) Since $q_{\mu\nu}$ is a tensor (Eq. 12), Eq. (D1) implies that $\mathbf{Q}^{\rm eq}$ also transforms as a tensor under rotations. Therefore without loss of generality we can choose the coordinate system in the $\Omega$-integral so that $\mathbf{Q}^{\rm eq}$ is diagonal232323Thus systems with $Q^{\rm eq}_{xx}\not=Q^{\rm eq}_{yy}$ are an example of spontaneously broken symmetry.. The whole system is free to change its principal axes perpendicular to $\bm{L}$. (iv) Using similar arguments to those in Appendix B we may then show that $\bm{\gamma}$ lies along one of these axes, which we may choose to be the $z$-axis.

Given these results, all of the elements of $\mathbf{Q}^{\rm eq}$ are either zero or linear combinations of the parameters $Q$ and $W$ defined in Eqs. (3)–(36). It is therefore useful to replace the nine coordinates $Q_{\mu\nu}$ by nine new coordinates $\bm{P}=(Q,W,X_{0}\ldots,X_{3},Y_{1},\ldots,Y_{3})$ defined by

[TABLE]

The Jacobian $|\partial\bm{P}/\partial\mathbf{Q}|=\frac{1}{12}$. In the new coordinates

[TABLE]

Here $q=\cos^{2}\theta-\frac{1}{3}$ and $w=\sin^{2}\theta\cos 2\phi$, as in Eqs. (35) and (36). It is straightforward to confirm that the extrema of $G$ occur at $X_{0}=X_{1}=\cdots=Y_{3}=0$. An extremum is a minimum if the Hessian $\partial^{2}G/\partial P_{\mu}\partial P_{\nu}$ is positive definite (i.e., all of its eigenvalues are positive) and a maximum if it is negative definite; otherwise the extremum is a saddle point. It is straightforward to confirm that all off-diagonal elements of the Hessian are zero in the 7 rows and columns corresponding to $(X_{0},X_{1},\ldots,Y_{3})$ and that

[TABLE]

Thus $G$ is a minimum at the extrema as a function of these seven coordinates.

Appendix E Analytic results for the equilibria

Here we derive analytic expressions for the thermal equilibria. We define the moment generating function

[TABLE]

where $q=\cos^{2}\theta-{\textstyle\frac{1}{3}}$, $w=\sin^{2}\theta\cos 2\phi$, and $s=\cos\theta$, the integral is over the unit sphere, and the parameters $\kappa_{1}$, $\kappa_{2}$, and $c$ are defined by Eq. (53). Once $Z_{0}(\kappa_{1},\kappa_{2},c)$ is known, all statistical quantities follow straightforwardly242424We suppress the $N\ln N$ constants from $S$, $F$, and $G$.:

[TABLE]

[TABLE]

and for the $\omega TN$-ensemble we have

[TABLE]

Eqs. (E9) follow from Eq. (53) assuming $\left\langle q\right\rangle\neq 0$ and $\left\langle w\right\rangle\neq 0$. In the special case $\left\langle q\right\rangle=\left\langle w\right\rangle=0$, self-consistency requires that $c=0$. In this case $Z_{0}=4\pi$ is a trivial solution for arbitrary $T$, and $E=L=0$. This is the isotropic distribution $f(\bm{n})=N/(4\pi)$, which is an equilibrium for any temperature and is shown by the curve $O_{3}D$ in Figure 3.

E.1. Axisymmetric equilibria with zero total angular momentum

Systems that are axisymmetric must have two of the diagonal elements of $Q_{\mu\nu}$ (Eq. 3) equal. This requires that either $\left\langle w\right\rangle=0$ or $\left\langle w\right\rangle=\pm 3\left\langle q\right\rangle$. Since the total angular momentum $\bm{L}=0$ we are free to choose the $z$-axis to be aligned with any one of the three eigenvectors of $Q_{\mu\nu}$ so without loss of generality we can assume that $\left\langle w\right\rangle=0$. Then $\kappa_{2}=0$ by Eq. (53) and $c=0$ by Eq. (55). Eq. (E1) simplifies to

[TABLE]

Note that this is real for either positive or negative $\kappa_{1}$. Eqs. (E2) and (E9) yield the parametric solution

[TABLE]

As $\kappa_{1}$ increases from $-\infty$ to $\infty$, $\left\langle q\right\rangle$ increases monotonically from $-\frac{1}{3}$ to $\frac{2}{3}$. Both limits $\kappa_{1}\to\pm\infty$ correspond to $T\to 0^{+}$ so the two limits are the zero-temperature states $O_{2}$ and $O_{1}$ in Figure 3. As $\kappa_{1}$ increases from $-\infty$ the temperature grows monotonically from $T=0$ to $kT_{A}/(JN)=0.148556$ at $\kappa_{1}=2.178289$, then it decreases monotonically for higher $\kappa_{1}$, approaching zero at $\kappa_{1}\rightarrow\infty$. At $T_{A}$, $\left\langle q\right\rangle_{A}=0.21573$, Point $B$ in Figure 3 has $\kappa_{1}\rightarrow 0$, $\left\langle q\right\rangle_{B}=0$, and $T_{B}=\frac{2}{15}$.

E.2. Axisymmetric, rotating equilibria

Next consider systems with non-zero angular momentum. Since $Z_{0}$ is an even function of $\kappa_{2}$, $\kappa_{2}=0$ implies $\left\langle w\right\rangle=0$ (see Eq. E3). These states are axisymmetric around the angular-momentum axis with $0\leq\left\langle s\right\rangle<1$ and $\left\langle q\right\rangle\neq 0$. We evaluate $Z_{0}$ (E1) at $\kappa_{2}=0$

[TABLE]

which is real for all $-\infty<\kappa_{1}<\infty$. The quantities $\left\langle s\right\rangle$, $\left\langle q\right\rangle$, and $T$ may be obtained from Eqs. (E4), (E2), and (E9), which yields the parametric solution

[TABLE]

For any fixed $-\infty<\kappa_{1}<\infty$, Eq. (E16) defines a monotonic function of $0<c<\infty$ mapping onto $L/(Nl)\equiv\left\langle s\right\rangle\in[0,1)$. For any fixed $0<\left\langle s\right\rangle<1$, the temperature assumes all values from $-\infty$ to $\infty$ since the expression for $T$ is singular at $\kappa_{1}\rightarrow 0$ (state $D$). In the canonical and microcanonical ensembles with fixed $\left\langle s\right\rangle$, $\left\langle q\right\rangle$ is a monotonically increasing function of $\kappa_{1}$, but the temperature has three local extrema as a function of $\kappa_{1}$ for $\left\langle s\right\rangle<0.13714$, one local maximum for $0.13714<\left\langle s\right\rangle<1/\sqrt{3}$ and no extrema for $1/\sqrt{3}<\left\langle s\right\rangle<1$. Note that $T\rightarrow 0$ at $\kappa_{1}\rightarrow\infty$ (state $O_{1}$), $\kappa_{1}\rightarrow-\infty$ (state $O_{2}$), and at a finite $\kappa_{1}$ for $\left\langle s\right\rangle<1/\sqrt{3}$ (state $O_{3}$). Note that $\left\langle q\right\rangle=0$ at $O_{3}$ for $\left\langle s\right\rangle<1/\sqrt{3}$ due to Eq. (E17).

In the $\omega TN$-ensemble, $\omega=kcT(\kappa_{1},c)/l$ and Eq. (E18) defines a monotonic function of $c$ for fixed $\kappa_{1}>0$, which can be inverted numerically to find $c(\kappa_{1},\omega)$ for the equilibrium states between $O_{1}$ and $D$ for fixed $\omega$. Substituting in Eqs. (E17) and (E18), $\left\langle q\right\rangle$ is a monotonic function of $\kappa_{1}$ for fixed $0<\omega<\infty$, while $T$ has two local extrema as a function of $\kappa_{1}$ for $0<l\omega/JN<0.057285$ and no local extrema for $l\omega/JN>0.057285$. The heat capacity is negative between the two local maxima and positive otherwise.

Next we provide asymptotic expressions near $O_{2}$, $O_{1}$, and $D$. We utilize the asymptotics of the error function

[TABLE]

In particular as $\kappa_{1}\rightarrow-\infty$ and $c\rightarrow\infty$ with fixed $c/\kappa_{1}$ (state $O_{2}$ of Figure 7)252525The expressions below contain factors of $\sigma=1+(c/2\kappa_{1})^{2}$. The factor $c/(2\kappa_{1})$, which is proportional to $\omega$ defined in Eq. (20), is retained because at fixed angular momentum $\left\langle s\right\rangle$, $c\propto\kappa_{1}$ as $\kappa_{1}\to-\infty$ according to Eq. (E23).

[TABLE]

which imply that

[TABLE]

and for $\kappa_{1}\rightarrow\infty$ and $0\leq\lim_{\kappa_{1}\rightarrow\infty}c<\infty$ (state $O_{1}$ of Figure 7)

[TABLE]

As $\kappa_{1}\to 0$ (state $D$ of Figures 15, 16, and 20) the temperature diverges,

[TABLE]

These results show that for axisymmetric states $\left\langle q\right\rangle$ generally increases from $\left\langle s\right\rangle^{2}-1/3$ to $2/3$ as $\kappa_{1}$ changes from $-\infty$ to $\infty$, the endpoints being the $O_{2}$ and $O_{1}$ states (see Figure 7). The temperature approaches zero at both of these endpoints. For very small negative $\kappa_{1}$, the temperature is positive for $\left\langle s\right\rangle<1/\sqrt{3}$ and negative for $\left\langle s\right\rangle>1/\sqrt{3}$. For intermediate values of $\kappa_{1}$ (these are not shown by the asymptotics), we find that $\left\langle q\right\rangle$ increases monotonically as a function of $\kappa_{1}$ for fixed $\left\langle s\right\rangle$ and the temperature assumes all values between $-\infty$ and $\infty$ with up to three local maxima as a function of $\kappa_{1}$, i.e., one with $\left\langle q\right\rangle<0$ for $\left\langle s\right\rangle<1/\sqrt{3}$ and two with $\left\langle q\right\rangle>0$ for $\left\langle s\right\rangle<\left\langle s\right\rangle_{\rm cr}=0.13714$. In the limit $\left\langle s\right\rangle\rightarrow 0$, two of the local maxima of $T(\kappa_{1})$ approach points $B$ and the third approaches point $A$ in Figure 3.

E.3. Non-axisymmetric rotating equilibria

For arbitrary $\kappa_{2}$ one of the integrals over the two polar angles may be evaluated in $Z_{0}$ Eq. (E1), which gives (see Eq. 50)

[TABLE]

The zero-temperature non-axisymmetric equilibria correspond to $\kappa_{2}\rightarrow\infty$. We use

[TABLE]

for $z\equiv\kappa_{2}(1-s^{2})\rightarrow\pm\infty$. We obtain

[TABLE]

The asymptotics near state $O_{5}$ in Figure 17 may be obtained from Eq. (E34) in the limit $\Delta\kappa\equiv\kappa_{2}-\kappa_{1}\rightarrow\infty$ and $c\rightarrow\infty$ such that

[TABLE]

is fixed. Laplace’s method262626 For any twice differentiable function $g(s)$ with a unique minimum $\lambda$ in $[a,b]$, the following integral (if it exists) approaches

$\displaystyle\int_{a}^{b}f(s)e^{-xg(s)}ds\overset{x\rightarrow+\infty}{\longrightarrow}\frac{\sqrt{2\pi}}{\sqrt{ag^{\prime\prime}(\lambda)}}e^{-xg(\lambda)}\left[f(\lambda)+\mbox{O}(x^{-1})\right]\,.$

(E36)

We apply Eq. (E36) in Eq. (E34) with $x=\Delta\kappa$, $f(s)=(1-s^{2})^{-1/2}$, and $g(s)=(s-\lambda)^{2}$.

yields

[TABLE]

Substituting in Eqs. (E2)–(E4) gives

[TABLE]

Substituting in the self-consistency equation (E9) gives a relation between $\Delta\kappa$ and $\kappa_{2}$

[TABLE]

We may now eliminate $\Delta\kappa$. Next solve for $\lambda$ using Eq. (E38) and substitute back into Eqs. (E9) and (E39)–(E40) to get the asymptotics near state $O_{5}$ parameterized by $\kappa_{2}$ for any given $\left\langle s\right\rangle$:

[TABLE]

The $O_{5}$ state corresponds to the limit $\kappa_{2}\rightarrow\infty$.

The asymptotics near $O_{4}$ may be obtained from Eq. (E34) in the limit that $\Delta\kappa$ approaches a finite value while $\kappa_{2}\rightarrow\infty$. In practice, we find numerically that $0\leq\Delta\kappa<1$ for $0\leq\left\langle s\right\rangle<0.52$, $\lim_{O_{4}}\Delta\kappa=2\left\langle s\right\rangle^{2}$ for $0\leq\left\langle s\right\rangle\ll 0.5$, and for any fixed $0<\left\langle s\right\rangle<2^{-1/2}$, $c$ approaches a finite value. We derive an analytic approximation for the asymptotic behavior for $\left\langle s\right\rangle<0.52$ by expanding $Z_{0}$ Eq. (E34) in $\Delta\kappa$ around 0 and using the identity (46),

[TABLE]

where $I_{n}\equiv I_{n}(c)$ is the modified Bessel function272727These equations become inaccurate for $s\gtrsim 0.52$. For $0.65<s<2^{-1/2}$, we find numerically that $\Delta\kappa>10$, which implies that accurate analytic expressions exist in this regime (not shown), which may be derived with the Laplace method as in Eq. (E36).. Similarly, from Eqs. (46), (E34), and (E2)–(E4), we get

[TABLE]

We may substitute in the self-consistency equation (E9) to get a relation between $\Delta\kappa$ and $\kappa_{2}$

[TABLE]

Now eliminate $\Delta\kappa$ from Eqs. (E48)–(E50)

[TABLE]

[TABLE]

Therefore at zero temperature, the $O_{2}$ and $O_{4}$ equilibria depend on $\left\langle s\right\rangle$, while $O_{1}$ and $O_{5}$ are independent of $\left\langle s\right\rangle$. The latter series of equilibria depend on $\left\langle s\right\rangle$ for $T\neq 0$ at first and second beyond leading order in $T$ for $O_{5}$ and $O_{1}$, respectively. In the $T\rightarrow 0$ limit, Eqs. (E43)–(E46) show that the $O_{5}$ order parameters are in the range $-\frac{1}{3}\leq\left\langle q\right\rangle\leq\frac{1}{6}$ and $1\geq\left\langle w\right\rangle\geq\frac{1}{2}$ depending on $\left\langle s\right\rangle$ as long as it satisfies $0\leq\left\langle s\right\rangle\leq 2^{-1/2}$. Outside of this range, $\Delta\kappa<0$ for all $\kappa_{2}$, hence the $\Delta\kappa\rightarrow\infty$ assumption cannot be satisfied, and there is no $O_{5}$ state. The asymptotics near $O_{4}$ (Eqs. E52)–(E55) show that $\left\langle q\right\rangle=\frac{1}{6}$ and $\left\langle w\right\rangle=\frac{1}{2}$, independent of $\left\langle s\right\rangle$.

In the $\left\langle s\right\rangle\rightarrow 0$ limit, the non-axisymmetric equilibria near $O_{4}$ and $O_{5}$ (Eqs. E52–E55 and E43–E46) reduce to the axisymmetric asymptotics near $O_{2}$ and $O_{1}$ (Eqs. E20–E22 and E26–E28), respectively, in a rotated coordinate system (see Eq. 3)

[TABLE]

The energy of non-axisymmetric equilibria is bounded between

[TABLE]

where

[TABLE]

and $E_{G}\leq E_{O_{3}}$ for all $L$ where $E_{O_{3}}$ is the upper energy bound in Eq. (44) of the main text (see Figure 19).