THE STATISTICS OF MULTI-PLANET SYSTEMS

The distribution of inclinations in multi-planet systems provides fundamental insights into planet formation. The small inclinations of the planets in the solar system—the largest is 7°, for Mercury—strongly suggest that they formed from a disk. However, we should not be surprised if other multi-planet systems have larger inclinations, for several reasons: (1) the rms inclinations in the asteroid and Kuiper Belts are substantially larger, 12° and 16° respectively; (2) in most astrophysical disks, the rms eccentricity and inclination are correlated, and the eccentricities of extrasolar planets are much larger than those of solar-system planets (mean of 0.25 for exoplanets with periods greater than 10 days, compared to 0.06); (3) a number of dynamical mechanisms can excite inclinations, including Kozai–Lidov oscillations, planet–planet scattering, and resonance sweeping; (4) measurements of the Rossiter–McLaughlin effect in transiting systems (e.g., Winn et al. 2010) show a broad distribution of obliquities (angle between the spin axis of the host star and the orbit axis of the planet) and some processes that excite obliquities do so by exciting inclinations; (5) most extrasolar planetary systems have quite different configurations from the solar system (e.g., giant planets 100 × closer to the star), so they may form by quite different mechanisms; (6) there are still serious theoretical obstacles to the formation of planets from a circumstellar disk, and several authors have suggested that some or all planets may be formed by other mechanisms, more similar to star formation (e.g., Black 1997; Papaloizou & Terquem 2001; Ribas & Miralda-Escudé 2007; Abt 2010).

There is only fragmentary evidence that most extrasolar planetary systems have small relative inclinations.

• 1.  
  The mutual inclination of planets B and C in the system surrounding the pulsar B1257+12 is less than ∼13° (Konacki & Wolszczan 2003); however, this result is only marginally relevant to planetary systems around main-sequence stars since pulsar planets must have had a very different history.
• 2.  
  Using radial-velocity (RV) and astrometric data, Bean & Seifahrt (2009) estimate that the mutual inclination between the planets GJ 876 b and c is $5\mbox{$.\!\!^\circ $}0{+ 3\mbox{$.\!\!^\circ $}9\atop -2\mbox{$.\!\!^\circ $}3}$. Using RVs and dynamical modeling of the planet–planet interactions Correia et al. (2010) conclude that the mutual inclination is ≲ 2°, while Baluev (2011) finds that the same quantity is between 5° and 15°. The large scatter among these results means that they should be used with caution.
• 3.  
  McArthur et al. (2010) find from astrometric and RV measurements that the mutual inclination of υ And c and d is 30° ± 1°, much larger than in GJ 876 but still small enough to suggest formation from a disk.
• 4.  
  Dynamical fits to the transit timing of two planets in the Kepler-9 system yield an upper limit to the mutual inclination of ∼10° (Holman et al. 2010). However, this system is likely to have smaller inclinations than a typical multi-planet system, since it was discovered in a transit survey and such surveys are far more likely to detect systems with small inclinations rather than large ones.
• 5.  
  Lissauer et al. (2011a) studied the six-planet system Kepler-11 and concluded that the absence of transit duration changes in Kepler-11e implies that its inclination relative to the mean orbital plane of the other planets is less than 2°; once again, this small inclination is unlikely to be typical because the probability that two or more planets will transit is a strong function of their mutual inclination. Lissauer et al. also concluded that the mean mutual inclination of the planets in Kepler-11 was small from a different argument: they used Monte Carlo simulations to assess the probability that a system of six planets with a given rms inclination would exhibit six transits to a randomly placed observer and found that this probability declined from its maximum (at zero rms inclination) by e^−1/2 when the mean inclination was 1°–2°. However, this argument is fallacious: the probability that a given star with six planets would show six transits is different from the probability that one star from the Kepler sample of over 100,000 stars would show six transits (this is sometimes called the prosecutor's fallacy: if the probability of a false match between DNA samples is 1 in 100,000, and crime-scene DNA is matched to someone in a database of 100,000 men, the probability that the suspect is innocent is of order unity, not 1 in 100,000 as a zealous prosecutor might argue).

As noted above, if one planet in a two-planet system transits its host star as viewed from Earth, the probability that the second planet will also transit is higher if the mutual inclination of the two planetary orbits is small (e.g., Ragozzine & Holman 2010). This argument suggests that the numbers of 1-planet, 2-planet,..., N-planet systems detected in a large transit survey contain information about both the multiplicity function—the fraction of host stars containing 0, 1, 2, ..., N planets—and the inclination distribution. The challenge is to disentangle the two distributions to distinguish thick systems with many planets from thin systems with few planets.

The first attempt to do this was made by Lissauer et al. (2011b), who modeled the number of multiple-planet systems detected in the first four months of data from the Kepler survey (Borucki et al. 2011)—115 with two transiting planets, 45 with three, 8 with four, and 1 each with five and six. Lissauer et al. used a variety of simple models for the distribution of the number of planets per system. They found that their models produced too few systems in which a single transiting planet was observed, but that the best-fit models typically had rms inclinations ≲ 10°.

The purpose of this paper is to develop a general formalism that relates the intrinsic properties of multi-planet systems to the properties of the multi-planet systems that are detected in transit or other surveys (Sections 2 and 3), and to apply this formalism to the Kepler planet survey (Section 4) and to RV surveys (Section 5). Previous analyses have used Monte Carlo simulations to explore these problems, but our calculations are mostly analytic or semi-analytic and do not employ Monte Carlo methods.

1.1. Preliminaries

First we introduce some notation. (1) The Kepler team uses the term planet "candidate" to denote a possible planet that has been discovered through transits but not yet been confirmed by RV measurements. Morton & Johnson (2011) estimate that 90%–95% of the Kepler planet candidates are real planets, so for the remainder of this paper we will simply assume that all the Kepler planet candidates are real and delete the word "candidate." (2) We must constantly distinguish between the number of planets in a system and the number of transiting planets in that system. We use the contraction "tranet" to denote "transiting planet." Thus one could have, for example, a two-tranet, three-planet system (Ragozzine & Holman 2010 call this a "double-transiting triple system"). (3) We distinguish two types of selection effects that limit a planet sample. Every survey has a set of detection thresholds, determined by the parameters of the survey that limit the properties of the planets that it can detect (maximum orbit period, minimum reflex RV, minimum transit depth, etc.). A survey selection effect is a limitation on the number of detectable planets due to the detection thresholds. A geometrical selection effect is a limitation arising from the orientation of the planetary system—in particular, the planet must cross in front of the stellar disk to be detectable in a transit survey.²

We assume that the stars in a survey may have 0, 1, ..., K planets and denote the number of stars in the survey with k planets by N_k. Thus ∑^K_{k = 0}N_k is the total number of stars in the survey. The vector $\mathbf {N}=(N_0,N_1,\ldots,N_K)$ is called the multiplicity function.

Because of survey and geometric selection effects, only a fraction of these planets will be detected in the survey. Let the survey selection matrix element S_km be the probability that a system containing m planets has k of them that pass the survey selection criteria. Similarly, let the geometric selection matrix G_jk be the probability that j of these k planets pass the geometric selection criteria. Then, the expected number of systems that the survey should detect with j tranets is

$$\begin{equation} \overline{n}_j=\sum _{k=j}^K G_{jk}\sum _{m=k}^KS_{km}N_m,\quad \mbox{or}\quad \overline{\mathbf {n}}=\mathbf {G}\cdot \mathbf {S}\cdot \mathbf {N}. \end{equation} \tag{ 1 }$$

We call $\overline{\mathbf {n}}$ the observable multiplicity function. Clearly,

$$\begin{equation} G_{mn}=S_{mn}=0\hbox{ for }m>n,\;\; G_{00}=S_{00}=1,\;\; G_{mn},S_{mn}\ge 0. \end{equation} \tag{ 2 }$$

Moreover, since the number of detectable planets in an n-planet system must be between 0 and n, we have

$$\begin{equation} \sum _{m=0}^nG_{mn}=\sum _{m=0}^nS_{mn}=1. \end{equation} \tag{ 3 }$$

Thus $\mathbf {G}$ and $\mathbf {S}$ are (K + 1) × (K + 1) upper-triangular stochastic matrices. For physical reasons $\mathbf {G}$ and $\mathbf {S}$ must commute (Equation (1) should not depend on whether we consider the survey selection effects or the geometric selection effects first). We have confirmed that this is indeed true for the selection matrices that we derive below.

1.2. Separability

Let $\mathbf {w}$ represent all of the orientation-independent properties of a planet that determine its detectability (planetary mass and radius; stellar mass, radius, distance, and luminosity; orbital period, etc.) and let $f(\mathbf {w}_1,\ldots,\mathbf {w}_n)$ represent the probability distribution of these parameters for an n-planet system. Thus $\int d\mathbf {w}_1\cdots d\mathbf {w}_n\,f(\mathbf {w}_1,\ldots,\mathbf {w}_n)=1$.

A natural approximation for describing multi-planet systems is that the n-planet distribution function is separable, that is,

$$\begin{equation} f(\mathbf {w}_1,\ldots,\mathbf {w}_n)=\prod _{m=1}^n f(\mathbf {w}_m), \qquad \int d\mathbf {w}\,f(\mathbf {w})=1. \end{equation} \tag{ 4 }$$

We describe the evidence on the validity of this approximation in Section 3.2.

Let $\Theta ^{A}(\mathbf {w})$ be the probability that a planet with properties $\mathbf {w}$ is detected in the survey labeled by A if its host star is on the target list for this survey and the orientation of the observer is correct (we assume that whether or not a planet can be detected is independent of the presence or absence of other planets in the same system, which is a reasonable first approximation). Thus, the function $\Theta ^{A}(\mathbf {w})$ describes the survey selection effects for A, but not the geometric selection effects. The probability that a planet is detected, ignoring geometric selection effects, is then

$$\begin{equation} W^{A}=\int f(\mathbf {w})\Theta ^{A}(\mathbf {w})\,d\mathbf {w}. \end{equation} \tag{ 5 }$$

If the survey target list contains N^A_m stars with m planets, then using the separability approximation (4) the expected number of systems in which k planets will be detected is

$$\begin{equation} \overline{n}^{A}_k=\sum _{m=k}^K S_{km}(W^{A})N^{A}_m, \quad 0\le k\le K; \end{equation} \tag{ 6 }$$

where the survey selection matrix $\mathbf {S}$ is a (K + 1) × (K + 1) matrix whose entries are given by the binomial distribution,

$$\begin{equation} S_{km}(W)\equiv \frac{m!}{k!(m-k)!}W^k(1-W)^{m-k}, \quad 0\le k\le m\le K, \end{equation} \tag{ 7 }$$

and zero otherwise. Note that $\mathbf {S}(1)$ is the unit matrix. A useful identity is

$$\begin{equation} \mathbf {S}(A)\cdot \mathbf {S}(B)=\mathbf {S}(AB), \end{equation} \tag{ 8 }$$

which in turn implies

$$\begin{equation} \mathbf {S}^{-1}(W)=\mathbf {S}(W^{-1}). \end{equation} \tag{ 9 }$$

Although the physical motivation (6) for the definition of $\mathbf {S}(W)$ requires 0 ⩽ W ⩽ 1, the matrix is well defined for all values of W.

Assuming separability, it is straightforward to show that the conditional probability distribution of the parameters $\mathbf {w}_m$, given that k planets are detected, is (cf. Equation (4))

$$\begin{equation} f(\mathbf {w}_1,\ldots,\mathbf {w}_k)=\prod _{m=1}^k f(\mathbf {w}_m). \end{equation} \tag{ 10 }$$

Thus, a separable distribution is still separable after survey selection effects are applied, so long as the selection effects depend only on the properties of an individual planet.

The factor W (Equation (5)) is difficult to determine reliably since (1) we do not have good models for the distribution $f(\mathbf {w})$ of the planetary parameters, (2) in most cases the survey selection effects $\Theta (\mathbf {w})$ are not known accurately, and (3) in many cases the target list from which a given sample of exoplanets was detected is not even known (the Kepler survey is an exception to the last two limitations). However, useful results can be obtained without an explicit evaluation of W. Suppose, for example, we have two surveys A and B that examine populations of target stars with similar characteristics, then the ratio of the number of m-planet systems in the target populations of the two surveys should be independent of m, so N^B_m = cN^A_m where c is a constant given by the ratio of the number of target stars in B and A. Equation (6) can then be written as

$$\begin{equation} \overline{\mathbf {n}}^{A}=\mathbf {S}(W^{A})\mathbf {N}^{A}, \quad \overline{\mathbf {n}}^{B}=c\,\mathbf {S}(W^{B})\mathbf {N}^{A}. \end{equation} \tag{ 11 }$$

Applying Equations (8) and (9), we have

$$\begin{equation} \overline{\mathbf {n}}^{B}=c\,\mathbf {S}(f^{BA})\overline{\mathbf {n}}^{A}, \end{equation} \tag{ 12 }$$

where f^BA ≡ W^B/W^A = 1/f^AB. Thus, the observable multiplicity function $\overline{\mathbf {n}}^{B}$ of survey B is directly related to that of survey A by a matrix that depends only on a single parameter f ^BA (the normalization constant c is known, since it is just the ratio of the number of target stars in the two surveys). Thus, the expected number of 1-, 2-, ..., planet systems in survey A can be determined from the number of 1-, 2-, ..., planet systems in B as functions of a single variable f ^BA. Alternatively, since f ^BA is not easy to determine from the data, we can choose $\overline{n}_1^{B}$ as the independent variable and plot $\overline{n}_2^{B},\overline{n}_3^{B},\ldots$ and f ^BA as functions of $\overline{n}_1^{B}$; such plots can be directly compared with the actual numbers n^B_k observed in survey A as a test of the separability approximation. In practice, we must use $\mathbf {n}^{A}$ (the actual multiplicity function of the systems detected in survey A) rather than the mean multiplicity function $\overline{\mathbf {n}}^{A}$ (the expected multiplicity function given survey A's properties) on the right side of Equation (12) but these should not be very different so long as n^A_k ≫ 1. Equation (11) is well defined whether f ^BA is larger or smaller than unity, but if f^BA > 1 the statistical errors will be amplified and it is likely that some of the predicted values of $\overline{n}_k^{B}$ will be negative, which is unphysical. Thus, if the separability approximation is valid, the observable multiplicity function of deep surveys can be used to predict the observable multiplicity function of shallow surveys (but not vice versa).

Example. To illustrate this procedure, we examine the Kepler catalog of Borucki et al. (2011), trimmed by 20% as described at the start of Section 4 to produce a more homogeneous set of target stars. This is catalog A. All of the planets in this catalog are detected with a signal-to-noise ratio (S/N) of at least 7. We construct a sequence of shallower catalogs (catalogs "B") by gradually increasing the minimum S/N up to values exceeding 100 at which point only a handful of multi-planet systems are left. The relation (12) implies that apart from statistical fluctuations the numbers of multiple-planet systems n^B_k, k = 1, ..., are functions only of f ^BA and the known $\mathbf {n}^{A}$, which approximates the observable multiplicity function $\overline{\mathbf {n}}^{A}$. Hence by eliminating f ^BA in favor of n^B₁, the number of k-planet systems in any survey B can be predicted as a function of the number of one-planet systems in that survey. These predictions for k = 2, 3, 4 are shown in the upper left panel of Figure 1 as solid lines, along with approximate estimates of the 1σ confidence bands (dashed lines, given by n_k ± n^1/2_k). The actual numbers of multi-planet systems after S/N cuts on the Kepler data are shown as open circles.

Zoom In Zoom Out Reset image size

To obtain a quantitative statistical test of the consistency between the predictions and the data, we use the Kolmogorov–Smirnov (K-S) statistic D ≡ max |n_k(predicted) − n_k(observed)|. Standard tables of the significance as a function of D are not directly applicable since n_k and n₁ are cumulative distributions of a third parameter, the S/N, rather than being directly related; however, the significance of a given value of D can be determined by Monte Carlo simulations. The p-values (probability of observing deviations D at least as extreme as those seen, given the null hypothesis) are 0.81, 0.52, and 0.15 for k = 2, 3, 4, indicating that the observed numbers of multi-planet systems are consistent with the separability approximation.

The upper right panel of Figure 1 shows a similar comparison for a sequence of catalogs based on cuts at increasing planet radius rather than S/N. The results are consistent with separability to within the statistical errors (p-values of 0.99, 0.43, and 0.83 for k = 2, 3, 4).

The lower panels of Figure 1 show similar results for RV surveys. The "A" catalog consists of 240 FGK dwarf stars hosting one or more planets (see Equation (46) for more detail), and the cuts are based on K (semi-amplitude of the RV curve) on the left and orbital period on the right. The predictions are consistent with separability for cuts based on K (p-values 0.47 and 0.11 for k = 2, 3) and for the cut based on period for k = 3 (p = 0.28), and marginally consistent for the cut based on period with k = 2 (p = 0.05).

Possibly the relatively poor performance of the separability approximation for cuts based on period reflects the observational strategy of most RV surveys: when observers find one planet around a star they tend to observe the star more frequently and for a longer interval, thereby enhancing the chance of finding a second planet even if its period is long or its mass is low.

These results do not establish that the separability approximation is correct but they do demonstrate that separability provides an effective tool for removing survey selection effects and converting the observable multiplicity function between surveys.

Throughout this paper we shall assume that tranets are in circular orbits. Moorhead et al. (2011) estimate that the mean eccentricity of planets discovered in the Kepler survey is only 0.1–0.25, so this assumption should not cause significant errors. We shall also assume that a transit occurs when the line of sight to the center of the planet intersects the stellar disk. This assumption should be approximately correct so long as the planetary radius is much smaller than the stellar radius (the median ratio of planetary radius to stellar radius in the Kepler survey is only 0.026).

Let R_⋆ be the radius of the star, a the semi-major axis of a planet in a circular orbit, and ≡ R_⋆/a. Consider a system containing n planets with semi-major axes specified by ₁, ..., _n. Let g_mn(₁, ..., _n) be the probability that a randomly oriented observer will detect m tranets in this system.

One planet. First consider the case n = 1. We define three unit vectors: $\hat{{\bf o}}$ points toward the observer, $\hat{{\bf n}}$ is normal to the planetary orbit, and $\hat{{\bf z}}$ is normal to the reference plane from which inclinations i are measured. Thus $\hat{{\bf z}}\,{\cdot}\, \hat{{\bf n}}=\cos i$ and $\hat{{\bf o}}\,{\cdot}\, \hat{{\bf n}}=\cos \gamma$. If the planet's size is negligible, it transits if and only if $|\hat{{\bf o}}\,{\cdot}\, \hat{{\bf n}}|<\epsilon$ or |cos γ| < so

$$\begin{equation} g_{11}(\epsilon)=1-g_{01}(\epsilon)=\frac{\int _{|\cos \gamma |<\epsilon }\sin \gamma \, d\gamma }{\int \sin \gamma \,d\gamma }=\epsilon. \end{equation} \tag{ 13 }$$

Two planets. Let h(w) = 1 if |w| < 1 and zero otherwise. Then transits occur if and only if h(⁻¹cos γ) is unity. We can expand this function as a series of Legendre polynomials P_ℓ(cos γ):

$$\begin{equation} h(\epsilon ^{-1}\cos \gamma)=\sum _{\ell =0}^\infty b_\ell (\epsilon) P_\ell (\cos \gamma). \end{equation} \tag{ 14 }$$

Multiplying this equation by P_n(cos γ) and integrating with respect to cos γ from −1 to 1, then using the orthogonality property ∫¹ _{− 1}dx P_ℓ(x)P_n(x) = 2δ_nℓ/(2ℓ + 1) we have

$$\begin{eqnarray} b_\ell (\epsilon)&=&\frac{2\ell +1}{2}\int _{-\epsilon }^\epsilon dx\,P_\ell (x)\nonumber\\ &=& \left\lbrace \begin{array}{@{}ll}\epsilon, & \ell =0, \\ P_{\ell +1}(\epsilon)-P_{\ell -1}(\epsilon), & \ell \hbox{ even, } \ell >0 \\ 0, & \ell \hbox{ odd.} \end{array}\right. \end{eqnarray} \tag{ 15 }$$

Now let (θ, ϕ) be the polar coordinates for $\hat{{\bf o}}$ relative to the polar axis $\hat{{\bf z}}$ and $(\Omega -\textstyle {\frac{1}{2}}\pi,i)$ be the polar coordinates for $\hat{{\bf n}}$. Then the addition theorem for spherical harmonics gives

$$\begin{equation} h(\epsilon ^{-1}\cos \gamma)=4\pi \sum _{\ell =0}^\infty \frac{b_\ell (\epsilon)}{2\ell +1}\sum _{m=-\ell }^\ell Y_{\ell m}^\ast (\theta,\phi)Y_{\ell m}(i,\Omega -\textstyle {\frac{1}{2}}\pi). \end{equation} \tag{ 16 }$$

Let the probability distribution of planetary inclinations be $q(i|\mbox{\boldmath $\kappa $})di$, where $\mbox{\boldmath $\kappa $}$ is a set of free parameters describing the inclination distribution, which we may vary to fit the observations. Then the probability of a transit of a single planet, given the observer orientation x ≡ cos θ, is

$$\begin{eqnarray} u(x|\epsilon,\mbox{\boldmath $\kappa $})&=&\int \frac{di\, d\Omega }{2\pi }q(i|\mbox{\boldmath $\kappa $})h(\epsilon ^{-1}\cos \gamma)\nonumber\\ &=& \displaystyle 4\pi \sum _{\ell =0}^\infty \frac{b_\ell (\epsilon)}{2\ell +1}\int di\,q(i|\mbox{\boldmath $\kappa $})Y_{\ell 0}^\ast (\theta,0)Y_{\ell 0}(i,0)\nonumber \\ &=& \displaystyle \sum _{\ell =0}^\infty Q_\ell (\mbox{\boldmath $\kappa $})\,b_\ell (\epsilon) P_\ell (x), \end{eqnarray} \tag{ 17 }$$

where

$$\begin{equation} Q_\ell (\mbox{\boldmath $\kappa $})\equiv \int _0^\pi di\,q(i|\mbox{\boldmath $\kappa $})P_\ell (\cos i), \quad Q_0=1. \end{equation} \tag{ 18 }$$

If a system contains two planets, the probability that both transit for a random orientation of the observer is

$$\begin{eqnarray} g_{22}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $})&=& {\frac{1}{2}}\int _{-1}^1 dx\, u(x|\epsilon _1,\mbox{\boldmath $\kappa $})u(x|\epsilon _2,\mbox{\boldmath $\kappa $}) \nonumber \\ &=&{\frac{1}{2}}\displaystyle \sum _{\ell,n=0}^\infty b_\ell (\epsilon _1)b_n(\epsilon _2)Q_\ell (\mbox{\boldmath $\kappa $}) Q_n(\mbox{\boldmath $\kappa $})\nonumber\\ &&\times\int _{-1}^1 dx\,P_\ell (x)P_n(x) \nonumber \\ &=& \displaystyle \sum _{\ell =0}^\infty \frac{Q_\ell ^2(\mbox{\boldmath $\kappa $})}{2\ell +1}b_\ell (\epsilon _1)b_\ell (\epsilon _2). \end{eqnarray} \tag{ 19 }$$

Moreover, the probability that one and only one of the two planets transits is

$$\begin{eqnarray} g_{12}(\epsilon _1,\epsilon _2,\kappa)&=& {\frac{1}{2}}\int _{-1}^1 dx\,\{u(x|\epsilon _1,\mbox{\boldmath $\kappa $})[1-u(x|\epsilon _2,\mbox{\boldmath $\kappa $})]\nonumber\\ &&{+}\,[1-u(x|\epsilon _1,\mbox{\boldmath $\kappa $})]u(x|\epsilon _2,\mbox{\boldmath $\kappa $})]\} \nonumber \\ &=&\displaystyle g_{11}(\epsilon _1,\mbox{\boldmath $\kappa $})+g_{11}(\epsilon _2,\mbox{\boldmath $\kappa $})-2g_{22}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $}) \nonumber\\ \end{eqnarray} \tag{ 20 }$$

and the probability that no planets transit is

$$\begin{eqnarray} g_{02}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $})&=& \displaystyle 1-g_{12}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $}) -g_{22}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $}) \nonumber \\ &=& \displaystyle 1-g_{11}(\epsilon _1,\mbox{\boldmath $\kappa $})-g_{11}(\epsilon _2,\mbox{\boldmath $\kappa $})+g_{22}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $}). \nonumber\\ \end{eqnarray} \tag{ 21 }$$

For example, if the planets are distributed isotropically then $q(i)di=\textstyle {\frac{1}{2}}\sin i\,di$, Q_ℓ = δ_ℓ0, and g₂₂(₁, ₂) = ₁₂. If the planets have zero inclination, it can be shown that

$$\begin{equation} g_{22}(\epsilon _1,\epsilon _2)=\sum _{\ell =0}^\infty \frac{b_\ell (\epsilon _1)b_\ell (\epsilon _2)}{2\ell +1}=\hbox{min\,}(\epsilon _1,\epsilon _2), \end{equation} \tag{ 22 }$$

although this result is derived more easily in other ways.

Three or more planets. These results can be extended to any number of planets³:

$$\begin{eqnarray} g_{mn}(\epsilon _1,\ldots,\epsilon _n,\mbox{\boldmath $\kappa $})&=&\textstyle {\frac{1}{2}}\int _{-1}^1 dx\, \sum _{P_n} \prod _{i=1}^m u(x|\epsilon _{p_i},\mbox{\boldmath $\kappa $})\nonumber\\ &&\times\prod _{j=m+1}^n [1-u(x|\epsilon _{p_j},\mbox{\boldmath $\kappa $})], \end{eqnarray} \tag{ 23 }$$

where P_n is the set of all permutations (p₁, ..., p_n) of the numbers 1, ..., n, and m ⩽ n. For example,

$$\begin{eqnarray} g_{23}(\epsilon _1,\epsilon _2,\epsilon _3,\mbox{\boldmath $\kappa $})&=& g_{22}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $})+ g_{22}(\epsilon _2,\epsilon _3,\mbox{\boldmath $\kappa $})\nonumber\\ &&+g_{22}(\epsilon _3,\epsilon _1,\mbox{\boldmath $\kappa $}) - 3g_{33}(\epsilon _1,\epsilon _2,\epsilon _3,\mbox{\boldmath $\kappa $}) \nonumber \\ g_{13}(\epsilon _1,\epsilon _2,\epsilon _3,\mbox{\boldmath $\kappa $}) &=& g_{11}(\epsilon _1,\mbox{\boldmath $\kappa $}) + g_{11}(\epsilon _2,\mbox{\boldmath $\kappa $}) + g_{11}(\epsilon _3,\mbox{\boldmath $\kappa $})\nonumber\\ &&- 2g_{22}(\epsilon _1,\epsilon _2,\mbox{\boldmath $\kappa $})-2g_{22}(\epsilon _2,\epsilon _3,\mbox{\boldmath $\kappa $}) \nonumber \\ && -2g_{22}(\epsilon _3,\epsilon _1,\mbox{\boldmath $\kappa $}) +3g_{33}(\epsilon _1,\epsilon _2,\epsilon _3,\mbox{\boldmath $\kappa $}). \nonumber\\ \end{eqnarray} \tag{ 24 }$$

The geometric selection matrix G_mn (Equation (1)) is simply $\langle g_{mn}(R_\star /a_1,R_\star /a_2,\ldots,R_\star /a_l,\mbox{\boldmath $\kappa $})\rangle$, the average of the geometric selection factor over the joint distribution of stellar radius R_⋆ and planetary semi-major axis a for the survey. To evaluate $G_{mn}(\mbox{\boldmath $\kappa $}),$ we use the separability approximation (4) with respect to = R_⋆/a. Thus,

$$\begin{equation} G_{mn}(\mbox{\boldmath $\kappa $})=\int g_{mn}(\epsilon _1,\ldots,\epsilon _n,\mbox{\boldmath $\kappa $})\prod _{k=1}^n f(\epsilon _k)d\log \epsilon _k, \end{equation} \tag{ 25 }$$

where f()dlog  represents the probability distribution of (i.e., the semi-major axis distribution) as modified by the survey selection effects.

With this parameterization and Equations (17) and (23) it is straightforward to show that $G_{mn}(\mbox{\boldmath $\kappa $})$ is given by the binomial distribution,

$$\begin{eqnarray} G_{mn}(\mbox{\boldmath $\kappa $})&=&\frac{n!}{2 m!(n-m)!}\int _{-1}^1 dx\, U^m(x|\mbox{\boldmath $\kappa $})[1-U(x|\mbox{\boldmath $\kappa $})]^{n-m}\nonumber\\ &=&{\frac{1}{2}}\int _{-1}^1 dx\,S_{mn}[U(x|\mbox{\boldmath $\kappa $})], \end{eqnarray} \tag{ 26 }$$

where S_mn is given by Equation (7),

$$\begin{equation} U(x|\mbox{\boldmath $\kappa $})\equiv \frac{\int f(\epsilon)u(x|\epsilon,\mbox{\boldmath $\kappa $})\,d\log \epsilon }{\int f(\epsilon)\,d\log \epsilon }=\sum _{\ell =0}^\infty Q_\ell (\mbox{\boldmath $\kappa $}) B_\ell P_\ell (x), \end{equation} \tag{ 27 }$$

and

$$\begin{equation} B_\ell \equiv \int f(\epsilon) b_\ell (\epsilon)\,d\log \epsilon \quad \mbox{with}\quad \int f(\epsilon)\,d\log \epsilon =1. \end{equation} \tag{ 28 }$$

Since B_ℓ does not depend on the unknown parameters $\mbox{\boldmath $\kappa $}$ of the inclination distribution, it can be evaluated once and for all at the start of any optimization procedure. It is straightforward to show that the relations (2) are satisfied by these formulae, and that the matrices $\mathbf {G}$ and $\mathbf {S}$ commute. In numerical work, we typically truncate infinite series such as Equation (27) at ℓ = ℓ_max = 50, but for very thin disks it may be necessary to include terms of higher ℓ.

We pointed out in Equation (10) that most survey selection effects preserve the separability approximation. This result does not generally hold for geometric selection effects. To illustrate this, consider the simple case of a population of stars containing two planets, with zero relative inclination. Write the probability distribution of = R_⋆/a of two-planet systems as f(₁)f(₂)dlog ₁dlog ₂ (after survey selection effects but before geometric selection effects). Then using Equation (22) it is evident that the probability distribution of two-tranet systems is

$$\begin{equation} dp_2(\epsilon _1,\epsilon _2)=f(\epsilon _1)f(\epsilon _2)\hbox{min\,}(\epsilon _1,\epsilon _2)\,d\log \epsilon _1\,d\log \epsilon _2, \end{equation} \tag{ 29 }$$

which is not separable. Nevertheless, in a wide range of systems, geometric selection effects do approximately conserve separability, if the rms inclination is not too small and the planets are not too close to the host star. A specific example is given in Section 3.2.

3.1. The Inclination Distribution

In this paper, we model the probability distribution of the inclinations dp = q(i|κ)di as a Fisher distribution,

$$\begin{equation} q(i|\kappa)=\frac{\kappa }{2\sinh \kappa }\exp (\kappa \cos i)\sin i. \end{equation} \tag{ 30 }$$

The parameter κ is related to the mean-square value of sin i through

$$\begin{equation} \langle \sin ^2i\rangle =\int di\sin ^2i\,q(i|\kappa)= 2\frac{\coth \kappa }{\kappa }-\frac{2}{\kappa ^2}. \end{equation} \tag{ 31 }$$

When κ ≪ 1, the Fisher distribution approaches an isotropic distribution, $\lim _{\kappa \rightarrow 0}q(i|\kappa)=\frac{1}{2}\sin i$, while for κ ≫ 1 it approaches the Rayleigh distribution, $\lim _{\kappa \rightarrow \infty }q(i|\kappa)=(2i/s^2)\*\exp (-i^2/s^2),$ where s = (2/κ)^1/2 is the rms inclination and $\frac{1}{2}\pi ^{1/2}s=0.8862s$ is the mean inclination. The Rayleigh distribution is commonly used to model the inclination distribution of asteroids, Kuiper Belt objects, stars in the Galactic disk (where it is known as the Schwarzschild distribution), etc. As κ → −∞, the Fisher distribution approaches a retrograde Rayleigh distribution.

For the Fisher distribution, Equation (18) becomes

$$\begin{equation} Q_\ell (\kappa)=\sqrt{\frac{\pi \kappa }{2}}\frac{I_{\ell +1/2}(\kappa)}{\sinh \kappa }, \end{equation} \tag{ 32 }$$

where I denotes a modified Bessel function.

3.2. Validity of the Separability Approximation

The separability approximation described in Section 1.2 can only be approximately valid—for example, it is inconsistent with the observational finding that planets tend to be concentrated near mutual orbital resonances, and with the theoretical finding that planets separated by less than a few Hill radii are unstable. Nevertheless, its success in accounting for survey selection effects in the example at the end of Section 2 suggests that it may be sufficiently accurate for a preliminary analysis of the statistics of multi-planet systems.

There is limited direct evidence on the accuracy of the separability approximation for multi-planet systems. First, consider RV surveys, in which there are no geometric selection effects. The most important survey selection effects depend only on the properties of an individual planet so an RV survey of a separable parent distribution should lead to a separable detected distribution (Equation (10)).

The separability approximation implies that the semi-major axis distributions in single- and multi-planet systems should be the same. These distributions are shown in the bottom left panel of Figure 2; as usual we restrict ourselves to FGK dwarfs. The distributions look quite similar, and a K-S test confirms this impression (p = 0.49); thus, the semi-major axis distribution appears to be consistent with separability. Wright et al. (2009) reached a different conclusion, probably because they included ∼8 planets discovered by transit surveys in their sample (J. Wright 2011, private communication).

Zoom In Zoom Out Reset image size

Separability also implies that the mass distributions in single- and multi-planet systems should be the same. These distributions are shown in the top right panel of Figure 2. The distributions look different—there is an excess of low-mass planets in multi-planet systems—and this is confirmed by a K-S test (p = 0.01). Wright et al. (2009) point out that this difference may be amplified by an unmodeled survey selection effect—stars hosting planets tend to be observed more frequently, thereby enhancing the chance of discovering additional low-mass planets. To check this possibility we can restrict the sample to planets that are easily detectable, with velocity semi-amplitude K > 5 m s⁻¹. With this restriction, the single- and multi-planet mass distributions are statistically indistinguishable (p = 0.47).

The eccentricity distributions in single- and multi-planet systems are clearly different (upper left panel of Figure 2): multi-planet systems have smaller eccentricities. The mean eccentricity is 0.28 in single-planet systems and only 0.20 in multi-planet systems. This result is not surprising since orbital stability favors low-eccentricity orbits when more than one planet is present. Thus, the eccentricity distribution is not consistent with separability, although this has little or no effect on the arguments in this paper since eccentricity does not strongly affect the observability of planets in either transit or RV surveys.

Note that separability does not require that the properties of the host stars in single- and multi-planet systems should be the same; thus, possible differences in metallicity or mass of the host star between single- and multi-planet systems are not tests of separability (see Wright et al. 2009 for discussion).

The evidence on separability from the Kepler survey is more difficult to interpret because geometric selection effects do not preserve separability (see discussion just before Section 3.1). Nevertheless, in a wide range of possible systems, separability implies that the semi-major axis distributions of single- and multiple-tranet systems in the Kepler survey should be quite similar. To illustrate this, the lower right panel of Figure 2 shows the probability distributions dp = p(log )dlog  of = R_⋆/a for the planets in the Kepler sample, with single-tranet systems in black and multi-tranet systems in red. The two distributions are almost the same according to a K-S test (p = 0.13). If the single-tranet systems arise mostly from single-planet systems their probability distribution should be dp = f() dlog , where f() is given by Equation (35)—the additional factor of accounts for geometric selection effects (Equation (13))—and this curve, plotted in green, offers a good fit to the data. If the multiple-tranet systems arise mostly from two-planet systems we can also calculate their expected distribution; this depends on the rms inclination and ranges from the green curve if the inclination distribution is isotropic to the magenta curve for flat systems (cf. Equation (29)). The distribution of observed multi-tranet systems in red is roughly consistent with either curve, implying that the semi-major axis distribution of multi-tranet systems in the Kepler survey is consistent with separability and insensitive to the rms inclination.

Latham et al. (2011) have shown that Kepler systems with multiple tranets are less likely to include a giant planet (larger than Neptune) than systems with a single tranet. We confirm using a K-S test that the distributions of radii in the single- and multiple-tranet systems are different at a high significance level. This difference is similar to the one described above for RV surveys—more low-mass planets in multiple systems. Its interpretation is not straightforward since geometric selection effects can lead to differences of this kind if the masses and semi-major axes of planets are correlated.

These comparisons show that separability can successfully account for some but not all of the properties of multi-planet systems in the RV and Kepler planet samples. Further work, and better data, are needed to assess the role of this approximation in describing the properties of multi-planet systems. Nevertheless, these preliminary comparisons are encouraging enough that in the remainder of this paper we shall use separability to explore the inclination distribution in exoplanet systems.

4.1. Properties of the Survey

The Kepler survey has a complex set of survey selection effects, which we do not attempt to model. The constraints on the multiplicity function that we derive therefore apply to the population of planets in radius, semi-major axis, etc., that Kepler detects, whatever that population may be (for a discussion of selection effects and completeness in the Kepler catalog see Howard et al. 2011 and Youdin 2011). If we denote the multiplicity function of this population by $\mathbf {N}$ and the observable multiplicity function of the Kepler survey by $\overline{\mathbf {n}}$ then Equation (1) becomes

$$\begin{equation} \overline{\mathbf {n}}=\mathbf {G}\cdot \mathbf {N}. \end{equation} \tag{ 33 }$$

The validity of this equation requires only the plausible approximation that the probability that Kepler will detect a given transiting planet around a given star is independent of whether it detects other transits around the same star.

To produce a more homogeneous sample, we trim the catalog of Borucki et al. (2011) to include only stars with effective temperatures between 4000 and 6500 K and surface gravity log g > 4.0 (roughly equivalent to FGK dwarfs), and to Kepler magnitudes between 9.0 and 16.0; this trimming leaves 124,613 stars from the original sample of 153,196. We also restrict the catalog to planets with orbital period less than 200 days and radius less than 2 Jupiter radii; this leaves 1092 planets from the original sample of 1235. The numbers of stars with 0, 1, 2, ... tranets in the trimmed catalog are

$$\begin{equation} \begin{array}{@{}l@{\ \ }l@{\ \ }l@{\ \ }l@{}} n_0=1.237\times 10^5, & n_1=737, & n_2=104, & n_3=37, \\ n_4=7, & n_5=1, & n_6=1, & n_k=0\hbox{ for }k>6.\end{array} \end{equation} \tag{ 34 }$$

We need to determine the function f(), where f() dlog  is the fraction of planets in the range dlog  given the intrinsic distribution of planets and the survey selection effects for Kepler. As usual = R_⋆/a is the ratio of stellar radius to planetary semi-major axis; the stellar radius is determined from the host-star mass and surface gravity and the semi-major axis is determined from the host-star mass and the planetary orbital period. Figure 3 shows data points for f(), constructed by adding a contribution of ⁻¹ (to account for geometric selection effects) from each tranet to the corresponding bin, then normalizing so that the integral over log  is unity. The distribution can be adequately fit by the parameterization

$$\begin{equation} f(\epsilon) = 0.656\frac{(\epsilon /\epsilon _0)^{0.5}}{1+(\epsilon /\epsilon _0)^{3.6}}, \quad \epsilon _0=0.055,\quad \mbox{for }\epsilon >0.004 \end{equation} \tag{ 35 }$$

and zero for < 0.004. The sharp decline for ≳ 0.1 is due to an absence of planets with semi-major axis ≲ 0.04 AU (Borucki et al. 2011), while the cutoff at ≲ 0.004 is due to the limited timespan of the Kepler data.

Zoom In Zoom Out Reset image size

4.2. Statistical Method

The probability that the survey actually detects {n₀, n₁, ..., n_K} stars having 0, 1, ..., K planets is

$$\begin{equation} P(\mathbf {n}|\mathbf {N},\mbox{\boldmath $\kappa $})=\prod _{k=0}^K \frac{\overline{n}_k^{n_k}\exp (-\overline{n}_k)}{n_k!}, \end{equation} \tag{ 36 }$$

where $\mathbf {n}=(n_0,n_1,\ldots,n_K)$ and $\overline{n}_k$ is related to $\mathbf {N}$ by Equation (33).

Estimating the multiplicity function $\mathbf {N}$ and the inclination distribution parameters $\mbox{\boldmath $\kappa $}$ from $\mathbf {n}$ is a straightforward but challenging problem in statistics and optimization. This problem can be attacked with a variety of methods (linear programming, minimum χ², maximum likelihood, Bayesian analysis using a Markov Chain Monte Carlo algorithm, etc.), and we have experimented with several of these. In this paper, we have usually chosen maximum likelihood, as a reasonable compromise between generality, computation time, and clarity of interpretation.

The log of the likelihood of a given observational result $\mathbf {n}$ is

$$\begin{eqnarray} \log P(\mathbf {n}|\mathbf {N},\mbox{\boldmath $\kappa $})&=&\sum _{k=0}^K n_k\log \left[\sum _{l=k}^K G_{kl}(\mbox{\boldmath $\kappa $})N_l\right]\nonumber\\ &&-\sum _{k=0}^K\sum _{l=k}^K G_{kl}(\mbox{\boldmath $\kappa $})N_l-\sum _{k=0}^K\log n_k!. \end{eqnarray} \tag{ 37 }$$

Note that the second term on the right can be simplified to ∑_lN_l using Equation (3). We then maximize log P with respect to $\mathbf {N}$ and $\mbox{\boldmath $\kappa $}$, subject to the constraints N_k ⩾ 0, k = 0, ..., K.

4.3. Results

The top panel of Figure 4 shows the maximum likelihood as a function of the rms inclination and the maximum number of planets per system, K, for 6 ⩽ K ⩽ 40. The minimum allowed value of K is 6, since Kepler has found one system with six tranets. The maximum-likelihood models with a given K are connected to form solid lines, and the families with K = 10, 20, 30, and 40 are colored for emphasis. There are occasional small dips in the lines when the optimization algorithm (a quasi-Newton algorithm from the Numerical Algorithms Group) converged on a local rather than global maximum. The figure shows that:

• 1.  
  The highest likelihood is for razor-thin systems, with near zero rms inclination. However, the preference for zero rms inclination has only marginal statistical significance: systems exist at all rms inclinations—even isotropic systems—with likelihood only exp (− 0.73) smaller than the razor-thin solutions.
• 2.  
  Systems with large rms inclinations are only consistent with the data if a fraction of them contain a large number of planets. At the 3σ level (likelihood smaller than the maximum by exp (− 4.5), marked by a horizontal dashed line in the figure), the maximum rms inclination is related to the maximum number of planets by

  $$\begin{equation} \langle \sin ^2i\rangle ^{1/2}\le \left\lbrace \begin{array}{@{}ll}0.15+0.037(K-6), & K<24 \\ (\frac{2}{3})^{1/2}\mbox{ (isotropic)}, & K\ge 24.\end{array}\right. \end{equation} \tag{ 38 }$$

Zoom In Zoom Out Reset image size

It is possible, of course, that even the maximum-likelihood model does not fit the data well. To explore this possibility, we have calculated the standard Pearson χ² statistic,

$$\begin{equation} \chi ^2=\sum _{k=0}^K\frac{(n_k-\overline{n}_k)^2}{\overline{n}_k}=\sum _{k=0}^K\frac{(n_k-\sum _l G_{kl}N_l)^2}{\sum _l G_{kl}N_l}. \end{equation} \tag{ 39 }$$

The distribution of the χ² statistic is not straightforward to interpret, since $\overline{n}_k\lesssim 1$ for many k and since the number of degrees of freedom is not well defined. Nevertheless, it is probably reasonable to expect that there is a good fit to the data if χ² ≲ 5. The values of χ² for the maximum-likelihood solutions in the top panel of Figure 4 are shown in the bottom panel of that figure. There are satisfactory models with all rms inclinations, but as before such models require that some systems contain many planets if the rms inclination is large.

It is instructive to examine the isotropic solution with K = 30 in more detail (the behavior of the isotropic solutions with K > 30 is qualitatively similar). The fraction of stars with k-planet systems is

$$\begin{equation} \frac{N_k}{\sum _{l=0}^K N_l}=\left\lbrace \begin{array}{@{}ll}0.944 & k=0, \\ 0.0065 & k=1, \\ 0 & k=2, \\ 0.0452 & k=3, \\ 0 & k=4,\ldots,29, \\ 0.0043 & k=30. \end{array}\right. \end{equation} \tag{ 40 }$$

Thus, in this solution, about half of the planets are contained in three-planet systems, and the other half in a small population (<0.5%) of stars with many-planet systems. This multiplicity function and inclination distribution are neither unique nor particularly plausible but they are consistent with the Kepler data.

Figure 5 shows the fraction of stars in the Kepler sample with 0-, 1-, 2-, 3-, ..., planet systems, as a function of the assumed rms inclination. The results are for K = 30 but are qualitatively similar for larger values of K. Our initial attempts to construct this figure were unsuccessful because the appearance of the figure is very sensitive to cases when the optimization algorithm settles on a local maximum of the likelihood. To avoid this difficulty, we re-cast the optimization as a problem in linear programming: we demanded that each $\overline{n}_k$ should lie within the 90% confidence interval determined by the Poisson distribution (36), and from these solutions we chose the one with the minimum total number of planets ∑^K_{k = 1}N_k. This specifies a unique solution, if one exists.

Zoom In Zoom Out Reset image size

At the smallest inclinations (〈sin ²i〉^1/2 < 0.05), the solution contains a mix of 1-, 2-, 3-, 4-, and 8- or 9-planet systems. As the rms inclination increases, the mixture becomes strongly dominated by 1-planet and n-planet systems where n varies monotonically with the rms inclination—for example, n = 12 when 〈sin ²i〉^1/2 ≃ 0.3. We caution that these results should not be regarded as a prediction of the Kepler multiplicity function for a given rms inclination because they are likely to depend on our arbitrary choice of cost function.

The need for many-planet systems is straightforward to understand. Consider the extreme case of an isotropic distribution. Then κ = 0 and q(i|κ = 0) = 1/2sin i; thus, Q_ℓ(κ = 0) = δ_ℓ0 from Equation (18) and the orthogonality properties of the Legendre polynomials. Thus, U(x|κ = 0) = B₀ (Equation (27)) and using Equation (26)

$$\begin{equation} G_{mn}(\mbox{\boldmath $\kappa $})=\frac{n!}{m!(n-m)!}B_0^m(1-B_0)^{n-m}. \end{equation} \tag{ 41 }$$

If all systems contain n planets, the ratio of the number of m-tranet systems to the number of (m + 1)-tranet systems is

$$\begin{equation} \frac{G_{mn}}{G_{m+1,n}}=\frac{m+1}{n-m}\frac{1-B_0}{B_0}, \quad n\ge m+1. \end{equation} \tag{ 42 }$$

Using Equations (28) and (35) we find that B₀ = 0.0321 for the Kepler survey. From Equation (34) we find n₁/n₂ = 7.1 ± 0.7. For comparison, the ratio G_1n/G_2n is less than 7.1 + 0.7 = 7.8 only for n ⩾ 9; thus, any population dominated by systems with less than nine planets will overproduce 1-tranet systems relative to 2-tranet systems. Similarly, for the Kepler survey n₂/n₃ = 2.8 ± 0.5, and G_2n/G_3n > 2.8 + 0.5 = 3.3 unless n ⩾ 30.

The average number of planets per star from these solutions is shown in Figure 6. This result is insensitive to the rms inclination and the maximum number of planets per star (K), since it is given simply by the ratio of the total number of planets to the number of target stars, divided by the probability that a single randomly oriented planet will transit (Youdin 2011). Mathematically,

$$\begin{equation} \langle \mbox{number of planets per star}\rangle =\frac{\sum _{k=1}^K kn_k}{B_0\sum _{k=0}^Kn_k}=0.274. \end{equation} \tag{ 43 }$$

Zoom In Zoom Out Reset image size

The large open circles in Figure 6 show the probability that a system with one, two, or three tranets has additional planets. Typically the fraction of one-tranet systems with additional planets is 0.2–0.5, without a strong dependence on rms inclination. For two or three tranets the probability that there are additional unseen planets is substantially higher. The additional planets may be detectable by transit timing variations (Ford et al. 2011).

The observable multiplicity function in transit surveys depends strongly on the inclination distribution, but the observable multiplicity function in RV surveys does not. Thus, a comparison of the observable multiplicity functions of transit and RV surveys can offer a powerful probe of the inclination distribution. The principal obstacle to making this comparison is that the masses and orbital periods of the planets detected through these two observational techniques are quite different, as illustrated in Figure 7, and the multiplicity functions in these two regions of parameter space are likely to be different. In this section, we use the separability approximation and the methods of Section 2 to overcome this obstacle.

Zoom In Zoom Out Reset image size

Suppose that we wish to combine the Kepler survey with a RV survey (or a set of such surveys). The surveys yield n^Kep_k and n^RV_k systems containing k planets. We assume that both surveys have similar target star populations (we cull the list of target stars in both cases to include only FGK dwarfs), with multiplicity function $\mathbf {N}$ for Kepler and $c\,\mathbf {N}$ for the RV survey, where c < 1 is a constant to be determined. Let $\mathbf {S}(W^\mathrm{Kep})$ and $\mathbf {S}(W^\mathrm{RV})$ be the survey selection functions. We assume that there are no geometric selection effects for the RV surveys (cf. footnote ³), and the geometric selection effects for the Kepler survey are represented by the matrix $\mathbf {G}(\mbox{\boldmath $\kappa $})$ (Equation (26)). The generalization of Equation (36) for the likelihood is

$$\begin{eqnarray} P(\mathbf {n}^\mathrm{Kep},\mathbf {n}^\mathrm{RV}|\mathbf {N},\mbox{\boldmath $\kappa $})&=&\prod _{k=0}^K\frac{\big(\overline{n}_k^\mathrm{Kep}\big)^{n_k^\mathrm{Kep}}\exp \big(-\overline{n}_k^\mathrm{Kep}\big)}{ n_k^\mathrm{Kep}!}\nonumber\\ &&\times\prod _{k=1}^K\frac{\big(\overline{n}_k^\mathrm{RV}\big)^{n_k^\mathrm{RV}}\exp \big(-\overline{n}_k^\mathrm{RV}\big)}{n_k^\mathrm{RV}!}, \end{eqnarray} \tag{ 44 }$$

where

$$\begin{equation} \overline{\mathbf {n}}^\mathrm{Kep}=\mathbf {G}(\mbox{\boldmath $\kappa $})\mathbf {S}(W^\mathrm{Kep})\mathbf {N}, \quad \overline{\mathbf {n}}^\mathrm{RV}=c\,\mathbf {S}(W^\mathrm{RV})\mathbf {N}. \end{equation} \tag{ 45 }$$

Note that the second product in Equation (44) starts at k = 1 since it is difficult to determine how many stars have been unsuccessfully examined for planets by RV methods (see further discussion below). We then maximize the likelihood (44) over N₀, N₁, ..., N_K, W^Kep, W^RV, and c (as shown in Section 2, the likelihood actually depends only on the ratio W^RV/W^Kep).

We determine the observable multiplicity function for RV planets using all planets with FGK dwarf host stars in the exoplanets.org database (Wright et al. 2011) as of 2010 August,

$$\begin{eqnarray} n_1^\mathrm{RV}&=&162, \quad n_2^\mathrm{RV}=24, \quad n_3^\mathrm{RV}=7, \quad n_4^\mathrm{RV}=1,\nonumber\\ n_5^\mathrm{RV}&=&1, \quad n_k^\mathrm{RV}=0\hbox{ for }k>5, \end{eqnarray} \tag{ 46 }$$

for a total of 240 planets. The observable multiplicity function for Kepler planets is given in Equation (34). Figure 8 shows the maximum likelihood as a function of the rms inclination and the maximum number of planets per system, K (top), as well as χ² for these models (bottom). The plots are noisier than Figure 4, presumably because the optimization algorithm was less successful at finding the global maximum likelihood, but otherwise look similar. In particular, systems with large rms inclinations are consistent with the data if and only if a small fraction of them contain a large number of planets. Evidently adding data from RV surveys has not significantly tightened the constraints on the inclination distribution.

Zoom In Zoom Out Reset image size

We now show that adding information on the total number of target stars in the RV surveys does allow the inclination distribution to be determined. Figure 9 shows the expected numbers $\overline{n}_k^\mathrm{Kep}$ and $\overline{n}_k^\mathrm{RV}$ of k-tranet systems from the Kepler survey and k-planet systems from the RV surveys, as determined from the maximum-likelihood solutions described above. Each point corresponds to a given maximum number of planets (6 ⩽ K ⩽ 40) and rms inclination, and only solutions within 3σ of the global maximum likelihood are shown. The points with error bars, surrounded by circles for greater visibility, correspond to the observed numbers n^Kep_k and n^RV_k from Equations (34) and (46). Most of the expected values lie within the error bars of the corresponding observed value; this is no more than a confirmation that our optimization code is performing properly. The blue points show the total number of stars in the RV surveys, $\overline{n}_\mathrm{tot}^\mathrm{RV}=\sum _{k=0}^K \overline{n}_k^\mathrm{RV}$, as determined by the optimization code. The plot shows that $\overline{n}_\mathrm{tot}^\mathrm{RV}$ is tightly correlated with the rms inclination, so an accurate characterization of the total number of RV target stars would enable the determination of the rms inclination.

Zoom In Zoom Out Reset image size

This task is challenging given the heterogeneous surveys that have produced the RV planets known at the present time. We have used the following two distinct approaches.

• 1.  
  Cumming et al. (2008) have carefully examined selection effects in the Keck Planet Search, and derive the percentage of F, G, and K stars with a planet in various ranges of orbital period and mass. The sample of RV planets used in our analysis (Equation (46)) is not corrected for selection effects, but for sufficiently massive planets and sufficiently short orbital periods it should be complete. For example, for planets more massive than Jupiter, Msin γ > M_J, with orbital periods less than one year, P < 1 yr, the velocity semi-amplitude K_RV > 30 m s⁻¹ which is large enough to be detectable in most surveys. In this mass and period range our sample contains 46 planet-hosting stars and Cumming et al. (2008) estimate that the fraction of stars with planets is 0.019 ± 0.007, which implies n^RV_tot = 2400 ± 900. Decreasing the period range to P < 100 days gives n^RV_tot = 2500 ± 1200 (based on 21 host stars); lowering the mass cutoff to Msin γ > 0.5 M_J gives n^RV_tot = 1900 ± 500 (based on 63 host stars). This last estimate of n^RV_tot is probably low because the surveys we have used are not all complete at this level.
• 2.  
  We may estimate n^RV_tot using the tranet frequency derived from the Kepler mission. Once again, we restrict the Kepler sample to host stars that are F, G, and K dwarfs (4000 K <T_eff < 6500 K and log g > 4). We then carry out the following steps for a given orbital period P and velocity semi-amplitude K: (1) compute the corresponding mass M(P, K) = M_J(K/30 m s⁻¹)(1 yr/P)^1/3 assuming a circular orbit and a solar-mass host star; (2) find the number n^RV(P, K) of RV planets with period less than P and mass greater than M(P,K); (3) find all Kepler tranets in the same period and mass range, using an empirical mass–radius relation found by fitting mass and radius measurements from transiting planets in the range 0.1–10M_J (see Figure 10) to a log–quadratic relation

  $$\begin{eqnarray} \log R/R_J&=&0.087 +0.141\log M/M_{J}\nonumber\\ &&-0.171\left(\log M/M_{J}\right)^2; \end{eqnarray} \tag{ 47 }$$

  (4) compute the total number of Kepler planets in this range n^Kep(P, K) by counting each tranet as ⁻¹ planets to correct for geometric selection effects (Equation (13)); (5) estimate the total number of RV host stars as n^RV_tot = n_tot^Kepn^RV(P, K)/n^Kep(P, K). The results are shown in Figure 11 for K = 10, 15, 20, and 25 m s⁻¹. As the majority of RV surveys have reached precisions of ∼15 m s⁻¹ or better over the last decade, it is reassuring but not surprising that the estimates of n^RV_tot for K = 15, 20, and 25 m s⁻¹ are consistent. The rise in n^RV_tot at small periods is likely due to the known discrepancy in hot Jupiter frequency between transit and RV surveys (the frequency of hot Jupiters estimated from transit surveys is factor of ∼2 smaller than that derived from RV surveys, perhaps because the average metallicities are different; see Gould et al. 2006 and Howard et al. 2011).

Zoom In Zoom Out Reset image size

These independent approaches yield n^RV_tot ≃ 2500 ± 1000 and n^RV_tot ≃ 3000 ± 1000, respectively, which are consistent within the errors. The corresponding inclination ranges from Figure 9 are 0 < 〈sin ²i〉^1/2 < 0.08 and 0.02 < 〈sin ²i〉^1/2 < 0.09 which correspond to an rms or mean inclination range of 0°–5° (as shown in Section 3.1, for a Rayleigh distribution the rms inclination is only larger than the mean inclination by 12%, which is much less than the uncertainty).

If the separability approximation is valid, our results should be insensitive to cuts made on the Kepler planet candidates. To check this, we have repeated the analysis for the Kepler sample examined by Lissauer et al. (2011b), who imposed a period cut 3 days < P < 125 days, a radius cut 1.5 R_⊕ ⩽ R ⩽ 6 R_⊕, and a signal/noise cut S/N ⩾ 16, which reduced the number of planets by almost 40%. We find the mean inclination for this sample to be 0°–4°, not significantly different from the estimate in the preceding paragraph.

Although the range of rms inclinations is tightly constrained by this analysis, the multiplicity function is not. For example, within 1σ of the maximum-likelihood model (Δlog P ⩽ 0.5) we have found models that have no 1-planet systems (67% have no planets, 29% have 2 planets, and 4% have 13 planets) and others that have no zero-planet systems (93% have 1 planet, 2% have 6 planets, and 5% have 25 planets).

A by-product of this analysis is the ratio W^RV/W^Kep (Equation (45)), which measures the relative sensitivity of the RV and Kepler surveys. This ratio varies smoothly from 0.5 for razor-thin systems to 0.2 for 〈sin ²i〉^1/2 = 0.1, independent of the maximum number of planets in the model. In other words, 20%–50% of the Kepler planets could have been detected in RV surveys. If this ratio can be determined independently by fitting models of the period, radius, and mass distributions it will provide a constraint on the rms inclination that does not require estimating the total number of RV target stars.

A weak link in these arguments is the assumption that the population of FGK dwarf stars is the same in the Kepler and RV surveys. However, we note that our two approaches to estimating n^RV_tot, one using only RV surveys and one comparing the Kepler and RV surveys, yield similar answers, which suggests that the estimate of the rms inclination that we derive from this answer is insensitive to differences between the host stars of the Kepler and RV surveys.

It is interesting to compare this estimate of the mean inclination to the mean eccentricity for Kepler planets. Restricting our sample to planets with period P > 10 days (to avoid the effects of tidal circularization) and minimum mass between 0.01 and 0.1 Jupiter masses, the mean eccentricity of planets discovered in RV surveys is 0.15 (we have also excluded planets with a reported eccentricity of zero, which may include cases in which no eccentricity was fit). These results are roughly consistent with estimates of the mean eccentricity of Kepler planets from transit timing: Moorhead et al. (2011) find that the mean eccentricity is between 0.13 and 0.25 at a p-value of 0.05. We have

$$\begin{equation} \frac{\langle i\rangle }{\langle e\rangle }=0.35\frac{\langle i\rangle }{3^\circ }\frac{0.15}{\langle e\rangle }. \end{equation} \tag{ 48 }$$

Theoretical studies of eccentricity and inclination growth in planetesimal disks (e.g., Ida et al. 1993) find 〈i〉/〈e〉 = 0.45–0.5, somewhat larger than this value. A possible explanation is that the eccentricities may have been systematically overestimated. Zakamska et al. (2011) find that the typical bias due to measurement errors is Δe ∼ 0.04 in RV catalogs, and the bias in the present sample is likely to be higher since the S/N is low for low-mass planets. Possibly a similar bias is present in the Kepler measurements of the eccentricity distribution.

The Kepler survey can measure transit timing variations of a minute or less in favorable cases (Ford et al. 2011). These variations can be used to detect and characterize additional planets. Given the rms inclination of 0–0.09 radians that we have derived, roughly 20%–30% of the single-tranet Kepler systems are expected to have additional planets (Figure 6), and many of these may be detectable by transit timing variations. Ford et al. (2011) estimate that ∼10%–20% of suitable Kepler tranets show evidence of transit timing variations, and this number is likely to increase as the survey duration grows. Figure 6 also shows that the fraction of two- or three-tranet systems with additional planets is substantially higher and strongly dependent on the rms inclination. A preliminary analysis by Ford et al. (2011) estimates that the fraction of two- and three-tranet systems showing transit timing variations is similar to the fraction for 1-tranet systems, 10%–20%. This result is difficult to reconcile with any of our models, whatever the rms inclination may be.

Zoom In Zoom Out Reset image size

We have described a methodology for analyzing the multiplicity function—the fraction of host stars containing a given number of planets—in RV and transit surveys. Our approach is based on the approximation of separability that the probability distribution of planetary parameters in an n-planet system is the product of identical one-planet distributions (Section 1.2). Exoplanet surveys show that separability is not precisely satisfied but the departures from this approximation are small enough that it provides a powerful tool for the study of multi-planet systems. Using this approximation, we have shown how to relate the observable multiplicity function in surveys with different sensitivities, so long as they examine populations of potential host stars with similar properties (Section 2). We have also shown how to derive the multiplicity function from transit surveys (Section 3) assuming a given form for the inclination distribution (the Fisher distribution, Section 3.1). Our principal conclusions are as follows.

• 1.  
  At present, the Kepler data alone (Borucki et al. 2011) do not constrain the inclination distribution of multi-planet systems. Models with all rms inclinations—from razor-thin to spherical—are able to reproduce the observable multiplicity function in the Kepler sample. This conclusion is somewhat different from the study of Lissauer et al. (2011b), who found that their simulations tended to underpredict the number of systems with a single transiting planet (tranet for short), and in addition found evidence for a separate group of nearly coplanar systems. These conclusions may reflect the restricted, though plausible, range of models for the multiplicity function considered by Lissauer et al. (2011b). Their estimated upper limit to the mean inclination of 10° is consistent with our conclusions below.
• 2.  
  Systems with large rms inclinations are only consistent with the Kepler data if at least some of them contain a large number of planets. The relation between rms inclination and maximum number of planets is given by Equation (38).
• 3.  
  In our models, the percentage of one-tranet systems with additional planets is 20%–30%, and for two- or three-tranet systems this percentage is even higher (Figure 6). These fractions can be probed observationally using transit timing variations.
• 4.  
  The rms inclination can be constrained by combining estimates of the observable multiplicity function from Kepler and RV surveys, but only after estimating the effective number of stars that have been examined in RV surveys. We have made two estimates, one using Kepler data and one without; these are consistent, and yield 〈sin ²i〉^1/2 ⩽ 0.09, corresponding to mean inclinations in the range 0°–5°.
• 5.  
  Although the range of rms inclinations is tightly constrained by this analysis, the multiplicity function is not: the data are well fit by many pathological models containing no zero-planet systems, no one-planet systems, etc.

This research was supported in part by NASA grant NNX08AH83G and has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. Work by S.D. was performed under contract with the California Institute of Technology (Caltech) funded by NASA through the Sagan Fellowship Program. We acknowledge helpful conversations with Dan Fabrycky, Debra Fischer, Matt Holman, Boaz Katz, Darin Ragozzine, and Jason Wright.