Pairwise Tidal Equilibrium States and the Architecture of Extrasolar Planetary Systems

TL;DR

This paper investigates the tidal equilibrium states of multi-planet systems, proposing that planetary formation tends to produce configurations with nearly equal masses and circular, coplanar orbits, aligning with observed system architectures.

Contribution

It introduces a model for pairwise tidal equilibrium states in multi-planet systems, explaining their observed regular spacing and mass distribution as a minimum energy configuration.

Findings

Observed systems are close to pairwise equilibrium states.

Systems do not reside in a global minimum energy state.

Surface density profile of such systems scales as r^{-2}.

Abstract

Current observations indicate that the planet formation process often produces multiple planet systems with nearly circular orbits, regular spacing, a narrow range of inclination angles, and similar planetary masses of order $m_{\rm p}\sim10M_\oplus$. Motivated by the observational sample, this paper determines the tidal equilibrium states for this class of extrasolar planetary systems. We start by considering two planet systems with fixed orbital spacing and variable mass ratios. The basic conjecture explored in this paper is that the planet formation process will act to distribute planetary masses in order to achieve a minimum energy state. The resulting minimum energy configuration --- subject to the constraint of constant angular momentum --- corresponds to circular orbits confined to a plane, with nearly equal planetary masses (as observed). We then generalize the treatment to include multiple planet systems, where each adjacent pair of planets attains its (local) tidal equilibrium state. The properties of observed planetary systems are close to those expected from this pairwise equilibrium configuration. In contrast, observed systems do not reside in a global minimum energy state. Both the equilibrium states of this paper and observed multi-planet systems, with planets of nearly equal mass on regularly spaced orbits, have an effective surface density of the form $\sigma\propto r^{-2}$, much steeper than most disk models.

Equations

$$\mathbf{L}=h_1\hat{z}+h_2\left[\hat{z}\cos i+\hat{x}\sin i\right],$$

$$h_j^2=m_j^{2}G{M}_{\ast }{a}_j(1-e_j^2),$$

$$E=-\frac{G{M}_{\ast }{m}_1}{2a_1}-\frac{G{M}_{\ast }{m}_2}{2a_2}.$$

$$m_1+m_2=m_T=constant,$$

$$a_2=a_1(1+\Delta )=a_1\Lambda ,$$

$$\Delta =K{R}_H/a_1\mathrm{and}{R}_H={\left(\frac{m_T}{3M_{\ast }}\right)}^{1/3}{a}_1.$$

$$f=\frac{m_1}{m_T},1-f=\frac{m_2}{m_T},a=\frac{a_1}{R_{\ast }},a\Lambda =\frac{a_2}{R_{\ast }}.$$

$$E=-\frac{1}{a}\left[f+\frac{1-f}{\Lambda }\right]$$

$$\mathbf{L}=\sqrt{a}\left\{f(1-e_1^2{)}^{1/2}\hat{z}+(1-f)\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\left[\hat{z}\cos i+\hat{x}\sin i\right]\right\}.$$

$$\frac{\partial E}{\partial {\xi }_k}+{\lambda }_z\frac{\partial L_z}{\partial {\xi }_k}+{\lambda }_x\frac{\partial L_x}{\partial {\xi }_k}=0,$$

$$\frac{1}{a^2}\left[f+\frac{1-f}{\Lambda }\right]+\frac{{\lambda }_z}{2a^{1/2}}\left[f(1-e_1^2{)}^{1/2}+(1-f)\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\cos i\right]+\frac{{\lambda }_x}{2a^{1/2}}(1-f)\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\sin i=0.$$

$$-\frac{1}{a}\left[1-\frac{1}{\Lambda }\right]+{\lambda }_z\sqrt{a}\left[(1-e_1^2{)}^{1/2}-\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\cos i\right]-{\lambda }_x\sqrt{a}\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\sin i=0.$$

$$-{\lambda }_z\sqrt{a}f\frac{e_1}{(1-e_1^2{)}^{1/2}}=0.$$

$$-\sqrt{a}(1-f)\sqrt{\Lambda }\left[{\lambda }_z\cos i+{\lambda }_x\sin i\right]\frac{e_2}{(1-e_2^2{)}^{1/2}}=0.$$

$$\sqrt{a}(1-f)\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\left[-{\lambda }_z\sin i+{\lambda }_x\cos i\right]=0.$$

$$-{\lambda }_z{a}^{3/2}=2\left[f+\frac{1-f}{\Lambda }\right]{\left[f+(1-f)\sqrt{\Lambda }\right]}^{-1},$$

$$-{\lambda }_z{a}^{3/2}=\left[1-\frac{1}{\Lambda }\right]{\left[\sqrt{\Lambda }-1\right]}^{-1}.$$

$$f=\frac{\Lambda +\sqrt{\Lambda }-2}{3(\Lambda -1)}.$$

$$f\to \frac{1}{2}-\frac{\Delta }{24}+\frac{{\Delta }^2}{48}+\dots ,$$

$$\frac{m_1}{m_2}=\frac{f}{1-f}=\frac{\sqrt{\Lambda }+2}{2\sqrt{\Lambda }+1}.$$

$$L=\sqrt{a}\frac{2}{3}\frac{\Lambda +\sqrt{\Lambda }+1}{\sqrt{\Lambda }+1}.$$

$$E=-\frac{4}{9L^2}{\left[\frac{\Lambda +\sqrt{\Lambda }+1}{\sqrt{\Lambda }+1}\right]}^2\left[\frac{\Lambda +\sqrt{\Lambda }+1}{3\Lambda }\right].$$

$$\frac{\Delta E}{E}=\frac{4(\Lambda +\sqrt{\Lambda }+1{)}^3}{27\Lambda (\sqrt{\Lambda }+1{)}^2}-1=\frac{(\sqrt{\Lambda }-1{)}^2(2+\sqrt{\Lambda }{)}^2(1+2\sqrt{\Lambda }{)}^2}{27\Lambda (\sqrt{\Lambda }+1{)}^2}\to \frac{3{\Delta }^2}{16}+\dots .$$

$${\mathbb{H}}_{jk}=\frac{{\partial }^{2}E}{\partial {\xi }_j\partial {\xi }_k}.$$

$${\ell }_1\equiv f(1-e_1^2{)}^{1/2}\mathrm{and}{\ell }_2\equiv (1-f)\sqrt{\Lambda }(1-e_2^2{)}^{1/2}.$$

$$L^2=a\left[{\ell }_1^2+{\ell }_2^2+2{\ell }_1{\ell }_2\cos i\right],$$

$$E=-\frac{1}{L^2}\left[{\ell }_1^2+{\ell }_2^2+2{\ell }_1{\ell }_2\cos i\right]\left[f+\frac{1-f}{\Lambda }\right].$$

$$\frac{\partial E}{\partial f}=-2\left[({\ell }_1+{\ell }_2\cos i)(1-e_1^2{)}^{1/2}-({\ell }_2+{\ell }_1\cos i)\sqrt{\Lambda }(1-e_2^2{)}^{1/2}\right]\left[f+\frac{1-f}{\Lambda }\right]-\left[{\ell }_1^2+{\ell }_2^2+2{\ell }_1{\ell }_2\cos i\right]\left[1-\frac{1}{\Lambda }\right].$$

$${\partial^{2}E\over\partial f^{2}}\Bigg{|}_{0}=-2(1-\sqrt{\Lambda})^{2}\left[f+{1-f\over\Lambda}\right]-4(1-\sqrt{\Lambda})\left[f+(1-f)\sqrt{\Lambda}\right]\left[1-{1\over\Lambda}\right]\,.$$

$${\partial^{2}E\over\partial f^{2}}\Bigg{|}_{0}=2(\sqrt{\Lambda}-1)^{2}{\Lambda+\sqrt{\Lambda}+1\over\Lambda}\,,$$

Full text

††pagerange: Pairwise Tidal Equilibrium States and the Architecture of Extrasolar Planetary Systems–Pairwise Tidal Equilibrium States and the Architecture of Extrasolar Planetary Systems††pubyear: 2019

1 Introduction

With thousands of extra-solar planets now discovered, the observed collection of planetary systems displays an enormous degree of diversity. On the other hand, the multiple-planet systems detected through the Kepler mission (Borucki et al., 2010; Batalha et al., 2011) exhibit several instances of uniformity: The typical planetary mass in this sample is of order 10 $M_{\oplus}$ (Zhu et al., 2018), and planets within the same system tend to have similar masses (Millholland et al., 2017; Wang, 2017) and sizes (Weiss et al., 2018a, b). The planets are found to orbit in nearly the same plane (Fang & Margot, 2012; Tremaine & Dong, 2012) and are remarkably stable against oscillations in their orbital inclination angles (so that they stay aligned; Becker & Adams 2016, 2017). The orbital eccentricities are low, much smaller than those of the exoplanet sample taken as a whole (Van Eylen & Albrecht, 2015; Hadden & Lithwick, 2017). The angular momentum vectors of the stellar spins and planetary orbits are well-aligned in most cases (Winn & Fabrycky, 2015). These planetary systems generally have regularly spaced orbits with relatively small separations (Rowe et al., 2014), but are spread out enough to remain stable (Pu & Wu, 2015). Although the observed orbital proximity often results in period ratios that are close to small integer values, the systems are rarely found in bonafide low-order mean motion resonances (Fabrycky et al., 2014). These features of the observed planetary systems suggests that some organizing principle is operating during the process of planet formation, i.e., during the phase when planetary masses and orbits are determined. The goal of this work is to show how considerations of energy optimization — subject to constraints of mass and angular momentum conservation — provide a partial explanation for the regularities observed in the sample of extrasolar multi-planet systems.

A classic dynamical issue is to find the tidal equilibrium state for astrophysical systems containing both self-gravity and angular momentum (beginning with Darwin 1879, 1880). Given its origins, this type of constrained optimization procedure is sometimes known as a Darwin Problem. Any type of dissipation causes physical systems to evolve toward lower energy states, so that the lowest energy (equilibrium) state is the preferred configuration. These states are often called tidal equilibrium states because tidal interactions cause dissipation in many astrophysical systems of interest. For example, in binary star systems with both rotational and spin angular momentum, the existence of an equilibrium state requires both a minimum total angular momentum and a sufficiently large orbital component (Counselman, 1973; Hut, 1980), and tidal interactions provide a mechanism to reach equilibrium (Hut, 1981). These results for two-body systems have been applied to Hot Jupiters (Levrard et al., 2009), and have been generalized to include the quadrupole moment of the central star (Adams & Bloch, 2015) and a third body in hierarchical star-planet-moon systems (Adams & Bloch, 2016).

In the aforementioned applications, the components of the system include spin angular momentum, orbital angular momentum, and the orbital energy; the masses of the bodies are considered fixed. Moreover, the standard approach only allows for the inclusion of one orbit — that of the binary system. In the present case, however, the systems of interest include several planets, with orbits that are generally wide enough to be decoupled from the stellar rotation over the lifetime of the system. The spin angular momenta of the planets themselves are negligible. On the other hand, one quantity of interest is the mass ratio for adjacent pairs of planets. We thus need to solve for the optimal distribution of planetary masses. Note that the mass fractions will be determined through the process of planet formation, whereas the goal of this paper is to find the preferred (minimal energy) configuration.

This work poses a new type of Darwin problem that includes multiple orbits and allows the masses of the bodies to vary (subject to the constraint that the total mass $m_{\scriptstyle T}$ in planets is held constant). All of the angular momentum is contained in the planetary orbits and the total ${\bf L}$ is constrained to be constant. In this formulation, the orbital spacing, denoted as $\Lambda$, is also held fixed. For two planet systems, the relevant variables thus include the planet masses $m_{j}$, the orbital eccentricities $e_{j}$, the relative inclination angle $i$ of the orbits, and the semimajor axis $a$ of the inner planet (note that the semi-major axis of the outer planet is set by the fixed orbital spacing $\Lambda$). The problem is thus to optimize the system properties $(m_{1},m_{2},e_{1},e_{2},i,a)$ subject to the constraints of constant $({\bf L},m_{\scriptstyle T},\Lambda)$. With the two planet solution in hand, we then generalize to systems with more planets. If multiple planet systems evolve toward configurations where each adjacent pair of planets resides in a tidal equilibrium state, the resulting systems closely resemble those that are observed (with nearly equal mass planets and regularly spaced circular orbits in the same plane).

The goal of this paper is to provide a partial explanation for the properties of observed multi-planet systems by showing that they reside near optimized energy configurations. This approach does not describe the evolutionary path by which they attain such states. Instead, this calculation highlights a key physical principle — namely energy minimization — that is likely to constrain any scenario for planetary assembly. As a result, this work is complementary to previous studies of planet population synthesis (e.g., Mordasini et al. 2009; Ida & Lin 2010; Alessi et al. 2017; Mordasini 2018) that model the sequence of evolutionary steps, including rock agglomeration, core formation, gas accretion, planetary migration, and so on.

This paper is organized as follows. We first consider the case of a two planet system in Section 2. We find the resulting tidal equilibrium state and show that it is a minimum. This treatment is then generalized in Section 3 to include systems with multiple planets. We show that adjacent pairs of planets can attain minimum energy states, and that the resulting system architecture of the equilibrium configuration provides a good description of observed planetary systems. Using a sample of observed multi-planet systems, Section 4 considers how well the energy scales found in observed pairs conform to the theoretical expectations of pair-wise equilibrium states. The paper concludes in Section 5 with a summary of results and a discussion of their implications. We also show that no global energy minimum exists if the total planetary mass can be distributed among three or more bodies (Appendix A). Finally, we present an order of magnitude estimate for the expected planetary masses based on the concept of pebble accretion (Appendix B).

2 Two Planet Systems

Here we consider a system containing two planets. The goal is to find the lowest energy state of the two planet system, subject to conservation of planetary mass and total system angular momentum. The mass fractions of the planets are allowed to vary: The idea is to find the state where the planets optimize energy during their formation process.

The planetary spins are assumed to have negligible angular momentum, so that the total angular momentum of the system is constant and is determined by the orbital angular momentum of the two planets. Without loss of generality, we can take the angular momentum vector of the inner planet to point in the ${\hat{z}}$ direction, so that the total angular momentum can be written in the form

[TABLE]

where $h_{j}$ is the orbital angular momentum of each planet. The inclination angle $i$ is measured from the pole. Here we take the magnitude of the angular momenta to be given by the expression in the low mass limit, i.e.,

[TABLE]

where $M_{\ast}$ is the stellar mass, $m_{j}$ is the planetary mass, and $e_{j}$ is the orbital eccentricity (i.e., we assume that $m_{j}\ll{M_{\ast}}$). The energy of the system is given by

[TABLE]

Note that, in general, the orbits for each of the two planets are described by six orbital elements. This treatment considers only the energy and angular momentum of the orbits. The interactions between the planets are assumed to determine the orbital spacing (through the parameter $\Lambda$), but are not included in the energy budget for the optimization procedure. This correction to the energy is of order $m_{j}/[M_{\ast}(\Lambda-1)]\sim 10^{-4}$.

Here we invoke a constraint on the total mass of the planetary pair, such that

[TABLE]

but allow the individual masses to vary. Motivated by the tendency for planetary migration to push planetary orbits together, we enforce a second constraint on the orbital spacing. We thus consider systems with a fixed separation, corresponding to a given number $K$ of Hill radii $R_{H}$, so that

[TABLE]

where

[TABLE]

The minimum value of $K$ required for stability in two-planet systems is generally taken to be $K=2\sqrt{3}$ (Gladman, 1993), whereas larger values are required for systems with more planets. Numerical simulations suggest that $K\sim 10$ or larger is needed for long-term stability (Chambers et al., 1996), where the value depends on the number and masses of the planets (Obertas et al., 2017; Wu et al., 2019). For comparison, the observed separations of planetary orbits in the Kepler sample display a wide distribution with a broad peak in the range $K=10-20$ (e.g., compare Weiss et al. 2018a and Pu & Wu 2015). Accounting for variations in the masses of the host stars, this range in $K$ corresponds to spacing parameters in the approximate range $\Lambda\approx 1.2-1.8$ (see also Section 4.3 and Figure 4 below).

Next we define dimensionless quantities according to

[TABLE]

The energy can be written in terms of the energy scale $GM_{\ast}m_{\scriptstyle T}/2R_{\ast}$ and the angular momentum can be written in terms of the scale $m_{\scriptstyle T}(GM_{\ast}R_{\ast})^{1/2}$. After dividing out these scales, the dimensionless expressions for the energy and angular momentum become

[TABLE]

and

[TABLE]

2.1 Extremum of the Energy

For fixed orbital spacing (given by $\Lambda$), the system defined above is described by five variables: the semimajor axis $a$ of the inner planet, the mass fraction $f$ of the inner planet, the orbital eccentricities $e_{j}$, and the inclination angle $i$. The system is also subject to conservation of angular momentum. The problem is thus to find the values of the variables $(f,a,e_{1},e_{2},i)$ that minimize the energy $E$ subject to the constraint ${\bf L}$ = constant.

The first step is to to extremize the energy $E$ subject to the constraint that the angular momentum is constant using Lagrange multipliers $\lambda_{z}$ and $\lambda_{x}$ (one for each component of the angular momentum). For each variable $\xi_{k}$, we must enforce the optimization condition

[TABLE]

where the $\xi_{k}$ represent the variables $(f,a,e_{1},e_{2},i)$. We thus obtain five equations that specify the equilibrium state:

[1] For the semimajor axis $a$ we obtain

[TABLE]

[2] For the mass fraction $f$ we obtain

[TABLE]

[3] For the orbital eccentricity $e_{1}$ of the inner planet, the optimization condition takes the form

[TABLE]

[4] Similarly, for the orbital eccentricity $e_{2}$ of the outer planet, we obtain

[TABLE]

[5] Finally, for the inclination angle $i$, we find

[TABLE]

The optimization conditions can be solved to find $e_{1}=0=e_{2}$, $i=0$, and $\lambda_{x}=0$. The optimal solution thus corresponds to circular orbits in the same orbital plane. The remaining equations (11) and (12) can both be solved for $\lambda_{z}$:

[TABLE]

and

[TABLE]

Combining these two expressions, we can solve for the mass fraction $f$ to obtain

[TABLE]

In the limit $\Lambda\to 1$ ($\Delta\to 0$), we obtain the expansion

[TABLE]

which corresponds to the mass being (nearly) equally distributed between the two planets. Moreover, even for relatively large values of the separation parameter $\Delta$, the fraction $f$ remains close to 1/2. In the opposite limit $\Lambda\to\infty$, the mass fraction $f\to 1/3$. For completeness, we can also find the ratio of planetary masses:

[TABLE]

Although the main result of this analysis is that the planetary masses must be nearly equal when the system resides at its critical point, to higher order we find that the inner planet is somewhat smaller than the outer planet. This state of affairs is seen in the sample of extrasolar planetary systems: The radii of observed planets within a given system tend to increase (slowly) with their orbital separation (Kipping, 2018), e.g., the outer planet is larger for $\sim 65\%$ of the observed cases (Weiss et al., 2018a).

Given the solution (18) for the mass fraction, and the optimization conditions ($e_{1}=e_{2}=i=0$), the magnitude $L$ of the angular momentum is related to the semimajor axis according to

[TABLE]

The lowest energy state can then be written in the form

[TABLE]

Note that in the limit of large $\Lambda\gg 1$, this expression for the energy has the asymptotic form $E\propto-\Lambda$. In other words, the system can always evolve to a lower energy state by spreading out the planets, i.e., moving one planet inward and the other outward in such a way as to conserve angular momentum.

In the limits where either $f\to 0$ or $f\to 1$, the energy of the system becomes $E=-1/L^{2}$. The depth of the energy minimum for the equilibrium state (where $f$ is given by equation [18]), can thus be measured by the quantity

[TABLE]

The final expression shows the result in the limit where the spacing factor $\Lambda$ is close to unity so that $\Delta\ll 1$.

2.2 Second Variation

The above analysis defines the critical point of the system. To show that the critical point corresponds to a minimum of the energy, we must consider the second variation. As is well known, the critical point will be a minimum if and only if the eigenvalues of the Hessian matrix $\mathbb{H}$ are all positive (Hesse, 1872), where the matrix elements are given by

[TABLE]

This section shows that the eigenvalues are indeed positive, so that the critical point represents the energy minimum.

To evaluate the Hessian matrix, we explicitly incorporate conservation of angular momentum into the expression for the energy. For convenience, we define reduced dimensionless angular momenta for the orbits:

[TABLE]

The magnitude of the total angular momentum is then given by

[TABLE]

so that the energy takes the form

[TABLE]

This formulation of the problem explicitly conserves angular momentum and eliminates the semimajor axis $a$ (of the inner planet). The energy is thus a function of four variables $E=E(f,e_{1},e_{2},i)$.

We now re-calculate the derivatives and the second derivatives. In this case, all of the mixed derivatives are zero when evaluated at the critical point, so that the eigenvalues are given by the diagonal terms of the Hessian matrix. Notice also that the factor of $L^{2}$ in equation (27) represents a multiplicative constant, and can be ignored for purposes of testing for stability (we thus set $L=1$ for the following analysis).

The derivative of the energy with respect to the mass fraction $f$ takes the form

[TABLE]

Solving this equation at the critical point where $\partial{E}/\partial{f}$ = 0 (and where $e_{j}=0=i$ from before), we obtain the mass fraction of equation (18), as expected. The second derivative, evaluated at the critical point, can be written

[TABLE]

Using the solution for the mass fraction $f$ at the critical point, this expression simplifies to the form

[TABLE]

which is manifestly positive.

Similarly, for the inclination angle the first derivative has the form

[TABLE]

We thus obtain $i=0$ as before. The second derivative becomes

[TABLE]

This quantity if positive for all values of $f$, including the value at the critical point given by equation (18).

The derivatives for the eccentricities take the form

[TABLE]

and

[TABLE]

The derivatives vanish for $e_{1}=0$ and $e_{2}=0$. The corresponding second derivatives become

[TABLE]

and

[TABLE]

both of which are positive.

In summary, we find that the Hessian matrix is diagonal, with all four eigenvalues manifestly positive. As a result, the critical point is a minimum of the energy. In this minimum energy configuration, the planets have circular orbits ($e_{j}=0$) in the same plane ($i=0$), with a mass ratio close to unity. More specifically, the mass fraction $f$ of the inner planet varies slowly with the orbital separation and is given by equation (18).

The existence and depth of the energy minimum is illustrated in Figure 1, which shows the energy for a two planet system as a function of the mass fraction $f$ for the inner planet. In this case, the eccentricities and inclination angle are set to zero, the values appropriate for the equilibrium state. With the semimajor axis determined by the requirement of conservation of angular momentum, the energy is a function of the mass fraction $f$ only for fixed orbital spacing $\Lambda$. The collection of curves in the figure shows the dimensionless energy versus $f$ for orbital spacings in the range $\Lambda$ = 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, and 3.0 (from top to bottom). Each curve has a minimum corresponding to the lowest energy state. The locus of minima is marked by the red dotted curve, where the value of mass fraction $f$ for the minimum is given by equation (18).

3 Systems with Three or More Planets

The next step is to consider systems with more than two planets. However, the optimization procedure used for the two planet model cannot be immediately generalized for three planets. If we allow for the total mass to be distributed among three bodies, then no solution for the energy minimum exists. This result is shown in Appendix A. The equilibrium state requires two of the planets to occupy the same orbit, thereby reducing the system to having only two planets. This argument can be generalized to systems with more than three planets. As a result, multiple systems with $N\geq 3$ planets cannot reside in a configuration that represents a global minimum of energy (subject to constant angular momentum). Some additional constraint must be placed on the system. This section introduces the concept of pairwise equilibrium, where each pair of planets on adjacent orbits resides in the tidal equilibrium state for the two planet problem outlined in the the previous section.

We start with three planet systems, where the three bodies are considered to be the sum of two sets of two planets, where the middle planet belongs to both the inner pair and the outer pair. If the energy is minimized subject to conservation of angular momentum of both the inner pair and (separately) the outer pair, then a solution can be found. We denote this state of the system as pairwise equilibrium. For both planetary pairs, the orbital spacing is assumed to be fixed, i.e., determined by physics independent of the optimization procedure.

Note that each of the two pairs of planets has a fixed orbital spacing for purposes of determining the (individual) equilibrium states, but two possibilities exist: In the first case, the ratio $\Lambda$ of the semimajor axes is the same for both pairs. In the second case, the number $K$ of mutual Hill radii is the same. In this second case, the ratio $\Lambda$ could be different for the two planetary pairs because the masses are not necessarily the same ($m_{1}+m_{2}\neq m_{2}+m_{3}$). However, the difference between the two scenarios leads to corrections that are higher order in the parameter $\Delta$, so we consider the case of constant $\Lambda$ for simplicity.

In this case, for fixed orbital spacing $\Lambda$, the solution for the mass fraction $f$ will be the same for both pairs, and can be written in the form given by equation (18). For the stationary state, each planet has an orbital angular momentum $L_{k}$. The ratio of the angular momenta of the two inner planets is given by

[TABLE]

Similarly, the second pair of planets must obey the relation

[TABLE]

The total angular momentum is thus given by

[TABLE]

Given the form of the angular momentum expression, consistency also requires that

[TABLE]

This expression keeps track of the fact that the inner orbit for the outer problem must be the outer orbit of the inner problem (hence the factor of $\sqrt{\Lambda}$), and allows for the masses of the inner planet pair and the outer planet pair to be different. Using this second expression, we can also write the total angular momentum in the form

[TABLE]

Equating the two expressions (39) and (41) for the total angular momentum, we can derive a relationship between the functions $F$ and $G$, i.e.,

[TABLE]

For completeness, note that the above equation is quadratic in $F$ and has two roots. The second solution $F=-1$ allows for equality for any value of $G$, but is unphysical.

For a given spacing of orbits $\Lambda$, the mass ratio of the inner and outer pairs has the form

[TABLE]

where we have defined the parameter $\eta$ in the final equality. After some algebra, the mass fractions for the three planets can be written in the form

[TABLE]

We can gain some insight by writing these expressions to leading order in $\Delta$, where the mass fractions become

[TABLE]

Since $\Delta$ is small in practice (typically $\Delta\sim 1/2$, as discussed below), the mass ratios are close to unity.

This argument can be readily generalized for more than three planets. For constant orbital spacing $\Lambda$, each successive angular momentum is larger than the previous one by a factor of $F(\Lambda)$ from equation (37), so that $L_{n+1}=F(\Lambda)L_{n}$. Each successive pair of angular momenta will be larger than the previous pair by a factor of $G$ from equation (40). Provided that the system with $N$ planets has $G=F$ (as shown above for $N=3$), then it can be shown through a straightforward induction argument that the final pair in a system of $N+1$ planets will also have $G=F$. As a result, in a state of pairwise equilibrium, each planet is larger than its inner neighbor by a factor of $\eta$ given by equation (43). The mass fractions for the planets in an $n$ planet system can thus be written in the form

[TABLE]

4 Comparison to Observed Planetary Systems

This section uses the properties of the observed multiple planet systems, primarily discovered via the Kepler mission, to compare with the theoretical expectations for systems in pairwise equilibrium. After identifying the observed planetary sample of interest, we construct the equivalent surface density distribution for multiple-planet systems (Section 4.1). The resulting mass profile differs from the minimum mass solar nebular (MMSN), but is consistent with that predicted from pairwise equilibrium. We then assess the energy budget for planetary systems (Section 4.2) and show that the energy minimized through pairwise equilibrium states is generally larger than the self-gravity of the planets and the interaction energy of the planet-planet potential. Finally, we assess the distribution of orbital spacings found in observed planetary systems and derive an order-of-magnitude estimate for the expected range (Section 4.3).

A number of previous studies have made detailed statistical analyses of the Kepler multi-planet systems of interest (e.g., see Millholland et al. 2017; Weiss et al. 2018a, b, and references therein). For comparison with the theoretical expectations derived in the previous sections, we use a simplified version of this observational sample, which is constructed as follows: From the publicly available database111https://exoplanetarchive.ipac.caltech.edu (where the data are taken from a variety of sources, e.g., Batalha et al. 2011), we consider the 219 planetary systems with 3 or more planets. This data set thus includes a total of 777 planets. Most of these planets are discovered through transit observations, some are detected through radial velocity measurements, and some have both types of data. For this simplified treatment, the planets with measured radii but without reported masses are assigned masses through the empirical relation $m_{\rm p}\sim R_{\rm p}^{2}$. Similarly, for the planets with measured masses but no reported radii, we use the inverse of this relation to specify $R_{\rm p}$. Although this approximate relation is sufficient for the rough comparison used here, one should keep in mind that planets with a given radius can have a range of masses (e.g., Wolfgang et al. 2016). For the stars that do not have masses listed, we use the main-sequence relationship between stellar mass and effective temperature to assign a mass (which is required to estimate the semimajor axis from the measured period). We thus obtain the set of variables $(m_{\rm p},R_{\rm p},a,P_{\rm orb},M_{\ast})$ for each planet.

Although the observational sample includes all of the multi-planet systems that are currently detected (with $N_{\rm p}\geq 3$), the data are both heterogeneous and subject to observational bias. The sample includes planetary systems orbiting stars with a range of masses. However, it is likely that planet formation proceeds differently around solar-type stars and low-mass stars (e.g., Laughlin et al. 2004). Unfortunately, with only $\sim 200$ systems, not enough data exist to construct separate samples for different stellar masses. The current collection includes planets detected by both transit observations and radial velocity measurements, and the estimates for planetary masses are subject to well-known uncertainties. In addition, planets are more readily detected in close orbits, so that the sample is likely to be missing planets with longer periods. In the present application, however, we are primarily interested in the properties of planetary pairs within the same system. These pairs have the same stellar host and are almost always detected by the same method, so that the aforementioned complications are partially mitigated.

4.1 Surface Density for Observed Multiple Planet Systems

If planetary systems achieve or approach the minimum energy states derived in the previous sections, then their corresponding surface density distributions must have a particular form. Here the surface density is that produced by smearing out the planetary masses into a smooth profile (e.g., see Chiang & Laughlin 2013). This subsection constructs the mass profile and surface density distribution predicted by pairwise equilibrium states, which have nearly equal mass planets and regularly spaced orbits. Such a configuration has a surface density profile $\sigma\sim a^{-2}$, as shown below. We also include a comparison to the surface density profile constructed from observed multi-planet systems (see Figure 2).

For a given orbital spacing $\Lambda$, the planetary masses for a pairwise equilibrium state are nearly equal. As a result, at the location in the disk corresponding to the $N$th planet, the semimajor axis and enclosed mass are given by

[TABLE]

where $a_{1}$ and $m_{1}$ are the semimajor axis and mass of the first (innermost) planet. For purposes of finding a benchmark mass distribution for comparison, we consider $N$ to be a continuous variable, so that

[TABLE]

The surface density corresponding to this mass distribution can be found from the defining relation

[TABLE]

This surface density distribution is thus significantly steeper than the usual MMSN where $\sigma\sim r^{-3/2}$ (compare with Weidenschilling 1977; Hayashi 1981; Desch 2007) and steeper than the result $\sigma\sim r^{-1}$ often inferred for circumstellar disks associated with newly formed stars (e.g., Andrews et al. 2009; Pérez et al. 2012). This profile is also steeper than that of the minimum mass extrasolar nebula (Chiang & Laughlin, 2013), which was constructed from the entire sample of exoplanets, not just those found in multiple planet systems.

It is important to keep in mind that equation (49) corresponds to the effective surface density distribution for the rocky component of planets at their current locations. If planetary systems evolve toward lower energy states during the process of planet formation by distributing the available mass in an optimal manner (see Sections 2 and 3), then they will attain nearly equal masses and the final, cumulative mass distribution of equation (48).

Figure 2 compares the mass profile expected from a pairwise equilibrium state with that constructed from the observed bodies in multiple planet systems. The blue dashed curve shows the mass distribution $M(a)\sim\ln(a)$ from equation (49), which results from equal mass planets with uniform spacing (where we specify the inner radius $a_{1}$ = 0.04 AU). The red curves shows the cumulative mass distribution for the observed sample. This observational profile was constructed as follows: The planets in the sample were first ordered according to their semimajor axes. Since we are interested in the mass profile of the rocky component, planets with masses larger than 15 $M_{\oplus}$ were assumed to have only 15 $M_{\oplus}$ of rocky material. The mass profile is then given by the straightforward expression

[TABLE]

where $(a_{k},m_{k})$ are the masses and semimajor axes of the observed planets (subject to $m_{k}\leq 15M_{\oplus}$) and ${\cal H}(x)$ is the Heaviside step function (Abramowitz & Stegun, 1972). Only planets with periods less than 100 days are included. Both the theoretical and observed mass distributions are then normalized to unity at the outermost value of $a$. Notice that the observed distribution is roughly consistent with that expected from pairwise equilibrium. One interesting point of departure occurs for small radii a\,\raise 1.29167pt\hbox{<}\kern-7.5pt\lower 3.01385pt\hbox{\sim}\,0.04 AU, where the observational profile shows a weak power-law tail; this portion of the distribution contains an excess of $\sim 7\%$ of the total mass. Both the observed mass distribution and that predicted by pairwise equilibrium are clearly distinct from the power-law form of the MMSN. This latter profile is depicted in Figure 2 by the dotted line, which has the form $M(a)\sim{a^{1/2}}$, corresponding to the surface density distribution $\sigma(a)\sim{a^{-3/2}}$.

The similarity between the observed mass profile and the theoretical construction of equation (49) is not unexpected: This agreement means that the observed sample of multi-planet systems displays nearly equal mass planets with regularly spaced orbits as claimed previously (e.g., Rowe et al. 2014; Millholland et al. 2017; Weiss et al. 2018a, b). If we enforce regular orbital spacing, then the condition of pairwise equilibrium requires nearly equal planetary masses, as observed.

Nonetheless, the observational profile of Figure 2 is subject to a number of uncertainties: Given the limited number of planets detected in multi-planet systems (777 in this sample), we have included planets orbiting host stars of all masses. The surface density profiles could display some dependence on stellar mass, and that trend is not captured by this compilation. This construction also implicitly assumes that the planetary systems are complete out to periods $P=100$ days. Additional planets, as yet undetected, could also be present in these systems, and would alter the inferred mass distribution. More specifically, biases in planet detection favor the discovery of short-period planets. The apparent absence of additional long-period planets could contribute to the inferred steepness of the surface density distribution found here.

4.2 Energy Scales

The underlying idea of pairwise equilibrium is that planetary systems can attain lower energy states by optimizing their mass ratios, where the critical point corresponds to nearly equal masses (specifically, the mass fraction given by equation [18]). The energy available through this optimization, presumably approached during the epoch of planet formation, is illustrated in Figure 1. However, the analysis thus far has not included the self-gravity of the planets themselves or the potential energy due to gravitational interactions between planets. This section assesses the relative sizes of these additional energy contributions and shows that they are sub-dominant for the observed multi-planet systems of interest.

The energy available due to re-arrangement of planetary mass is given by equation (23), where this energy is scaled relative to the orbital energy. The dimensionless self-gravity of the planets, again scaled relative to their orbital energy, is given by

[TABLE]

where $\alpha_{\rm g}$ is a dimensionless constant of order unity and depends on the interior structure of the planet. Finally, the (dimensionless) energy of interaction between planets in adjacent orbits can be written in the form

[TABLE]

where $\beta$ is another dimensionless factor of order unity (see Murray & Dermott 1999 for further detail).

The energy scales for adjacent planetary pairs residing in observed multi-planet systems are depicted in Figure 3. All of the energies are scaled to the orbital energy and are thus dimensionless. The upper collection of points shows the energy available from pairwise equilibrium (using equation [23]) as a function of the summed self-gravity of the two constituent planets (see equation [51]). The dark blue points show the results for pairs of rocky planets, whereas the cyan points correspond to planets with higher masses, namely $(m_{1}+m_{2})>30M_{\oplus}$. The solid black curve marks the line of equality. Notice that the majority of the blue points fall above the solid curve, so that the pairwise energy dominates over that of self-gravity. In contrast, more cyan points fall below the solid curve than above it. This result suggests that the self-gravity of the planets is relatively unimportant during the formation process for rocky planets (blue points), but self-gravity plays a more important role in the assembly of the largest planets (cyan points). Moreover, the neglect of self-gravity represents a valid approximation for rocky planets in multiple systems, whereas it should be included in systems containing larger Jovian planets. The analysis of this paper is thus limited to the former type of planetary system.

This treatment has ignored both the energy and angular momentum of the planetary spins. Although the rotation rates of exoplanets are generally not known, the rotational speeds are constrained to be less than break up velocity. As a result, the rotational kinetic energy must be smaller than that of self-gravity, which in turn is much less than the pairwise equilibrium energy considered in this paper, as shown in Figure 3. The neglect of spin in the energy budget is thus justified.

Figure 3 also shows the interaction energy between pairs of planets (from equation [52]) as a function of self-gravity. The red points show the results for rocky planets, whereas the yellow points correspond to the larger planets $(m_{1}+m_{2}>30M_{\oplus})$. The red points fall well below the line of equality. Moreover, the red points (self-gravity) fall even farther below the blue points (pairwise energy). This latter result can be interpreted as follows: The leading order expression for the pairwise energy has the form $E_{\rm pair}\approx 3\Delta^{2}/16$, and that for the interaction energy becomes $E_{\rm int}\approx{m_{\rm p}}/(M_{\ast}\Delta)$. These two (dimensionless) energies are equal when the planetary separation decreases to the value

[TABLE]

In other words, the spacing of adjacent orbits must be as small as $K=2^{4/3}$ for the interaction energy to dominate. Of course, given that these systems experience a large number of dynamical time scales, planetary interactions can affect system stability for wider spacing. As discussed previously, observed systems have $K\approx 10-20$, and K\,\raise 1.29167pt\hbox{>}\kern-7.5pt\lower 3.01385pt\hbox{\sim}\,10 is required for long-term stability. The finding that the pairwise energy is much greater than the interaction energy is thus expected from considerations of orbital stability and from the observed spacing of planetary orbits.

4.3 Orbital Spacing

The discussion thus far has considered the orbital spacing $\Lambda$ to be fixed. In this section, we consider the observed distribution of orbital spacings and place approximate constraints on the expected range of values. One criterion is that the planets must not be too close together in order to remain stable. In addition, if the planets are spaced as tightly as possible subject to stability constraints — within a few Hill radii of each other — the disk will not have enough mass to form them at their final locations. As a result, we can obtain a rough estimate for orbital spacing based on considerations of the total mass budget in the disk. Note that this order of magnitude estimate does not include possibile redistribution of mass within the disk. Such rearrangement could take place by moving raw material through the disk before the epoch of planet formation and/or by changing the orbital elements of assembled planets through migration.

The mass scale for rocky planets is observed to be of order $m_{\rm p}\approx 10M_{\oplus}$ (Zhu et al., 2018) and is remarkably consistent within planetary systems. This mass scale can be derived from considerations of pebble accretion taking place within the circumstellar disks that gives birth to planets, as outlined in Appendix B (see also Lambrechts et al. 2014; Bitsch et al. 2015 and references therein). Since the total mass of the disk is finite, the total number of rocky planets that can be produced is limited. The rocky mass within a disk can be written in the form

[TABLE]

where $Z\sim 0.01$ is dust to gas mass ratio and we have assumed that the maximum disk mass is about 10% of the stellar mass. This latter fraction can be understood in terms of disk stability (Shu et al., 1990) and is consistent with the upper limits on disk masses in observed systems (Hartmann, 2007). Under the assumption that all of the heavy elements in the disk can be utilized for planet formation, the maximum number ${\cal N}$ of rocky planets is given by

[TABLE]

Combining this result with the mass profile from equation (49), we can solve for the spacing parameter to find

[TABLE]

where $a_{\rm out}$ and $a_{\rm in}$ represent the outer and inner edges of the disk. For typical disk parameters $a_{\rm out}$ = 100 AU and $a_{\rm in}$ = 0.05 AU, we find $\Lambda\approx 1.3-1.6$.

Figure 4 shows the observed distribution of orbital spacings for the planetary pairs in our comparison sample. For the sake of definiteness, this compilation includes only pairs with total mass $m_{1}+m_{2}\leq 30M_{\oplus}$, roughly corresponding to superearth planets and thereby excluding Jovian planets. The observed spacing distribution (blue histogram) shows a peak near $\Lambda\sim 1.4$. The vertical green lines in the figure show the benchmark values $\Lambda$ = 1.3 and 1.6 derived above on the basis of the available mass in rocky material. About 54% of the pairs in the observational sample fall in this interval ($1.3\leq\Lambda\leq 1.6$), whereas about 77% of the pairs fall within the wider range $1.2\leq\Lambda\leq 1.8$. The observed distribution of orbital spacings in these systems is roughly consistent with that expected from considerations of the mass supply. Nonetheless, the distribution displays a tail corresponding to wider spacings, where $\sim$13% of the pairs have $\Lambda>2$. Such large values of $\Lambda$ could be explained by missing planets, planetary migration, and/or the rearrangement of disk material before planets are formed (see, e.g., Hansen & Murray 2013). In any case, the full spacing distribution requires additional explanation and should be considered in future work.

5 Conclusion

This paper formulates and solves a new type of Darwin problem that determines the lowest energy (tidal equilibrium) states for multi-planet systems subject to a collection of constraints. The optimized system architectures have properties that are roughly consistent with those of observed planetary systems. This section presents a summary of our specific results (Section 5.1) and a discussion of several unresolved issues (Section 5.2).

5.1 Summary of Results

We first consider two planet systems subject to constraints of constant angular momentum and fixed orbital spacing (Section 2). The minimum energy state of the system corresponds to planets with nearly equal masses on circular orbits ($e_{j}=0$) in the same plane ($i=0$). The mass fraction $f\approx 1/2$, so that $m_{1}\approx{m_{2}}$, and varies slowly with orbital spacing as given by equation (18). These properties are realized in the observed sample of extrasolar planetary systems found by the Kepler mission.

For systems with more than two planets, where the total mass can be freely distributed among the planets, no (global) minimum energy state exists (Appendix A). The system can evolve to a lower energy state by consolidating the planets, thereby reducing their number. Observed multiple planet systems do not reside in states corresponding to global energy minima, so that some additional consideration must act to define their properties.

For multiple planet systems, we introduce the concept of pairwise equilibrium, where each pair of planets in adjacent orbits attains the minimum energy state appropriate for the pair (Section 3), subject to conservation of angular momentum and fixed orbital spacing. With each pair in this type of equilibrium, the resulting planetary system has properties much like those observed in the Kepler sample: circular orbits confined to a single plane, with regular orbital spacing and nearly uniform planetary masses. To higher order, the planetary mass is a slowly increasing function of semimajor axis (see equations [43] and [44]).

The optimal configurations found in this paper can be compared to observed planetary systems (Section 4). The pairwise equilibrium states have nearly equal mass planets on regularly spaced orbits. This type of system has an effective surface density of the form $\sigma\propto$ $r^{-2}$ and mass profile $M(r)\propto\ln(r)$, which differs substantially from the usual form $\sigma\propto{r^{-3/2}}$ of the MMSN or the form $\sigma\propto{r^{-1}}$ inferred for circumstellar disks. The rocky component of observed multi-planet systems have similar surface density and mass profiles for orbital periods $P<100$ days (see Figure 2). These profiles arise from nearly equal mass planets on regularly spaced orbits, a configuration that has been suggested in previous studies (Rowe et al., 2014; Pu & Wu, 2015; Millholland et al., 2017; Wang, 2017; Weiss et al., 2018a, b).

In the multiple-planet systems of interest, the typical mass $m_{\rm p}\sim 10M_{\oplus}$ can be understood as the isolation mass produced through pebble accretion (Appendix B). Given this planetary mass scale, the energy difference that can be realized by systems approaching a tidal equilibrium state is larger than both the self-gravity of the individual bodies and the interaction energy between planets (see Figure 3). As a result, during the process of planet formation, considerations of pairwise equilibrium are energetically important. For Jovian planets, however, the self-gravity plays a more dominant role. This apparent dichotomy is consistent with the previous finding that intra-system uniformity of planetary masses breaks down in the presence of giant planets (Wang, 2017).

5.2 Discussion

This paper has shown that the minimum energy configurations available to multi-planet systems require nearly equal planetary masses, along with circular orbits in the same plane. Although this finding is consistent with the observed properties of planetary systems, it provides only a partial explanation, and a number of issues remain unresolved:

The optimization procedure of this work shows that the critical state for two planet systems has the mass fraction $f\approx 1/2$ given by equation (18). The determination of planetary masses takes place (by definition) during the epoch of planet formation. However, the existence and properties of the minimum energy state do not explain how the process of planet formation realizes this goal. On a related note, the assumption of constant mass for planetary pairs is an idealization. Moreover, in planetary systems containing giant planets, the self-gravity of the planets dominates the pairwise effects considered in this paper, and observations show that such systems no longer display the same regular features (e.g., Wang 2017). As a result, some process — yet to be determined — drives some forming planetary systems to produce giant planets while other systems produce multiple superearths.

Another unresolved issue concerns the orbital spacing of the planets. This work shows that planetary systems prefer nearly equal masses for fixed orbital spacing, but does not provide an explanation for the size of the observed spacing. Previous work shows that orbits cannot be too close without rendering the systems dynamically unstable (e.g., Gladman 1993; Pu & Wu 2015), thereby providing a lower limit on orbital spacing. Considerations of statistical mechanics can produce consistent spacing distributions (Tremaine, 2015), but require particular choices for the relevant parameters (see also Mogavero 2017; Pakter & Levin 2018). Using limits on the available disk mass and radius, this paper constructs estimates for the orbital spacing parameters that are roughly consistent with those observed (Section 4.3). The spacing could be tighter than these estimates, however, if some migration process moves planets together after they are formed, or transports rocky material inward before it becomes incorporated into planets. This work also indicates that equal mass planets on more widely separated orbits are energetically preferred. Nonetheless, we would like to know the details concerning how the process of planet formation produces the observed distribution of orbital spacing (see Figure 4).

Planets with equal masses and regular spacing, as found in both observed systems and the minimum energy configurations of this paper, correspond to an effective surface density profile $\sigma(r)\sim$ $r^{-2}$, which is steeper than that of both observed circumstellar disks and the MMSN. As a result, some process must redistribute the rocky mass within these systems, either as raw material before (perhaps during) planet formation or later as planetary bodies with their final masses. The timing and mechanism(s) of this redistribution process also remains unresolved. In addition, the observed mass distribution shows a power-law tail inside the nominal inner radius of the disk (see Figure 2). The manner in which this region is populated by planets poses another interesting problem.

The model of pairwise equilibrium developed here applies to compact planetary systems, such as those observed by Kepler. However, we currently have limited sensitivity to additional planets on wider orbits. If these planetary systems can reach the energy-optimized states considered here, then the masses of additional exterior planets would slowly increase with semimajor axis (for fixed orbital spacing; see Section 3; see also Jiang et al. 2015). Reaching such equilibrium states requires the planet formation process to communicate across the entire planetary system. Future developments in planet formation theory and extended observations of planetary systems are necessary to determine if this is the case. On a related note, the observed systems could have inner planets that are closely spaced, with a wide gap separating a collection of additional planets. In this context, the two (separated) sets of planets could reach separate states of pairwise equilibrium. If the outer planets are sufficiently massive, however, then their self-gravity would dominate over the pairwise energy, and the systems are less likely to be well-ordered.

Finally, we note that the Galilean satellites also display nearly uniform masses, regular orbital spacing, low eccentricities, and coplanar orbits aligned with the equatorial plane of Jupiter. This remarkable regularity suggests that these moons represent a miniature solar system, whose formation requires an explanation (e.g., Lunine & Stevenson 1982). The pairwise equilibrium solutions of this paper lead to similar configurations and could thus operate on smaller mass and size scales.

The open issues outlined above, along with the diversity of the planetary systems that are ultimately produced, highlight the complicated nature of the planet formation process. In the face of such complexity, one way forward is to identify some type of organizing principle that allows us to understand and/or constrain the possible outcomes. The pairwise minimum energy configurations found in this paper provide one avenue toward this end.

Acknowledgments: We would like to thank Konstantin Batygin, Juliette Becker, Tony Bloch, Greg Laughlin, Erik Petigura, Darryl Seligman, and Chris Spalding for useful discussions. We also thank an anonymous referee for constructive comments. This work was supported by the University of Michigan.

Appendix A Absence of Equilibrium States for Three Planet Systems

This Appendix considers the possible equilibrium state for a three planet system in which the mass can be distributed among the three planets and no pairwise constraints are applied. In this case, no equilibrium state for the three planet system exists. Equivalently, the equilibrium state requires the middle planet to have the same orbit as the inner planet, thereby reducing the system to two bodies.

For three planets, the mass fractions can be written in the form

[TABLE]

where $m_{\scriptstyle T}$ is the total mass of the three planets. Here we specialize to the case of circular orbits, all confined to the same orbital plane. The semimajor axes of the orbits are then given by

[TABLE]

where the factors $\Lambda$ and $\Gamma$ are considered fixed. For this reduced problem, the dimensionless energy and angular momentum take the form

[TABLE]

and

[TABLE]

To minimize the energy subject to conservation of angular momentum, we must find the critical point of the function $F=E+\lambda L$, where $\lambda$ is the Lagrange multiplier. The derivatives for the mass fractions have the form

[TABLE]

and

[TABLE]

Solving both equations for the quantity $\lambda a^{3/2}$, we find

[TABLE]

The second equality determines the required relation between the separation factors $\Lambda$ and $\Gamma$, and can be simplified to the form

[TABLE]

This expression implies that $\Lambda=1$. In other words, any solution for the lowest energy state of the three planet system requires the middle planet to have the same orbit as the inner planet, thereby producing a two planet system. For three separate planets, which requires $\Lambda>1$, no tidal equilibrium state exists.

This trend is illustrated in Figure 5, which shows the iso-energy contours for a three planet system as a function of the mass fractions $f$ and $g$ of the inner two planets (the mass fraction of the third planet is determined by conservation of mass to be $[1-f-g]$). The mass fraction $g$ of the middle planet approaches the $f$-axis for the minimum energy state, i.e., the mass of the middle planet vanishes. In this $g\to 0$ limit, the mass fraction $f$ of the inner planet approaches the value appropriate for the two planet system (from Section 2 and equation [18]). Note that the mass fraction of the outer planet is given by $1-f-g$, so that the physically relevant part of the diagram is confined to the region defined by $f+g\leq 1$.

Appendix B Characteristic Mass Scale for Rocky Planet Formation

This Appendix derives a characteristic mass scale for planet formation based on the idea that planets accumulate mass through pebble accretion (e.g., Lambrechts et al. 2014; Bitsch et al. 2015) which eventually shuts down when the pressure gradient produced by the planet becomes sufficiently large. As shown below, this benchmark mass scale falls in the range $m_{\rm p}=4-8M_{\oplus}$, consistent with planetary masses typically found in observed multi-planet systems.

The concept behind the isolation mass for pebble accretion is that the mass of the growing planet (or planetary core) becomes large enough to affect the pressure of the surrounding gas, and thereby provides a pressure gradient that inhibits further growth. In order of magnitude, application of the Bernoulli equation indicates that the planet will influence the flow of gas parcels by increasing their speed, and this increase will become significant when

[TABLE]

where $m_{\rm p}$ is the mass of the growing planet (core), $b$ is the impact parameter of the fluid parcel trajectory, and $c_{\rm s}$ is the sound speed in the gas.

The planet receives pebbles from within its Hill sphere, so we take the impact parameter to be comparable to the Hill radius. This choice allows us to write the stopping condition in the form

[TABLE]

where $\Omega_{\rm k}$ is the Keplerian rotation rate of the disk and where $\alpha_{\rm p}$ is a dimensionless factor of order unity. Solving for the planetary mass yields

[TABLE]

Next we want to write the isolation mass in terms of the scale height $H$ of the disk. The standard definition has the form

[TABLE]

so that

[TABLE]

Using the typical value $H/r\sim 0.05$, we can write this expression in the form

[TABLE]

For the observational sample of planets, the typical stellar mass $M_{\ast}\approx 0.5-1M_{\odot}$, so that the benchmark planetary mass is about $m_{\rm p}\approx 4-8M_{\oplus}$. This value is roughly consistent with the mass scale of planets observed in multiple systems (Zhu et al., 2018).