The statistical mechanics of planet orbits

# The statistical mechanics of planet orbits

Scott Tremaine11affiliationmark: Institute for Advanced Study, Princeton, NJ 08540
###### Abstract

The final “giant-impact” phase of terrestrial planet formation is believed to begin with a large number of planetary “embryos” on nearly circular, coplanar orbits. Mutual gravitational interactions gradually excite their eccentricities until their orbits cross and they collide and merge; through this process the number of surviving bodies declines until the system contains a small number of planets on well-separated, stable orbits. In this paper we explore a simple statistical model for the orbit distribution of planets formed by this process, based on the sheared-sheet approximation and the ansatz that the planets explore uniformly all of the stable region of phase space. The model provides analytic predictions for the distribution of eccentricities and semimajor axis differences, correlations between orbital elements of nearby planets, and the complete N-planet distribution function, in terms of a single parameter, the “dynamical temperature”, that is determined by the planetary masses. The predicted properties are generally consistent with N-body simulations of the giant-impact phase and with the distribution of semimajor axis differences in the Kepler catalog of extrasolar planets. A similar model may apply to the orbits of giant planets if these orbits are determined mainly by dynamical evolution after the planets have formed and the gas disk has disappeared.

planets and satellites: dynamical evolution and stability — planets and satellites: formation — celestial mechanics

## 1 Introduction

It is always tempting to apply the powerful tools of statistical mechanics to macroscopic physical systems that exhibit some degree of regularity. The first hint of a role for statistical mechanics in planet-formation theory came from long-term integrations of the solar system. These showed that (i) the orbits of all the planets in the solar system are chaotic, with Liapunov or e-folding times of a few Myr (Sussman & Wisdom, 1988; Laskar, 1989; Sussman & Wisdom, 1992); (ii) the outer solar system, between Jupiter and Neptune, is “full” in the sense that there are almost no stable orbits for test particles in this region (Holman & Wisdom, 1993; Holman, 1997); (iii) there is a 1–2% probability that chaotic diffusion of Mercury’s eccentricity will lead to its loss—by ejection, collision with the Sun, or collision with another planet—within the next 5 Gyr (Laskar & Gastineau, 2009).

In the words of Laskar (1996), these findings lead to the speculation that “maybe there was some extra planet at the early stage of formation of the solar systembut this led to so much instability that one of the planetssuffered a close encounter or a collision with the other ones. This leads eventually to the escape of the planet and the remaining system gets more stable. In this case, at each stage, the system should have a time of stability comparable with its age.” If this hypothesis is correct, the current configuration of the solar system and other planetary systems might be determined, at least in part, by the statistics of orbital chaos. Moreover this evolution process might be approximately self-similar, in that the distributions of planet masses, eccentricities, inclinations, and semimajor axis differences in an ensemble of planetary systems would remain unchanged except for scale factors as the systems evolved.

Qualitatively similar ideas have emerged in the exoplanet community, where they are typically called the “packed planetary systems hypothesis”. The simplest version of this hypothesis (Barnes & Raymond, 2004; Fang & Margot, 2013) is that most planetary systems are “as tightly packed as possible”, that is, there is no room for additional planets on stable orbits. A bold prediction of this hypothesis is that if there are stable regions of phase space between known exoplanets then there must be undetected planets in these regions.

Current theories of terrestrial planet formation involve multiple stages and processes (for recent reviews see Kokubo & Ida 2012; Morbidelli et al. 2012; Haghighipour 2013; Raymond et al. 2013). As the gaseous protoplanetary disk cools, dust condenses and settles into the disk midplane; the dust particles then accumulate into planetesimals; and the largest planetesimals undergo runaway growth until they dominate the gravitational scattering of smaller planetesimals. At this stage a phase of self-regulated or oligarchic growth begins, in which the largest planetesimals—now called planetary embryos—all grow at similar rates. The oligarchic phase ends when the reservoir of small planetesimals is exhausted. At this point there are typically a few dozen embryos left. The surviving embryos gradually excite one another’s eccentricities until their orbits cross and they collide. In this last stage—variously called late-stage accretion, post-oligarchic growth, or the giant-impact phase—the number of surviving bodies slowly declines until we are left with a small number of planets on well-separated, stable orbits. Laskar’s hypothesis or the packed planetary systems hypothesis are roughly equivalent to the assumption that the giant-impact phase tends to produce an ensemble of planetary systems with statistically similar properties, an idea that has been advanced recently from different perspectives by Malhotra (2015), Pu & Wu (2015), and Volk & Gladman (2015).

The goal of this paper is to explore a simple model for the distribution of orbital elements in planetary systems after the giant-impact phase. The basic ansatz behind the model is that at the end of the giant-impact phase, the planets explore uniformly all of the stable phase space that is available to them. This ansatz neglects many physical processes that may play important roles in the late stages of planet formation (migration, gas drag, mean-motion and secular resonances, dynamical friction from a residual planetesimal population, hit-and-run, fragmenting, and catastrophic collisions, perturbations from exterior giant planets, etc.). The model does not describe the distribution of planetary masses following the giant-impact phase, only the distribution of orbits. We believe the model is useful because it yields clear predictions, with minimal free parameters, for several of the observables in multi-planet systems. It is, of course, no substitute for N-body simulations, but it can guide our interpretation of these simulations and does well at reproducing many of their results.

The initial stages of giant planet formation are believed to be similar to those of terrestrial planets, but the masses of giant planets are dominated by gas envelopes that are believed to accrete before the gaseous protoplanetary disk disappears, a few Myr after the birth of the star. It is possible that the distribution of giant planets continues to evolve over much longer times as the planets excite one another’s eccentricities. In this case the number of planets is whittled down mostly by ejection from the system rather than by collision as in the case of terrestrial planets (Chatterjee et al., 2008; Jurić & Tremaine, 2008)111The different outcomes arise because the Safronov number (eq. 34) is much larger for giant planets than terrestrial planets.. Despite this difference, our ansatz could also apply to giant planets if this late dynamical evolution is the primary mechanism that determines the distribution of their orbits.

## 2 The model

### 2.1 The planetary system

We shall work with a simplified model of a planetary system based on the sheared-sheet or Hill’s approximation (Spitzer & Schwarzschild, 1953; Hénon, 1969, 1970; Petit & Hénon, 1986; Binney & Tremaine, 2008). This simplification eliminates radial gradients in surface density, orbital angular speed, etc., which otherwise obscure the analysis.

In Hill’s approximation we focus on a small radial interval of the system centered on radius and of width . Within this interval there are orbiting masses (“planets”) with semimajor axes and eccentricities . The Hill radius of planet is

 rHi≡¯a(mi3M⋆)1/3 (1)

where is the mass of the host star. The orbital angular speed is . We assume that the planetary orbits are coplanar or nearly so. We index the planets so that their semimajor axes are ordered, . The positions of the planets are bounded by two fixed radii , , i.e., the pericenters and the apocenters . We ignore any interactions with planets outside this range.

We assume that the planet masses are small, , and (usually) that the number of planets . As grows we assume that since otherwise the dynamics is trivial: if then a system composed of planets on nearly circular orbits is stable, while if there will be frequent close encounters and collisions so the configuration is short-lived. This ordering condition can be satisfied either by letting the planetary masses shrink as for fixed radial interval (Hill’s approximation) or by fixing the planetary masses and letting the radial width grow as (the sheared-sheet approximation).

### 2.2 A simplified stability criterion

The criteria for long-term stability of multi-planet systems are not completely understood. However, N-body integrations of systems of several planets on nearly circular and coplanar orbits suggest that these systems are stable if the following approximate stability criterion is satisfied (Chambers et al., 1996; Yoshinaga et al., 1999; Zhou et al., 2007; Smith & Lissauer, 2009; Funk et al., 2010; Pu & Wu, 2015):

 ai+1−ai>Δcrit(t)rH;i,i+1whererH;i,i+1=ai+ai+12(mi+mi+13M⋆)1/3 (2)

is called the mutual Hill radius. This formula assumes that and . The stability parameter is a dimensionless function of the age of the system in units of the orbital period. Note that most of the simulations in the literature use equally spaced planets (either in semimajor axis or log semimajor axis) so they cannot determine whether (for example) the mean semimajor axis difference or the minimum semimajor axis difference determines stability. Our stability criterion (2) assumes that the minimum separation should be used in the stability criterion, as would be natural if stability is determined mostly by interactions with the nearest neighbor.

The most recent and comprehensive review of numerical determinations of is by Pu & Wu (2015). For typical observed exoplanetary systems, with , they estimate . However, their values are somewhat below those determined by other studies (see Figure 3 of Pu & Wu 2015) so we shall adopt a slightly more conservative estimate . For N-body simulations, which typically last for orbits, we use .

The form of equation (2) suggests that the minimum stable separation scales with mass as . This empirical finding is not a consequence of rigorous dynamical arguments. For example, the resonance overlap criterion (Wisdom, 1980; Deck et al., 2013) suggests the scaling ; however, the exponents and are sufficiently close that they are hard to distinguish in N-body experiments and the two scalings yield almost the same results for the purposes of this paper.

For eccentric orbits, the most important combination of orbital elements that determines the stability of adjacent planets is the distance between the apocenter of the inner orbit and the pericenter of the outer one (Pu & Wu, 2015; Petrovich, 2015). We therefore assume that the system is stable if

 ai+1−¯aei+1−ai−¯aei>hi (3)

where is some function of the masses and . This reduces to (2) in the limit of circular orbits if

 hi=Δcrit(t)¯a(mi+mi+13M⋆)1/3. (4)

An alternative criterion is given by Petrovich (2015), who has conducted extensive numerical experiments on the stability of systems of two planets with masses . The planets are on eccentric orbits and are followed for up to orbital periods. His empirical criterion for stability of closely spaced planets is equation (3) with

 hi=2.4¯a(max\,(mi,mi+1)M⋆)1/3+0.15¯a. (5)

We normally use equation (4) in our experiments, but we have also experimented with (5) and report results with both stability criteria.

The stability criterion (3) implies that between any pair of adjacent planets in a stable system there is an excluded length . In a system of planets there is a total excluded length (for consistency with the assumptions of §2.1 we take when determining and ). Obviously if this sum exceeds no stable planetary system exists. Thus an important parameter is the filling factor or packing fraction

 F=∑Ni=0hiΔa,0

### 2.3 The distribution of orbits in phase space

The phase space for coplanar Keplerian orbits has two degrees of freedom and actions , . In Hill’s approximation the actions can be taken to be and . The canonical phase-space volume element is .

We shall assume that the planets are uniformly distributed over the available phase space, that is, over the phase-space volume that satisfies the stability criterion (3). This assumption is reminiscent of the ergodic hypothesis, which is the basis of much of equilibrium statistical mechanics, but has a different interpretation in this case: it would apply, for example, if some process instantaneously and randomly re-distributed the planets throughout phase space in an ensemble of planetary systems and then only the systems having stable configurations survived. This assumption may well be incorrect: it is more likely that there is a slow diffusion of planets into the unstable region of phase space, which would reduce the phase-space density near the stability boundary. Nevertheless, the ergodic hypothesis is so simple and powerful that it is worthwhile to explore its implications before investigating more complex models with more free parameters. We shall call the models explored here “ergodic models”.

Given ergodicity, the -planet distribution function of semimajor axes and eccentricities is

 dp(a,e)=C(π2Ω2c¯a3)NH(a1−¯ae1−a0−h0)N∏i=1daide2iH(ai+1−¯aei+1−ai−¯aei−hi), (7)

where is the step function, is a normalizing constant, and we set .

### 2.4 Partition function

The partition function is the available volume in phase space and is given by the integral of the right side of equation (7) over the semimajor axes and eccentricities, with . To evaluate this integral we first carry out the integration over semimajor axes. For given eccentricities the total excluded length is ; the integral over depends only on this total and can be written as

 ∫Δa−G0dx1∫Δa−Gx1dx2⋯∫Δa−GxN−1dxN=1N!HN(Δa−G) (8)

where

 HN(x)={xNif x>00if x≤0. (9)

Note that is equal to the step function .

After doing the integral over the semimajor axes, the partition function becomes

 Z=(π2Ω2c¯a3)NN!N∏i=1∫de2iHN(Δa−∑Ni=0hi−2¯a∑Ni=1ei). (10)

In the sheared-sheet approximation, the integrals over van be taken from zero to infinity; then they can be done by induction, giving finally

 Z=π2NΩ2Nc¯aN(Δa)3N2N(3N)!(1−F)3N (11)

where the filling factor is defined in equation (6).

This result can be compared to the partition function for a one-dimensional gas of hard rods of length enclosed in a box of length (Tonks, 1936; Lieb & Mattis, 1966). In this case the filling factor is and the partition function is

 Z=LN(1−F)NN!. (12)

In both cases the partition function depends only on the filling factor. The main difference between the gas of hard rods and our model is that the exponent is replaced by , reflecting the extra degrees of freedom arising from the eccentricities.

### 2.5 The N-planet distribution function

The -planet distribution function (7) can be rewritten using the identity

 H(x)=limϵ→0+12πi∫∞−∞dss−iϵexp(isx) (13)

(for brevity, in future equations the limit is not written explicitly but is assumed to apply whenever appears). We have

 dp(a,e) =(3N)!¯a2N2(πi)N+1(Δa)3N(1−F)3N∫∞−∞ds0s0−iϵexp[is0(a1−¯ae1−a0−h0)] ×N∏i=1daide2i∫∞−∞dsisi−iϵexp[isi(ai+1−¯aei+1−ai−¯aei−hi)], (14)

in which we have replaced the normalization constant by so , and as usual , , .

The product of exponentials can be re-written as where

 Φ =s0(a1−¯ae1−a0−h0)+∑Ni=1si(ai+1−¯aei+1−ai−¯aei−hi) (15) =N∑i=1[ai(si−1−si)−¯aei(si−1+si)−sihi]+sN(¯a+12Δa)−s0(¯a−12Δa)−s0h0.

The characteristic function of the -planet distribution function is

 PN(k,p)=∫dp(a,e)exp[i∑Ni=1(kiai+piei)]. (16)

Carrying out the integrals over ,

 PN(k,p) =2N−1(3N)!¯a2NπiN+1(Δa)3N(1−F)3N∫∞−∞ds0s0−iϵexp[−is0(¯a−12Δa+h0)] ×N∏i=1∫∞−∞δ(ki+si−1−si)dsisi−iϵexp[−isihi+iδiNsN(¯a+12Δa)] ×∫∞0de2iexp[iei(pi−¯asi−1−¯asi)]. (17)

Next integrate over after multiplying by to ensure convergence:

 PN(k,p) =22NiN(3N)!¯a2N2πi(Δa)3N(1−F)3N∫∞−∞ds0s0−iϵexp[−is0(¯a−12Δa+h0)] ×N∏i=1∫∞−∞δ(ki+si−1−si)dsiexp[−isihi+iδiNsN(¯a+12Δa)](si−iϵ)(pi−¯asi−1−¯asi+iϵ)2. (18)

This result provides a complete description of the joint probability distribution of the semimajor axes and eccentricities of all planets.

### 2.6 The one- and two-planet distribution functions

For practical purposes, the one- or two-planet distributions are more useful than the full -planet distribution. Suppose we want to study the -planet distribution function, starting with planet . Then equation (18) implies that and this implies that . Moreover . Similarly and . Thus the -planet characteristic function is

 PK(kJ+1,⋯,kJ+K,pJ+1,⋯pJ+K) =22KiN(3N)!¯a2K2πi(Δa)3N(1−F)3N∫∞−∞dsJexp[isJ(12Δa−¯a−∑Jj=0hj)](sJ−iϵ)3J+1 (19) ×[J+K∏i=J+1∫∞−∞dsiexp(−isihi)δ(ki+si−1−si)(si−iϵ)(pi−¯asi−1−¯asi+iϵ)2]exp[isJ+K(¯a+12Δa−∑Nj=J+K+1hj)](sJ+K−iϵ)3(N−J−K).

For example let . Then

 P1(k,p) =2iN−1(3N)!¯a2π(Δa)3N(1−F)3Nexp[ik(¯a−12Δa+∑Jj=0hj)] ×∫∞−∞dsexp[isΔa(1−F)](s−iϵ)3(N−J)−2(s−k−iϵ)3J+1(p+k¯a−2s¯a+iϵ)2. (20)

with , , . The 1-planet distribution function is the inverse Fourier transform of the characteristic function,

 p1(a,e) =1(2π)2∫∞−∞dkdpexp[−i(ka+pe)]P1(k,p) =4(3N)!¯a2(3J)![3(N−J−1)]!(Δa)3N(1−F)3NH1(e)H3J[a−(¯a−12Δa+∑Jj=0hj+¯ae)] ×H3(N−J−1)[(¯a+12Δa−∑Nj=J+1hj−e¯a)−a]. (21)

If we integrate over semimajor axis, we have

 p1(e)=∫dap1(a,e)=12N(3N−1)¯a2(Δa)2(1−F−2e¯a/Δa)3N−2(1−F)3Ne,F<1. (22)

Notice that the result is independent of the planet number . Moreover it is independent of the excluded lengths associated with the planet in question or its neighbors; the distribution of eccentricities depends on these lengths only though the sum of the excluded lengths, which determines the filling factor through (6).

Similarly, we can evaluate the two-planet distribution function for adjacent planets:

 p2(aJ+1,eJ+1,aJ+2,eJ+2) =24(3N)!¯a4(3J)!(3N−3J−6)!(Δa)3N(1−F)3NH1(eJ+1)H1(eJ+2) ×H3J[aJ+1−¯aeJ+1−(¯a−12Δa+∑Jj=0hj)] ×H3N−3J−6[(¯a+12Δa−∑Nj=J+2hj)−aJ+2−¯aeJ+2] ×H0(aJ+2−¯aeJ+2−aJ+1−¯aeJ+1−hJ+1). (23)

We are interested in the dependence of the two-planet function on separation so we may integrate over keeping the separation fixed. We simplify the notation by setting , , , :

 p(a′−a,e,e′) =16(3N)!¯a4(3N−5)!(Δa)3N(1−F)3NH1(e)H1(e′)H0[a′−a−¯a(e+e′)−hJ+1] ×H3N−5[Δa(1−F)−(a′−a)−¯a(e+e′)+hJ+1]. (24)

Higher-order distribution functions are more complicated but could be useful; for example, the three-planet distribution function can be used to describe the distribution of semimajor axis differences if an intermediate planet is present but undetected.

## 3 Comparison to simulations and observations

We now re-cast the formulas of the preceding sections into simpler forms that can be more directly compared to numerical simulations, or to observations. We assume that the number of planets is large, so any single system should be considered as a subsystem of one with a much larger radial width . This approach is similar to the grand canonical ensemble of classical statistical mechanics. Formally we let , , while the filling factor . Thus, for example, a factor like .

The eccentricity distribution (eq. 22) is then

 p1(e)=eτ2exp(−eτ),τ≡Δa(1−F)6N¯a. (25)

The mean eccentricity is and we call the dynamical temperature since it parametrizes the level of non-circular motion in the planetary systems. The second of equations (25) shows that the dynamical temperature is related to the filling factor and the mean separation ; thus it is determined by the number and masses of the planets. Note also that the distribution (25) has a fatter tail at high eccentricities than the standard Rayleigh distribution, .

The two-planet distribution function (24) is

 p2(a′−a,e,e′) ×exp[−a′−a+¯a(e+e′)−h2¯aτ] (26)

where is the excluded length between the planets at and (eq. 4). Integrating over semimajor axes, the joint distribution in eccentricities is

 p2(e,e′)=ee′τ4exp(−e+e′τ). (27)

Thus the distribution of eccentricity of nearest neighbors is separable: the two-planet eccentricity distribution is the product of two one-planet distributions (eq. 25) and there is no correlation or anti-correlation between the eccentricities of adjacent planets.

The distribution in semimajor axis difference and total eccentricity is

 p2(a′−a,et) =∫∞0ede∫∞0e′de′δ(et−e−e′)p2(a′−a,e,e′) =e3t6¯aτ5H0(a′−a−¯aet−h)exp(−a′−a+¯aet−h2¯aτ). (28)

Integrating over the eccentricities gives

 p2(a′−a)=43¯aτD(a′−a−h2¯aτ) (29)

where . Taking means yields

 ⟨a′−a⟩=⟨h⟩+6¯aτ, (30)

From equation (6), for the filling factor can be estimated as

 F=⟨h⟩⟨a′−a⟩, (31)

and with this estimate equation (30) is equivalent to the second of equations (25).

Other authors have examined the relation between the semimajor axis differences of nearest neighbors and their mutual Hill radius, but typically by plotting the distribution of rather than (Fang & Margot, 2013; Hansen & Murray, 2013; Lissauer et al., 2014). In the ergodic model the second of these has a simpler interpretation: the distribution of is given by equation (29) so long as the dynamical temperature is similar for all systems in the sample, whereas deriving the distribution of from the ergodic model requires knowledge of the distribution of as well.

These formulas give the distribution of eccentricities and semimajor axis differences for a given value of the free parameter . Of course, different planetary systems may have different values of , so fitting the formulas to a large catalog of planets is only legitimate if the value of does not vary too much among the systems in the catalog. To avoid this difficulty, we may plot the normalized eccentricity,

 E≡¯a(e′+e)a′−a−h, (32)

which has the distribution

 p(E)=64E3(1+E)5,0≤E<1, (33)

independent of .

If the late stages of terrestrial planet formation are driven by long-term instabilities and collisions in all planetary systems, and all collisions result in perfect mergers (as assumed in most simulations), then we would expect the buildup of planets to proceed self-similarly, so the final filling factor would be similar in all planetary systems. If so, then at a given semimajor axis, age, and stellar mass (eq. 4), (eq. 31), (eq. 30), and the surface density where is the typical planet mass.

These scalings fail when the impact velocities become sufficiently large that the collisions are erosive (Schlichting, 2014). To explore the effect of this failure, we parametrize the collisional environment by the Safronov number,

 Θ=GmR⟨v2⟩, (34)

where and are the planetary mass and radius. Here is the mean-square velocity relative to the local circular speed, equal to where and are the mean-square eccentricity and inclination. Equation (25) gives so in the simple case where ,

 Θ=4ma27M⋆Rτ2. (35)

For the sample of Kepler planets analyzed in §3.2, so

 ⟨Θ⟩=0.0037τ2=(0.06τ)2. (36)

The typical collision becomes erosive when –3, where is the escape speed from the larger body and is the rms relative velocity at infinity (see Leinhardt & Stewart 2012 for a much more thorough analysis). Thus we expect that the giant-impact phase of planet formation should have –3, which in turn implies an upper limit –0.03 for the sample of Kepler planets. These arguments are roughly consistent with the results of Chambers (2013), who conducted N-body simulations without and with fragmentation in collisions and found mean eccentricities and 0.045 respectively, corresponding to and 0.02.

As discussed in the Introduction, the approach in this paper is to model the distribution of orbital elements of a set of planets of given masses, but not to model how the distribution of planetary masses is established. This approach is somewhat artificial since collisions, merging or erosive, affect both the mass and orbit distributions simultaneously. Nevertheless, the conclusion that the distribution of planetary masses affects the orbital distribution only through the dynamical temperature or filling factor (the two being related by eq. 25) appears to be a plausible and inevitable consequence of the ergodic model.

### 3.1 Comparison to simulations

We compare the predictions of the ergodic model to simulations of the late stages of terrestrial planet formation carried out by Hansen & Murray (2013)222Many other authors have also simulated the giant-impact stage of planet formation (Chambers & Wetherill, 1998; Agnor et al., 1999; Chambers, 2001; Kokubo et al., 2006; Raymond et al., 2006; Morishima et al., 2008; Raymond et al., 2009; Ida & Lin, 2010; Kokubo & Genda, 2010; Morishima et al., 2010; Fischer & Ciesla, 2014; Pfyffer et al., 2015) but the Hansen & Murray simulation offers the most direct comparison to our model.. They divide of material, distributed between and , into planets and follow the evolution of these planets for 10 Myr. Collisions are assumed to result in mergers. The calculation is repeated 100 times with randomly varying initial conditions to build up the statistics. Following Hansen & Murray, we do not include planets with in our analyses since these are mostly scattered objects that are not expected to fit the ergodic model.

Fitting equation (25) to the eccentricity distribution of 527 surviving planets in the simulations we obtain (1– error, or reduction in log likelihood of ). With this value of the theoretical eccentricity distribution (25) is a very good match to the distribution found in the simulations, as shown in Figure 1. Indeed, Hansen & Murray (2013) proposed the eccentricity distribution (25) as an empirical fitting formula for their results.

Equation (30) relates the mean of the relative semimajor axis differences of nearest-neighbor pairs to . In this case the relation depends on the stability criterion, i.e., on the excluded length . Using the stability criterion (3) with we obtain for the 426 pairs in the Hansen & Murray simulation; the alternative stability criterion (5) yields . These values for the dynamical temperature are roughly consistent with the value determined from the eccentricities, which is an encouraging test of the consistency of the ergodic model. Moreover the distribution of semimajor axis differences is fit well by the predicted distribution (29), as shown in Figure 2. The filling factor, as determined by equation (31), is using the stability criterion (3) and 0.36 using (5).

The distribution of normalized eccentricity (eq. 32) in the simulations is shown in Figure 3, along with the prediction of the ergodic model. In this case the predicted distribution does not match the simulations well. Part of the discrepancy is that 11% of the planet pairs have , and these would be unstable according to equation 4. The other major discrepancy is that there are fewer planets in the simulation with –0.7 than the model predicts. These probably reflect two oversimplifications of the ergodic model: (i) the stability criterion is not a sharp boundary, as assumed in the derivation; (ii) planets diffuse in eccentricity towards the stability boundary, so we expect their phase-space density to be lower near the boundary that the uniform density predicted by the ergodic model.

Figure 4 shows contour plots of the joint probability density for the reduced semimajor axis difference and total eccentricity of adjacent planets, derived from equation (28). The dynamical temperature is assumed to be . The red circles show the nearest-neighbor pairs in the Hansen & Murray (2013) simulation. In this plot, contours of constant are straight lines through the origin; the dashed line is . The deficit of planets at –0.7 seen in Figure 3 is evident, since the density of red circles does not rise as fast as the probability density near the line. This discrepancy is less visible in plots of the eccentricity or semimajor axis distribution because these are projections onto the vertical or horizontal axes.

### 3.2 Comparison to data

The final and most important step is to compare our predictions to the actual properties of exoplanets. There are several obstacles to accomplishing this task: (i) Generally, eccentricities are only available for planets whose orbits have been measured from radial velocities. (ii) The masses of most planets discovered by Kepler can only be estimated using an empirical mass-radius relation, and individual planet masses exhibit large deviations from the mean relation (Weiss & Marcy, 2014). (iii) There may be undiscovered planets in between the known members of a multi-planet system and these would affect the distribution of semimajor axis differences.

#### Kepler planets:

We have queried the NASA Exoplanet Archive for all systems containing more than one confirmed planet discovered by Kepler. This sample provides 932 planets in 362 distinct systems, containing 556 nearest-neighbor pairs. We estimate the planetary masses from the radii using the mass-radius relation from Weiss & Marcy (2014). We then compute the relative semimajor axis difference for each pair, , using . Using the estimated masses and the stability criterion (4) with we also compute the exclusion length and thus determine the distribution of , which ought to be described by equation (29).

The result is shown as the histogram in Figure 5. The general shape of the histogram is similar to that of the histogram in Figure 2 from the Hansen & Murray (2013) simulations. One obvious difference, however, is that the histogram derived from the Kepler data has a significant number of nearest neighbor pairs with negative values of (48 out of 556). These pairs should be unstable so their presence must be explained. One possibility is that our stability criterion is too conservative, but this is unlikely since the analogous distribution from the Hansen & Murray (2013) simulations contains no pairs with . Using the mass-radius relation from Fabrycky et al. (2014) or the alternative stability criterion (5) does not remove this difficulty: between 31 and 63 planet pairs in the sample are still apparently unstable. The most plausible explanation is that the unstable pairs arise because of scatter in the mass-radius relation.

To explore this possibility we shall assume that the distribution of errors in the excluded length (eq. 4) is Gaussian, with standard deviation . Since for fixed stellar mass, should be related to the error in , roughly as . To estimate the latter quantity, we use the sample of 65 exoplanets with measured masses and radii compiled by Weiss & Marcy (2014), and compare the observed values of with the predictions from the mean mass-radius relation derived in that paper. (The mean mass of the Weiss & Marcy sample, , is close to the mean for the planets in our sample, , suggesting that the two samples have similar properties.) We find . The mean excluded length in our sample is so we estimate from these arguments.

We may now convolve our predicted distribution of semimajor axis differences (29) with a Gaussian of width to find the distribution that would be obtained given realistic errors in the planet masses. We fit the resulting distribution to the data in Figure 5 using the measurement error and as fitting parameters, and find a best fit with , . The best-fit value for is remarkably close to the value obtained from the analysis of the Weiss & Marcy (2014) sample, ; thus the distribution of semimajor axis differences in the Kepler sample appears to be consistent with the ergodic theory once the scatter in the mass-radius relation is taken into account.

The fit to the data has per degree of freedom is 2.7, with most of the contribution to coming from a tail of planets with that is not present in the theoretical curve. The tail could arise from systems in which an intermediate planet was below the detection limit of the Kepler survey.

The filling factor is . Using the mass-radius relation from Fabrycky et al. (2014) instead of Weiss & Marcy (2014) changes this by less than the error bar, to . Using the stability criterion (5) reduces the filling factor somewhat, to .

Given the measured value for the Kepler semimajor axis differences, the mean eccentricity should be –0.06 (eq. 25), unless the eccentricities have been damped after the giant-impact phase is complete (Hansen & Murray, 2015). So far, we have only limited information on the distribution of eccentricities of the Kepler planets. (i) Hadden & Lithwick (2014) have estimated the eccentricities using transit timing variations333There is a typographical error in the abstract of this paper. The quantity is not the rms eccentricity; it is which equals the rms eccentricity divided by . and obtain ; for planets larger than the mean eccentricity is a factor of two smaller, or about . (ii) Moorhead et al. (2011) have estimated eccentricities from the distribution of transit durations and obtain –0.25. The larger value relative to Hadden & Lithwick may arise in part because transit timing variations can only be measured in multi-planet systems and such systems are expected to have smaller eccentricities (see below); or because planets with large transit timing variations are mostly near strong resonances, and such planets could have a different eccentricity distribution. (iii) Several authors have estimated the mean inclination of the Kepler planets. Tremaine & Dong (2012) find ; in most astrophysical disks, –2 times so this result implies . Similarly, Fang & Margot (2012) and Figueira et al. (2012) find , so ; Johansen et al. (2012) find , so ; and Fabrycky et al. (2014) find and so . A concern with all of these comparisons is that the multi-planet systems observed by Kepler may be biased towards low eccentricities because eccentricity and inclination are correlated and low-inclination systems are more likely to have multiple transits. We conclude that the observations so far are roughly consistent with the estimate of the mean eccentricity from the ergodic model, –0.06, but do not provide strong support for it.

The filling factor estimated from the Hansen & Murray (2013) simulations is , while the filling factor estimated from the semimajor axis differences of Kepler planets is 0.4–0.5 depending on the mass-radius relation and stability criterion. If the dynamics of the late stages of formation of the Kepler planets are faithfully modeled by the simulations, we might expect the two filling factors to be the same. The similarity of the two numbers is impressive but it is worthwhile to ask why they might differ. One intriguing possibility is that the Kepler planet masses have been overestimated; decreasing the masses by a factor of two would decrease the Kepler filling factor by about 0.1 and bring it to agreement with the filling factor in the simulations.

The values of the dynamical temperature are also similar: 0.06 in the Hansen & Murray (2013) simulations and 0.03 in the Kepler data. This difference cannot be accounted for by differences in the surface density or mass: the scalings described after equation (33) imply and , which would predict that should be 15% smaller in the simulations than in the data. A more likely cause of the difference is that the simulations did not allow for fragmentation: the arguments following equation (36) imply that if fragmentation is present the dynamical temperature cannot exceed 0.02–0.03, consistent with the Kepler data.

We have queried the Exoplanets Data Explorer (Han et al., 2014) for multiple-planet systems discovered by radial-velocity variations in their host stars. This sample provides 135 planets in 55 systems or 77 nearest neighbor pairs. The distribution of eccentricities is fit well by equation (25) with ; this value exceeds the upper limit derived after equation (36) by a factor of three or more but this is not a contradiction because the giant planets probably did not grow to their present masses by giant impacts. The distribution of semimajor axis differences, however, is not well-fit by the ergodic model. There are several likely reasons for this. (i) Almost one-third of the planet pairs violate the stability criterion (3) (we estimated the masses assuming the systems are edge-on, which gives a lower limit); this is partly because the criterion is not valid for planets in mean-motion resonances, which are common among giant planets. (ii) The ergodic model is based on the sheared-sheet approximation, which requires that the semimajor axis difference ; this is generally not true for planets larger than a Jupiter mass since the excluded length (eq. 4) satisfies . (iii) Giant planets formed by processes that differ from those of terrestrial planets, for example migration and the accretion of massive gas envelopes. As described at the end of §1, the ergodic model is only expected to apply to giant planets if their orbits are mostly determined by late-stage dynamical evolution.

Equation (25) suggests that in an ensemble of systems with the same filling factor the mean eccentricity should vary as (more precisely, inversely with the number of planets per unit radius). Limbach & Turner (2015) find that for known exoplanet systems the eccentricity decreases with multiplicity roughly as , in reasonable agreement with this prediction.

#### The solar system:

Figure 6 shows contour plots of the joint probability density for the reduced semimajor axis difference and total eccentricity of adjacent planets, derived from equation (28). The dynamical temperature is assumed to be in the left panel, the same as we found for the Kepler planets, and 0.03 on the right. The red circles show the actual values for the three nearest-neighbor pairs among the four terrestrial planets. For most of the weight of the probability density distribution lies at total eccentricities larger than those of the solar-system planet pairs. For the Mercury-Venus pair has too large a semimajor axis difference relative to the probability density distribution. These mismatches presumably reflect the well-known difficulty that simulations of the giant-impact phase produce eccentricities for the terrestrial planets of the solar system that are too large, unless there is a residual population of small planetesimals to damp the eccentricities (e.g., Morbidelli et al., 2012).

The outer planets of the solar system do not fit the ergodic model well, in part because Jupiter and Saturn are separated by only eight mutual Hill radii and so would nominally be unstable according to our crude stability criterion.

## 4 Discussion

We have described a simple model for the distribution of semimajor axes and eccentricities of planets. The model assumes that the last phase of terrestrial planet formation was the giant-impact phase, and is based on the simple ansatz that planets are uniformly distributed over the volume of phase space in which their orbits are stable for the lifetime of the planetary system (the “ergodic model”). Our model yields predictions, in terms of a single free parameter that we call the dynamical temperature, for the distribution of eccentricities (eq. 25), the distribution of separations of nearest neighbors (eq. 29), and more generally for any property derivable from the complete -planet distribution function. For example, the ergodic model predicts that in a given system the eccentricities should be independent of planetary mass. N-body simulations generally report only a weak negative correlation between eccentricity and mass (e.g., Chambers, 2013).

The ergodic model has many limitations. (i) It does not account for the likely influence of giant planets at larger radii on the formation of terrestrial planets. (ii) Different planetary systems may have different dynamical temperatures and fitting the data from a large ensemble of systems to the ergodic model will only work well if most systems have similar temperatures. (iii) Our analysis is based on the sheared-sheet or Hill’s approximation and hence does not work well when the planetary masses are large enough that the relative separations or the fractional excluded lengths (eq. 4) are of order unity or larger; typically this occurs for planets more massive than Jupiter. (iv) The stability criteria (4) and (5) are only approximate and probably no simple criterion perfectly separates stable from unstable orbits in a multi-planet system (Petrovich, 2015). (v) Our assumption that the orbits are distributed uniformly over the stable part of phase space is shaky, essentially because the neighboring unstable regions represent an absorbing barrier rather than a reflecting one. A better approximation would require following the diffusion of planet orbits through phase space, but implementing this would require a reliable model for the rate of this diffusion and how it depends on the orbits of the planet and its neighbors.

In this paper we examine only the distribution of orbital elements and not the distribution of planetary masses that emerges in the giant-impact phase, and any complete statistical model of this phase should predict both. An interesting question is whether the system described in §2.1, based on the sheared-sheet or Hill’s approximation, exhibits long-range order in the masses as ; in other words if the radial width but the average surface density is fixed, does the mass of the most massive planet grow as and if so what is the critical exponent ?

Although the ergodic model is simple and physically plausible, there are other approaches to the statistical mechanics of planet formation. An interesting alternative is due to Laskar (2000), who pointed out that in the secular approximation there is an important conserved quantity: the angular-momentum deficit,

 C=N∑i=1(GM⋆ai)1/2[1−(1−e2i)1/2cosIi], (37)

that is, the difference between the total angular momentum of the planets and the angular momentum that they would have on circular, coplanar orbits with the same semimajor axes. He hypothesized that the eccentricities and inclinations of planets evolve randomly, subject to conservation of the angular-momentum deficit, until there is a binary collision; in each collision the angular-momentum deficit is reduced; and collisions cease and the system becomes permanently stable once the angular-momentum deficit becomes too small to allow any more close encounters. Laskar’s model is not unique, since it requires an ad hoc prescription for the “random” evolution of the orbits; however, once this prescription is implemented it is straightforward to predict both the masses and the orbital elements of the planets in an ensemble of systems. Laskar’s model differs from ours in that it generally predicts a strong anti-correlation between eccentricity and mass.

One-dimensional models are powerful tools in statistical mechanics because they are often much simpler to solve than their three-dimensional analogs (Lieb & Mattis, 1966). It is remarkable that the giant-impact phase of planetary formation is most naturally modeled in one physical dimension, radius, basically because (i) the systems are nearly flat; (ii) orbital and apsidal motion effectively averages the orbital and collisional dynamics over the azimuthal angle. Although the ergodic model contains many simplifications, more general and accurate analyses of the statistical mechanics of the giant-impact phase could still be one-dimensional in this sense and therefore amenable to analytic treatments.

Exact models in statistical mechanics are mainly useful because they provide insight into the behavior of real systems, rather than because they are accurate representations of them. Similarly, the ergodic model presented here is mainly useful because it provides simple predictions for many properties of the 1–, 2–, and even –planet joint distribution of orbital elements. These predictions encapsulate some, though certainly not all, of the physics of the late stages of planet formation, and therefore should help to organize and interpret both observations of planet orbits and numerical simulations of planet formation.

This research was initially stimulated by conversations with Renu Malhotra. I thank Brad Hansen for comments on the manuscript, and for providing the results from his simulations. I am particularly grateful to Cristobal Petrovich for discussions, insight, and pointers to the literature, and for carrying out exploratory N-body integrations on my behalf.

## References

• Agnor et al. (1999) Agnor, C. B., Canup, R. M., & Levison, H. F. 1999, Icarus, 142, 219
• Barnes & Raymond (2004) Barnes, R., & Raymond, S. N. 2004, ApJ, 617, 569
• Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics (Princeton: Princeton University Press)
• Chambers (2001) Chambers, J. E. 2001, Icarus, 152, 205
• Chambers (2013) Chambers, J. E. 2013, Icarus, 224, 43
• Chambers & Wetherill (1998) Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304
• Chambers et al. (1996) Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261
• Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580
• Deck et al. (2013) Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, id. 129
• Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, id. 146
• Fang & Margot (2012) Fang, J., & Margot, J.-L. 2012, ApJ, 761, id. 92
• Fang & Margot (2013) Fang, J., & Margot, J.-L. 2013, ApJ, 767, id. 115
• Figueira et al. (2012) Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, id. A139
• Fischer & Ciesla (2014) Fischer, R. A., & Ciesla, F. J. 2014, Earth and Planetary Science Letters, 392, 28
• Funk et al. (2010) Funk, B., Wuchterl, G., Schwarz, R., Pilat-Lohinger, E., & Eggl, S. 2010, A&A, 516, id. A82
• Hadden & Lithwick (2014) Hadden, S., & Lithwick, Y. 2014, ApJ, 787, id. 80
• Haghighipour (2013) Haghighipour, N. 2013, Annual Review of Earth and Planetary Sciences, 41, 469
• Han et al. (2014) Han, E., Wang, S. X., Wright, J. T., et al. 2014, PASP, 126, 827
• Hansen & Murray (2013) Hansen, B.M.S., & Murray, N. 2013, ApJ, 775, id. 53
• Hansen & Murray (2015) Hansen, B.M.S., & Murray, N. 2015, MNRAS, 448, 1044
• Hénon (1969) Hénon, M. 1969, A&A, 1, 223
• Hénon (1970) Hénon, M. 1970, A&A, 9, 24
• Holman (1997) Holman, M. J. 1997, Nature, 387, 785
• Holman & Wisdom (1993) Holman, M. J., & Wisdom, J. 1993, AJ, 105, 1987
• Ida & Lin (2010) Ida, S., & Lin, D. N. C. 2010, ApJ, 719, 810
• Johansen et al. (2012) Johansen, A., Davies, M. B., Church, R. P., & Holmelin, V. 2012, ApJ, 758, id. 39
• Jurić & Tremaine (2008) Jurić, M., & Tremaine, S. 2008, ApJ, 686, 603
• Kokubo & Genda (2010) Kokubo, E., & Genda, H. 2010, ApJ, 714, L21
• Kokubo & Ida (2012) Kokubo, E., & Ida, S. 2012, Progress of Theoretical and Experimental Physics, 2012, id. 01A308
• Kokubo et al. (2006) Kokubo, E., Kominami, J., & Ida, S. 2006, ApJ, 642, 1131
• Laskar (1996) Laskar, J. 1996, in Dynamics, Ephemerides, and Astrometry of the Solar System, IAU Symposium 172, ed. S. Ferraz–Mello, B. Morando, and J.-E. Arlot (Dordrecht, NL: Kluwer), 75
• Laskar & Gastineau (2009) Laskar, J., & Gastineau, M. 2009, Nature, 459, 817
• Leinhardt & Stewart (2012) Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, id. 79
• Lieb & Mattis (1966) Lieb, E. H., & Mattis, D. C. 1966, Mathematical Physics in One Dimension (New York: Academic Press)
• Limbach & Turner (2015) Limbach, M. A., & Turner, E. L. 2015, PNAS, 112, 20
• Lissauer et al. (2014) Lissauer, J. J., Dawson, R. I., & Tremaine, S. 2014, Nature, 513, 336
• Malhotra (2015) Malhotra, R. 2015, arXiv:1502.05011
• Moorhead et al. (2011) Moorhead, A. V., Ford, E. B., Morehead, R. C., et al. 2011, ApJS, 197, id. 1
• Morbidelli et al. (2012) Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Annual Review of Earth and Planetary Sciences, 40, 251
• Morishima et al. (2008) Morishima, R., Schmidt, M. W., Stadel, J., & Moore, B. 2008, ApJ, 685, 1247
• Morishima et al. (2010) Morishima, R., Stadel, J., & Moore, B. 2010, Icarus, 207, 517
• Petit & Hénon (1986) Petit, J.-M., & Hénon, M. 1986, Icarus, 66, 536
• Petrovich (2015) Petrovich, C. 2015, in preparation
• Pfyffer et al. (2015) Pfyffer, S., Alibert, Y., Benz, W., & Swoboda, D. 2015, arXiv:1502.04260
• Pu & Wu (2015) Pu, B., & Wu, Y. 2015, arXiv:1502.05449
• Raymond et al. (2006) Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265
• Raymond et al. (2009) Raymond, S. N., O’Brien, D. P., Morbidelli, A., & Kaib, N. A. 2009, Icarus, 203, 644
• Raymond et al. (2013) Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., & Walsh, K. J. 2014, in Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, T. Henning (Tucson: University of Arizona Press). Also arXiv:1312.1689.
• Schlichting (2014) Schlichting, H. E. 2014, ApJ, 795, id. L15
• Smith & Lissauer (2009) Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381
• Spitzer & Schwarzschild (1953) Spitzer, L., Jr., & Schwarzschild, M. 1953, ApJ, 118, 106
• Sussman & Wisdom (1988) Sussman, G. J., & Wisdom, J. 1988, Science, 241, 433
• Sussman & Wisdom (1992) Sussman, G. J., & Wisdom, J. 1992, Science, 257, 56
• Tonks (1936) Tonks, L. 1966, Phys. Rev., 50, 955
• Tremaine & Dong (2012) Tremaine, S., & Dong, S. 2012, AJ, 143, id. 94
• Volk & Gladman (2015) Volk, K., & Gladman, B. 2015, arXiv:1502.06558
• Weiss & Marcy (2014) Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, id. L6
• Wisdom (1980) Wisdom, J. 1980, AJ, 85, 1122
• Yoshinaga et al. (1999) Yoshinaga, K., Kokubo, E., & Makino, J. 1999, Icarus, 139, 328
• Zhou et al. (2007) Zhou, J.-L., Lin, D. N. C., & Sun, Y.-S. 2007, ApJ, 666, 423
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters