The statistics of multi-planet systems
aa r X i v : . [ a s t r o - ph . E P ] J un The statistics of multi-planet systems
Scott Tremaine and Subo Dong Institute for Advanced Study, Princeton, NJ 08540, USA
ABSTRACT
We describe statistical methods for measuring the exoplanet multiplicityfunction—the fraction of host stars containing a given number of planets—fromtransit and radial-velocity surveys. The analysis is based on the approximation ofseparability—that the distribution of planetary parameters in an n -planet systemis the product of identical 1-planet distributions. We review the evidence thatseparability is a valid approximation for exoplanets. We show how to relate theobservable multiplicity function in surveys with similar host-star populations butdifferent sensitivities. We also show how to correct for geometrical selection ef-fects to derive the multiplicity function from transit surveys if the distribution ofrelative inclinations is known. Applying these tools to the Kepler transit surveyand to radial-velocity surveys, we find that (i) the Kepler data alone do not con-strain the mean inclination of multi-planet systems; even spherical distributionsare allowed by the data but only if a small fraction of host stars contain largeplanet populations ( &
1. Introduction
The distribution of inclinations in multi-planet systems provides fundamental insights intoplanet formation. The small inclinations of the planets in the solar system—the largestis 7 ◦ , for Mercury—strongly suggest that they formed from a disk. However, we shouldnot be surprised if extrasolar planetary systems have larger inclinations, for several rea-sons: (i) the rms inclinations in the asteroid and Kuiper belts are substantially larger, 12 ◦ and 16 ◦ respectively; (ii) in most astrophysical disks, the rms eccentricity and inclinationare correlated, and the eccentricities of extrasolar planets are much larger than those of Sagan Fellow • The mutual inclination of planets B and C in the system surrounding the pulsarB1257+12 is less than ∼ ◦ (Konacki & Wolszczan 2003); this result is only marginallyrelevant to planetary systems around main-sequence stars since pulsar planets musthave had a very different history. • Using radial-velocity and astrometric data, Bean & Seifahrt (2009) estimate that themutual inclination between GJ 876 b and c is 5 . ◦ +3 . ◦ − . ◦ . Using radial velocities anddynamical modeling of the planet-planet interactions Correia et al. (2010) concludethat the mutual inclination is . ◦ , while Baluev (2011) finds that the same quantityis between 5 and 15 ◦ . The large scatter among these results means that they shouldbe used with caution. • McArthur et al. (2010) find from astrometric and radial-velocity measurements thatthe mutual inclination of υ And c and d is 30 ◦ ± ◦ , much larger than in GJ 876 butstill small enough to suggest formation from a disk. • Dynamical fits to the transit timing of two planets in the Kepler-9 system yield anupper limit to the mutual inclination of ∼ ◦ (Holman et al. 2010). However, thissystem was discovered in a transit survey, and such surveys are far more likely to detectmulti-planet systems with small inclinations rather than large ones. • Lissauer et al. (2011a) studied the six-planet system Kepler-11 and concluded that theabsence of transit duration changes in Kepler-11e implies that its inclination relative 3 –to the mean orbital plane of other planets is less than 2 degrees ; once again, this resultis biased by the strong dependence of the probability that two or more planets willtransit on their mutual inclination.As noted above, if one planet in a two-planet system transits its host star as viewedfrom Earth, the probability that the second planet will also transit is higher if the mutualinclination of the two planetary orbits is small (e.g., Ragozzine & Holman 2010). This ar-gument suggests that the numbers of 1-planet, 2-planet, . . . , N -planet systems detected in alarge transit survey contain information about both the multiplicity function—the fraction ofhost stars containing 0 , , , . . . , N planets—and the inclination distribution. The challengeis to disentangle the two distributions to distinguish thick systems with many planets fromthin systems with few planets.The first attempt to do this was made by Lissauer et al. (2011b), who modeled thenumber of multiple-planet systems detected in the first four months of data from the Keplersurvey (Borucki et al. 2011)—115 with two transiting planets, 45 with three, 8 with four, andone each with 5 and 6. Lissauer et al. used a variety of simple models for the distributionof the number of planets per system. They found that none of their models fit the datawell, mostly because they produced too few systems in which a single transiting planet wasobserved, but that the best-fit models typically had mutual inclinations . ◦ .The purpose of this paper is to develop a general formalism that relates the intrinsicproperties of multi-planet systems to the properties of the multi-planet systems that aredetected in transit or other surveys ( § § §
4) and to radial-velocity surveys ( § First we introduce some notation. (i) The Kepler team uses the term planet “candidate” todenote a possible planet that has been discovered through transits but not yet been confirmed Lissauer et al. also concluded that the mean mutual inclination of the planets was 1–2 ◦ from MonteCarlo simulations of the probability that a randomly placed observer would see transits of all the planets;however, this conclusion is suspect since the probability that a random star with six planets would showsix transits is different from the probability that one star from the Kepler sample of ∼ ,
000 stars wouldshow six transits. survey selection effect is a limitation on the number of detectableplanets due to the detection thresholds. A geometrical selection effect is a limitation arisingfrom the orientation of the planetary system—in particular, the planet must cross in frontof the stellar disk to be detectable in a transit survey .We assume that the stars in a survey may have 0 , , . . . , K planets and denote thenumber of stars in the survey with k planets by N k . Thus P Kk =0 N k is the total number ofstars in the survey. The vector N = ( N , N , . . . , N K ) is called the multiplicity function.Because of survey and geometric selection effects, only a fraction of these planets willbe detected in the survey. Let the survey selection matrix element S km be the probabilitythat 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 planetspass the geometric selection criteria. Then the expected number of systems that the surveyshould detect with j tranets is n j = K X k = j G jk K X m = k S km N m , or n = G · S · N . (1)We call n the observable multiplicity function. Clearly G mn = S mn = 0 for m > n, G = S = 1 , G mn , S mn ≥ . (2)Moreover since the number of detectable planets in an n -planet system must be between 0and n , we have n X m =0 G mn = n X m =0 S mn = 1 . (3) There is also a geometrical selection effect in radial-velocity surveys, since the reflex velocity is propor-tional to sin γ where γ is the inclination of the planetary orbit to the line of sight. However, we can eliminatethis effect by working only with the minimum mass M sin γ where M is the planet mass; of course, for transitsurveys sin γ ≃ G and S are ( K + 1) × ( K + 1) upper-triangular stochastic matrices. For physicalreasons G and S should commute (eq. 1 should not depend on whether we consider thesurvey selection effects or the geometric selection effects first). We have confirmed that thecommutator [ G , S ] is indeed zero for the selection matrices that we derive below. Let w represent all of the orientation-independent properties of a planet and its host star thatdetermine its detectability (planetary mass and radius; stellar mass, radius, distance, andluminosity; orbital period, etc.) and let f ( w , . . . , w n ) represent the probability distributionof these parameters for an n -planet system. Thus R d w · · · d w n f ( w , . . . , w n ) = 1.A natural assumption for describing multi-planet systems is that the n -planet distribu-tion function is separable , that is, f ( w , . . . , w n ) = n Y m =1 f ( w m ) , Z d w f ( w ) = 1 . (4)This assumption can only be approximately valid—for example, it is inconsistent with theobservational 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 areunstable. Nevertheless, we argue that the separability assumption is sufficiently accurate toprovide a powerful tool for analyzing the statistics of multi-planet systems. We describe theevidence on its validity in §
2. Survey selection effects
Let Θ A ( w ) be the probability that a planet with properties w is detected in the surveylabeled by A if its host star is on the target list for this survey and the orientation of theobserver is correct (we assume that whether or not a planet can be detected is independentof the presence or absence of other planets in the same system, which is a reasonable firstapproximation). Thus the function Θ A ( w ) describes the survey selection effects for A, but notthe geometric selection effects. The probability that a planet is detected, ignoring geometricselection effects, is then W A = Z f ( w )Θ A ( w ) d w . (5) 6 –If the survey target list contains N Am stars with m planets, then using the separability as-sumption (4) the expected number of systems in which k planets will be detected is n Ak = K X m = k S km ( W A ) N Am , ≤ k ≤ K ; (6)where the survey selection matrix S is a ( K + 1) × ( K + 1) matrix whose entries are givenby the binomial distribution, S km ( W ) ≡ m ! k !( m − k )! W k (1 − W ) m − k , ≤ k ≤ m ≤ K, (7)and zero otherwise. Note that S (1) is the unit matrix. A useful identity is (Strum 1972,e.g.,) S ( A ) · S ( B ) = S ( AB ) , (8)which in turn implies S − ( W ) = S ( W − ) . (9)Although the physical motivation (6) for the definition of S ( W ) requires 0 ≤ W ≤
1, thematrix is well-defined for all values of W .With the assumption of separability it is straightforward to show that the conditionalprobability distribution of the parameters w m , given that k planets are detected, is (cf. eq.4) f ( w , . . . , w k ) = k Y m =1 f ( w m ) . (10)Thus a separable distribution is still separable after survey selection effects are applied, solong as the selection effects depend only on the properties of an individual planet.The factor W (eq. 5) is usually difficult to determine reliably since (i) we do not havegood models for the distribution f ( w ) of the planetary parameters; (ii) in most cases thesurvey selection effects Θ( w ) are not known accurately; (iii) in many cases the target listfrom which a given sample of exoplanets was detected is not even known (the Kepler survey isan exception to the last two limitations). However, useful results can be obtained without anexplicit evaluation of W . Suppose, for example, we have two surveys A and B that examinepopulations 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 Bm = cN Am where c is a constant given by the ratio of the number of target stars in B andA. Equation (6) can then be written n A = S ( W A ) N A , n B = c S ( W B ) N A . (11) 7 –Applying equations (8) and (9), we have n B = c S ( f BA ) n A (12)where f BA ≡ W B /W A = 1 /f AB . Thus the observable multiplicity function n B of survey Bis 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 targetstars in the two surveys). The parameter f BA can be eliminated if we plot n B , n B , . . . asfunctions of n B . In practice we must use the multiplicity function n A rather than n A onthe right side of equation (12) but these should not be very different so long as n Ak ≫ f BA is larger or smaller than unity, but if f BA > n Bk 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 observablemultiplicity 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 § n Bk , k = 1 , . . . , are functionsonly of f BA and the known n A , which approximates the observable multiplicity function n A .Hence by eliminating f BA in favor of n B , the number of k -planet systems in any survey Bcan be predicted as a function of the number of one-planet systems in that survey. Thesepredictions for k = 2 , , σ confidence bands (dashed lines). The actual numbers of multi-planet systemsafter SNR cuts on the Kepler data are shown as open circles. Within the statistical errorsthe predictions agree with the data for k = 2 and 3 and are marginally consistent for k = 4:using a Kolmogorov–Smirnov (KS) test , the p -value (probability of observing deviations atleast as extreme as those seen, given the null hypothesis) is 0.28, 0.27, and 0.06 respectively.The upper right panel of Figure 1 shows a similar comparison for a sequence of catalogsbased on cuts at increasing planet radius, rather than SNR. The results are consistent with The use of a KS test is not strictly applicable since n k and n are cumulative distributions of a thirdparameter, the SNR, rather than being directly related. However, the results should be approximately correctwhen n ≫ n k which is usually the case. K RV ). Open circles show n , n , n (numbers of 2, 3, and 4-planet systems)as a function of n . Solid and dashed curves show the predictions of equation (12) and the1– σ errors on the predictions. 9 –the separable model to within the statistical errors ( p -values of 0.72, 0.19, and 0.66 for k = 2 , , K RV (semi-amplitude of the radial-velocity curve) on theleft and orbital period on the right. The predictions are marginally consistent with the nullhypothesis ( p -values between 0.03 and 0.10) except for n as a function of the cut in K RV ,for which the null hypothesis is excluded.These results confirm that in many cases the separability approximation and equation(12) provide useful tools for removing survey selection effects and converting the observablemultiplicity function between surveys.
3. Geometric selection effects in transiting systems
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 only0.1–0.25, so this assumption should not cause significant errors. We shall also assume thata 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 muchsmaller than the stellar radius (the median ratio of planetary radius to stellar radius in theKepler 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 willdetect m tranets in this system. One planet
First consider the case n = 1. We define three unit vectors: ˆ o points towardsthe observer, ˆ n is normal to the planetary orbit, and ˆ z is normal to the reference plane fromwhich inclinations i are measured. Thus ˆ z · ˆ n = cos i and ˆ o · ˆ n = cos γ . If the planet’s sizeis negligible, it transits if and only if | ˆ o · ˆ n | < ǫ or | cos γ | < ǫ so g ( ǫ ) = 1 − g ( ǫ ) = R | cos γ | <ǫ sin γ dγ R sin γ dγ = ǫ. (13) 10 – Two planets
Let h ( w ) = 1 if | w | < h ( ǫ − cos γ ) is unity and we may write h ( ǫ − cos γ ) = ∞ X ℓ =0 b ℓ ( ǫ ) P ℓ (cos γ ) (14)where P ℓ is a Legendre polynomial. From the properties of these functions we have b ℓ ( ǫ ) = ǫ, ℓ = 0 ,P ℓ +1 ( ǫ ) − P ℓ − ( ǫ ) , ℓ even, ℓ > , ℓ odd. (15)Now let ( θ, φ ) be the polar coordinates for ˆ o relative to the polar axis ˆ z , and (Ω − π, i ) thepolar coordinates for ˆ n . Then h ( ǫ − cos γ ) = 4 π ∞ X ℓ =0 b ℓ ( ǫ )2 ℓ + 1 ℓ X m = − ℓ Y ∗ ℓm ( θ, φ ) Y ℓm ( i, Ω − π ) . (16)Let the probability distribution of planetary inclinations be q ( i | κ ) di , where κ is a set of freeparameters 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 u ( x | ǫ, κ ) = Z di d Ω2 π q ( i | κ ) h ( ǫ − cos γ ) =4 π ∞ X ℓ =0 b ℓ ( ǫ )2 ℓ + 1 Z di q ( i | κ ) Y ∗ ℓ ( θ, Y ℓ ( i, ∞ X ℓ =0 Q ℓ ( κ ) b ℓ ( ǫ ) P ℓ ( x ) . (17)where Q ℓ ( κ ) ≡ Z π di q ( i | κ ) P ℓ (cos i ) , Q = 1 . (18)If a system contains two planets, the probability that both transit for a random orien-tation of the observer is g ( ǫ , ǫ , κ ) = R − dx u ( x | ǫ , κ ) u ( x | ǫ , κ )= P ∞ ℓ,n =0 b ℓ ( ǫ ) b n ( ǫ ) Q ℓ ( κ ) Q n ( κ ) R − dx P ℓ ( x ) P n ( x )= ∞ X ℓ =0 Q ℓ ( κ )2 ℓ + 1 b ℓ ( ǫ ) b ℓ ( ǫ ) . (19) 11 –Moreover the probability that one and only one of the two planets transits is g ( ǫ , ǫ , κ ) = R − dx { u ( x | ǫ , κ )[1 − u ( x | ǫ , κ )] + [1 − u ( x | ǫ , κ )] u ( x | ǫ , κ )] } = g ( ǫ , κ ) + g ( ǫ , κ ) − g ( ǫ , ǫ , κ ) (20)and the probability that no planets transit is g ( ǫ , ǫ , κ ) =1 − g ( ǫ , ǫ , κ ) − g ( ǫ , ǫ , κ )=1 − g ( ǫ , κ ) − g ( ǫ , κ ) + g ( ǫ , ǫ , κ ) . (21)For example, if the planets are distributed isotropically then q ( i ) di = sin i di , Q ℓ = δ ℓ and g ( ǫ , ǫ ) = ǫ ǫ . If the planets have zero inclination, it can be shown that g ( ǫ , ǫ ) = ∞ X ℓ =0 b ℓ ( ǫ ) b ℓ ( ǫ )2 ℓ + 1 = min ( ǫ , ǫ ) , (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 : g mn ( ǫ , . . . , ǫ n , κ ) = R − dx P P n Q mi =1 u ( x | ǫ p i , κ ) Q nj = m +1 [1 − u ( x | ǫ p j , κ )] , (23)where P n is the set of all permutations ( p , . . . , p n ) of the numbers 1 , . . . , n , and m ≤ n . Forexample, g ( ǫ , ǫ , ǫ , κ ) = g ( ǫ , ǫ , κ ) + g ( ǫ , ǫ , κ ) + g ( ǫ , ǫ , κ ) − g ( ǫ , ǫ , ǫ , κ ) g ( ǫ , ǫ , ǫ , κ ) = g ( ǫ , κ ) + g ( ǫ , κ ) + g ( ǫ , κ ) − g ( ǫ , ǫ , κ ) − g ( ǫ , ǫ , κ ) − g ( ǫ , ǫ , κ ) + 3 g ( ǫ , ǫ , ǫ , κ ) . (24)The geometric selection matrix G mn (eq. 1) is simply h g mn ( R ⋆ /a , R ⋆ /a , . . . , R ⋆ /a l , κ ) i ,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 ( κ ) we use the separabilityassumption (4) with respect to ǫ = R ⋆ /a . Thus G mn ( κ ) = Z g mn ( ǫ , . . . , ǫ n , κ ) n Y k =1 f ( ǫ k ) d log ǫ k , (25) For n = 3 the functions g mn can be expressed as series in the Wigner 3- j symbols, but in practice it issimpler to evaluate the integral (23) numerically for any n >
12 –where f ( ǫ ) d log ǫ represents the probability distribution of ǫ as modified by the survey selec-tion effects.With this parametrization and equations (17) and (23) it is straightforward to showthat G mn ( κ ) is given by the binomial distribution, G mn ( κ ) = n !2 m !( n − m )! Z − dx U m ( x | κ )[1 − U ( x | κ )] n − m = R − dx S mn [ U ( x | κ )] (26)where S mn is given by equation (7), U ( x | κ ) ≡ R f ( ǫ ) u ( x | ǫ, κ ) d log ǫ R f ( ǫ ) d log ǫ = ∞ X ℓ =0 Q ℓ ( κ ) B ℓ P ℓ ( x ) , (27)and B ℓ ≡ Z f ( ǫ ) b ℓ ( ǫ ) d log ǫ with Z f ( ǫ ) d log ǫ = 1 . (28)Since B ℓ does not depend on the unknown parameters κ of the inclination distribution it canbe evaluated once and for all at the start of any optimization procedure. It is straightforwardto show that the relations (2) are satisfied by these formulae, and that the matrices G and S commute. In numerical work we typically truncate infinite series such as (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 sepa-rability assumption. This result does not generally hold for geometric selection effects. Toillustrate this, consider the simple case of a population of stars containing two planets, withzero relative inclination. Write the probability distribution of ǫ = R ⋆ /a of two-planet sys-tems as f ( ǫ ) f ( ǫ ) d log ǫ d log ǫ (after survey selection effects but before geometric selectioneffects). Then using equation (22) it is evident that the probability distribution of two-tranetsystems is dp ( ǫ , ǫ ) = f ( ǫ ) f ( ǫ )min ( ǫ , ǫ ) d log ǫ d log ǫ , (29)which is not separable. Only for isotropic distributions do geometric selection effects preserveseparability. In this paper we model the probability distribution of the inclinations dp = q ( i | κ ) di as aFisher distribution, q ( i | κ ) = κ κ exp( κ cos i ) sin i. (30) 13 –The parameter κ is related to the mean-square value of sin i through h sin i i = Z di sin i q ( i | κ ) = 2 coth κκ − κ . (31)When κ ≪ κ → q ( i | κ ) = sin i , while for κ ≫ κ →∞ q ( i | κ ) = (2 i/s ) × exp( − i /s ) where s = (2 /κ ) / is the rms inclination and π / s = 0 . s is the meaninclination. The Rayleigh distribution is commonly used to model the inclination distri-bution of asteroids, Kuiper-belt objects, stars in the Galactic disk (where it is known asthe Schwarzschild distribution), etc. As κ → −∞ the Fisher distribution approaches aretrograde Rayleigh distribution.For the Fisher distribution, equation (18) becomes Q ℓ ( κ ) = r πκ I ℓ +1 / ( κ )sinh κ (32)where I denotes a modified Bessel function. There is limited evidence on the accuracy of the separability approximation for multi-planetsystems. First consider RV surveys, in which there are no geometric selection effects. Themost important survey selection effects depend only on the properties of an individual planetso an RV survey of a separable parent distribution should lead to a separable detecteddistribution (eq. 10).Wright et al. (2009) compare 28 multi-planet systems and a much larger number ofsingle-planet systems detected by RV surveys. They find that (i) the eccentricities in multi-planet systems are smaller (mean eccentricity 0.22, compared to 0.30 in single-planet sys-tems); (ii) the logarithmic semi-major axis distribution in multi-planet systems is flatter,without the pileup of hot Jupiters between 0 . au and 0 . au and the enhancement out-side 1 au that are seen in single-planet systems; (iii) multi-planet systems exhibit an over-abundance of planets with minimum mass between 0.01 and 0.2 Jupiter masses. Thesedifferences are incompatible with separability and statistically significant ( p < . § § p -value 0.20; see also Figure 2),which is consistent with separability. Presumably the pileup of hot Jupiters at small semi-major axes seen in the RV surveys is less prominent in the Kepler sample because the typicalplanetary mass is much smaller, and the jump outside 1 au is not seen because Kepler is notsensitive to these orbital periods.Latham et al. (2011) have shown that Kepler systems with multiple tranets are less likelyto include a giant planet (larger than Neptune) than systems with a single tranet. We confirmusing a KS test that the distributions of radii in the single- and multiple-tranet systems aredifferent (maximum difference in the cumulative probability distribution of 0.20). However,the results at the end of §
4. Estimating the inclination distribution and the multiplicity function fromthe Kepler survey4.1. Properties of the survey
The Kepler survey has a complex set of survey selection effects, which we do not attemptto model. The constraints on the multiplicity function that we derive therefore apply tothe population of planets in radius, semi-major axis, etc. that Kepler detects, whateverthat population may be (for a discussion of selection effects and completeness in the Keplercatalog see Howard et al. 2011 and Youdin 2011). If we denote the multiplicity function ofthis population by N and the observable multiplicity function of the Kepler survey by n
15 –then equation (1) becomes n = G · N . (33)The validity of this equation requires only the plausible assumption that the probability thatKepler will detect a given transiting planet around a given star is independent of whether itdetects other transits around the same star.To produce a more homogeneous sample, we trim the catalog of Borucki et al. (2011) toinclude only stars with effective temperatures between 4000 and 6500 K and surface gravitylog g > . , , , . . . tranets are n = 1 . × , n = 737 , n = 104 , n = 37 ,n = 7 , n = 1 , n = 1 , n k = 0 for k > . (34)We need to determine the function f ( ǫ ), where f ( ǫ ) d log ǫ is the fraction of planets inthe range d log ǫ given the intrinsic distribution of planets and the survey selection effectsfor Kepler. As usual ǫ = R ⋆ /a is the ratio of stellar radius to planetary semi-major axis; thestellar radius is determined from the host-star mass and surface gravity and the semi-majoraxis is determined from the host-star mass and the planetary orbital period. Figure 2 showsdata points for f ( ǫ ) from single-tranet systems (red points) and from planets in multi-tranetsystems (blue points). The data points have been 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 distributions for single-tranetand multi-tranet systems are quite similar, and can be adequately fit by the parametrization f ( ǫ ) = 0 .
656 ( ǫ/ǫ ) . ǫ/ǫ ) . , ǫ = 0 . , for ǫ > .
004 (35)and zero for ǫ < . ǫ & . . . au (Borucki et al. 2011), while the cutoff at ǫ . .
004 is due to thelimited timespan of the Kepler data. 16 –Fig. 2.— The probability distribution of ǫ = R ⋆ /a , the ratio of stellar radius to planetarysemi-major axis, for tranets detected by Kepler. The differential probability distribution is dp = f ( ǫ ) d log ǫ . The data points for single- and multi-tranet systems are shown separately.The solid line shows the analytic fitting formula (35). 17 – The probability that the survey actually detects { n , n , . . . , n K } stars having 0 , , . . . , K planets is P ( n | N , κ ) = K Y k =0 n n k k exp( − n k ) n k ! (36)where n = ( n , n , . . . , n K ) and n k is related to N by equation (33).Estimating the multiplicity function N and the inclination distribution parameters κ from n is a straightforward but challenging problem in statistics and optimization. Thisproblem 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 most of these. In this paper we have usually chosen maximumlikelihood, as a reasonable compromise between generality, computation time, and clarity ofinterpretation.The log of the likelihood of a given observational result n islog P ( n | N , κ ) = K X k =0 n k log " K X l = k G kl ( κ ) N l − K X k =0 K X l = k G kl ( κ ) N l − K X k =0 log n k ! . (37)Note that the second term on the right can be simplified to P l N l using equation (3). We thenmaximize log P with respect to N and κ , subject to the constraint N k ≥ k = 0 , . . . , K . The top panel of Figure 3 shows the maximum likelihood as a function of the rms inclinationand the maximum number of planets per system, K , for 6 ≤ K ≤
40. The minimum allowedvalue is K = 6 since Kepler has found one system with six tranets. The maximum-likelihoodmodels 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 theoptimization algorithm (a quasi-Newton algorithm from NAG) converged on a local ratherthan global maximum. The figure shows that:(i) 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: systemsexist at all rms inclinations—even isotropic systems—with log likelihood only 0.73 smallerthan the razor-thin solutions. 18 –Fig. 3.— (top) The maximum likelihood of solutions for the multiplicity function of theKepler survey, as a function of rms inclination and maximum number K of planets persystem. Solid lines connect solutions with a given K , 6 ≤ K ≤
40; lines for K = 10 , , , σ (∆ ln L = 4 . h sin i i = 0. (bottom) Plots of χ (eq. 39) forthe maximum-likelihood models shown above. We estimate that models with χ . σ level (log likelihood smaller than the maximumby 4.5, marked by a horizontal dashed line on the figure), the maximum rms inclination isrelated to the maximum number of planets by h sin i i / ≤ (cid:26) .
15 + 0 . K − , K < ) / (isotropic) , K ≥ . (38)It is possible, of course, that even the maximum-likelihood model does not fit the datawell. To explore this possibility, we have calculated the standard Pearson χ statistic, χ = K X k =0 ( n k − n k ) n k = K X k =0 ( n k − P l G kl N l ) P l G kl N l . (39)The distribution of the χ statistic is not straightforward to interpret, since n k . k and since the number of degrees of freedom is not well-defined. Nevertheless it is probablyreasonable to expect that there is a good fit to the data if χ .
5. The values of χ for themaximum-likelihood solutions in the top panel of Figure 3 are shown in the bottom panelof that figure. There are satisfactory models with all rms inclinations, but as before suchmodels 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 (thebehavior of the isotropic solutions with K >
30 is qualitatively similar). The fraction ofstars with k -planet systems is N k P Kl =0 N l = . k = 0 , . k = 1 , k = 2 , . k = 3 , k = 4 , . . . , , . k = 30 . (40)Thus, in this solution, about half of the planets are contained in three-planet systems,and the other half in a small population ( < . , , , , . . . -planet sys-tems, as a function of the assumed rms inclination. The results are for K = 30 but arequalitatively similar for larger values of K . Our initial attempts to construct this figurewere unsuccessful, because the appearance of the figure is very sensitive to cases when the 20 –Fig. 4.— The fraction of stars in the Kepler sample containing k -planet systems, as a functionof the rms value of sin i . The curves are labeled by k for k ≤
13 and curves with 7 ≤ k ≤ n k must lie within the 90% confidence interval determined through equation (36). The costfunction minimized the total number of planets but the result is insensitive to this choice. 21 –optimization algorithm settles on a local maximum of the likelihood. To avoid this diffi-culty, we re-cast the optimization as a problem in linear programming: we demanded thateach 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 P Kk =1 N k . This specifies a unique solution, if one exists.At the smallest inclinations ( h sin i i / < .
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 stronglydominated by 1-planet and n i -planet systems where n i varies monotonically with the rmsinclination—for example, n i = 12 when h sin i i / ≃ .
3. We caution that these resultsshould not be regarded as a prediction of the Kepler multiplicity function for a given rmsinclination.The need for many-planet systems is straightforward to understand. Consider theextreme case of an isotropic distribution. Then κ = 0 and q ( i | κ = 0) = sin i ; thus Q ℓ ( κ = 0) = δ ℓ from equation (18) and the orthogonality properties of the Legendre poly-nomials. Thus U ( x | κ = 0) = B (eq. 27) and using equation (26) G mn ( κ ) = n ! m !( n − m )! B m (1 − B ) n − m . (41)If all systems contain n planets, the ratio of the number of m -tranet systems to the numberof ( m + 1)-tranet systems is G mn G m +1 ,n = m + 1 n − m − B B , n ≥ m + 1 . (42)Using equations (28) and (35) we find that B = 0 . n /n = 7 . ± .
7. For comparison the ratio G n /G n is less than 7 . . . n ≥
9; thus any population dominated by systems with less than 9 planets willoverproduce 1-tranet systems relative to 2-tranet systems. Similarly, for the Kepler survey n /n = 2 . ± .
5, and G n /G n > . . . n ≥ K ),since it is given simply by the ratio of the total number of planets to the number of targetstars, divided by the probability that a single randomly oriented planet will transit (Youdin2011). Mathematically, h number of planets per star i = P Kk =1 kn k B P Kk =0 n k = 0 . . (43)The large open circles in Figure 5 show the probability that a system with one, two,or three tranets has additional planets. Typically the fraction of one-tranet systems with 22 –additional planets is 0.2–0.5, without a strong dependence on rms inclination. For two orthree 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).
5. Combining Kepler and radial-velocity surveys
As described in the Introduction, a comparison of the observable multiplicity functions ofplanetary systems detected by radial velocities and by transits can offer a powerful probeof the inclination distribution. The principal obstacle to making this comparison is thatthe masses and orbital periods of the planets detected through these two observationaltechniques are quite different, as illustrated in Figure 6, and the multiplicity functions inthese two regions of parameter space are likely to be different. In this section we use theseparability approximation and the methods of § n Kep k and n RV k systems containing k planets. Weassume that both surveys have similar target star populations (we cull the list of target starsin both cases to include only FGK dwarfs), with multiplicity function N for Kepler and c N for the RV survey, where c < S ( W Kep ) and S ( W RV )be the survey selection functions. We assume that there are no geometric selection effectsfor the RV surveys (cf. footnote 2). The generalization of equation (36) for the likelihood is P ( n Kep , n RV | N , κ ) = K Y k =0 ( n Kep k ) n Kep k exp( − n Kep k ) n Kep k ! K Y k =1 ( n RV k ) n RV k exp( − n RV k ) n RV k ! (44)where n Kep = G ( κ ) S ( W Kep ) N , n RV = c S ( W RV ) . (45)Notice that the second product in equation (44) starts at k = 1 since it is difficult todetermine accurately how many stars have been unsuccessfully examined for planets byRV 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 §
2, the likelihood actually depends only onthe ratio W RV /W Kep ).We determine the observable multiplicity function for RV planets using all planets withFGK dwarf host stars in the exoplanets.org database (Wright et al. 2011) as of August 2010, n RV1 = 162 , n
RV2 = 24 , n
RV3 = 7 , n
RV4 = 1 , n
RV5 = 1 , n RV k = 0 for k > , (46)for a total of 240 planets. The observable multiplicity function for Kepler planets is givenin equation (34). Figure 7 shows the maximum likelihood as a function of the rms incli- 23 –Fig. 5.— The horizontal blue line, composed of ∼ ≤ K ≤ M for transiting planets are computed from radii R using M =( R/R ⊕ ) . M ⊕ (Lissauer et al. 2011b, for a more accurate relation see eq. 47) and masses forradial-velocity planets are minimum masses M sin γ . 25 –nation and the maximum number of planets per system, K (top), as well as χ for thesemodels (bottom). The plots are noisier than Figure 3, presumably because the optimizationalgorithm was less successful at finding the global maximum likelihood, but otherwise looksimilar. In particular, systems with large rms inclinations are consistent with the data if andonly if they contain a large number of planets. Evidently adding data from RV surveys hasnot significantly tightened the constraints on the inclination distribution.We now show that adding information on the total number of target stars in the RVsurveys does allow the inclination distribution to be determined. Figure 8 shows the expectednumbers n Kep k and n RV k of k -tranet systems from the Kepler survey and k -planet systems fromthe RV surveys, as determined from the maximum-likelihood solutions described above. Eachpoint 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 witherror 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 theerror bars of the corresponding observed value; this is no more than a confirmation that ouroptimization code is performing properly. The blue points show the total number of stars inthe RV survey, n RVtot = P Kk =0 n RV k , as determined by the optimization code. The plot showsthat n RVtot is tightly correlated with the rms inclination, so an accurate characterization of thetotal number of RV target stars would enable the determination of the rms inclination.This task is challenging given the heterogeneous surveys that have produced the RVplanets known at the present time. We have used two distinct approaches, which we nowdescribe.(i) Cumming et al. (2008) carry out a careful examination of selection effects in the KeckPlanet Search, and derive the percentage of F, G, and K stars with a planet in variousranges of orbital period and mass. The sample of RV planets used in our analysis (eq. 46) isnot corrected for selection effects, but for sufficiently massive planets and sufficiently shortorbital periods it should be complete. For example, for planets more massive than Jupiter, M sin γ > M J , with orbital periods less than one year, P < K RV >
30 m s − , large enough to be detectable in most surveys. In this mass and periodrange our sample contains 46 planet-hosting stars and Cumming et al. (2008) estimate thatthe fraction of stars with planets is 0 . ± . n RVtot = 2400 ± P <
100 d gives n RVtot = 2500 ± M sin γ > . M J gives n RVtot = 1900 ±
500 (based on 63 host stars). Thislast estimate of n RVtot is probably low because the surveys we have used are not all completeat this level.(ii) We may estimate n RVtot using the tranet frequency derived from the Kepler mission. 26 –Fig. 7.— As in Figure 3, except the data include both the Kepler transit survey and radial-velocity surveys. 27 –Fig. 8.— The expected numbers of 0,1,2,3 tranet systems from the Kepler survey and of1,2,3 planet systems from RV surveys, as predicted by our models. The observed numbersare shown as error bars surrounded by circles. Also shown is the total number of targets inthe RV surveys as predicted by our models (blue points). 28 –
50 100 150 200200040006000 P (days)10 m/s15 m/s20 m/s25 m/s
Fig. 9.— The estimated number of host stars in RV surveys, as determined by comparisonwith the Kepler survey. The curves and associated error bars show the number of RV hoststars as estimated by comparing the number of RV and Kepler planets with period lessthan P and mass exceeding that required to induce a given velocity semi-amplitude K RV atperiod P . The observed number of Kepler planets is multiplied by g − (eq. 13) to correct forgeometric selection effects, and the conversion between radius and mass is given by equation(47). Results are shown for four semi-amplitudes, K RV = 25 , , ,
10 m s − ; the plot at thesmallest semi-amplitude is low because the RV surveys are incomplete at this level. 29 –Once again, we restrict the Kepler sample to host stars that are F, G, and K dwarfs(4000 K < T eff < g > P and velocity semi-amplitude K RV : (i) compute the corresponding mass M ( P, K RV ) = M J ( K RV /
30 m s − )(1 yr /P ) / assuming a circular orbit and a solar-mass hoststar; (ii) find the number n RV ( P, K RV ) of RV planets with period less than P and massgreater than M ( P, K RV ); (iii) find all Kepler tranets with mass greater than M ( P, K RV )and period less than P , using an empirical mass-radius relation found by fitting mass andradius measurements from transiting planets in the range 0.1–10 M J (see Figure 10) to alog-quadratic relationlog R/R J = 0 .
087 + 0 .
141 log
M/M J − .
171 (log
M/M J ) ; (47)(iv) compute the total number of Kepler planets in this range n Kep ( P, K RV ) by countingeach tranet as ǫ − planets, to correct for geometric selection effects (eq. 13); (v) estimate thetotal number of RV host stars as n RVtot = n Keptot n RV ( P, K ) /n Kep ( P, K ). The results are shownin Figure 9 for K RV = 10 , , ,
25 m s − . As the majority of RV surveys have reachedprecisions of ∼
15 m s − or better over the last decade, it is reassuring but not surprisingthat the estimates of n RVtot for K RV = 15 , ,
25 m s − are consistent. The rise in n RVtot at smallperiods is likely due to the known discrepancy in hot Jupiter frequency between transit andRV surveys (the frequency of hot Jupiters estimated from transit surveys is factor of ∼ n RVtot ≃ ± n RVtot ≃ ± < h sin i i / < .
08 and 0 . < h sin i i / < .
09 which correspond to anrms or mean inclination range of 0–5 ◦ (as shown in § § < P <
125 d, a radius cut 1 . R ⊕ ≤ R ≤ R ⊕ , and a signal/noise cut SNR ≥
16, which reduced the number of planets to 63%of our sample. We find the mean inclination for this sample to be 0–4 ◦ , not significantlydifferent from the estimate in the preceding paragraph.Although the range of rms inclinations is tightly constrained by this analysis, themultiplicity function is not. For example, within 1– σ of the maximum-likelihood model(∆ log P ≤ .
5) we have found models that have no 1-planet systems (67% have no planets, 30 –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 (eq. 45), which measures therelative sensitivity of the RV and Kepler surveys. This ratio varies smoothly from 0.5 forrazor-thin systems to 0.2 for h sin i i / = 0 .
1, independent of the maximum number ofplanets in the model. In other words 20–50% of the Kepler planets could have been detectedin 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 doesnot require estimating the total number of RV target stars.A weak link in these arguments is the assumption that the population of FGK dwarf starsis the same in the Kepler and RV surveys. One sign that these populations are different is thehigher frequency of hot Jupiters found in RV surveys, as mentioned above. However, we notethat our two approaches to estimating n RVtot , one using only RV surveys and one comparingthe Kepler and RV surveys, yield similar answers, which suggests that the estimate of therms inclination that we derive from this answer is insensitive to differences between the hoststars of the Kepler and RV surveys.It is interesting to compare this estimate of the mean inclination to the mean eccentricityfor Kepler planets. Restricting our sample to planets with minimum mass between 0.01 and0.1 Jupiter masses and period
P >
10 d (to avoid the effects of tidal circularization), themean eccentricity of planets discovered in RV surveys is 0.15 (we have also excluded planetswith 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 planetsfrom transit timing: Moorhead et al. (2011) find that the mean eccentricity is between 0.13and 0.25 at a p -value of 0.05. We have h i ih e i = 0 . h i i ◦ . h e i . (48)Theoretical studies of eccentricity and inclination growth in planetesimal disks (e.g., Ida et al.1993) find h i i / h e i = 0 . e ∼ .
04 in RV catalogs, and thebias in this sample is likely to be higher since the SNR is low for low-mass planets. Possiblya 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 favorablecases (Ford et al. 2011). These variations can be used to detect and characterize additionalplanets. Given the rms inclination of 0–0.09 radians that we have derived, roughly 20–30% 31 –
Fig. 10.— The masses and radii of confirmed transiting exoplanets. The green solid line isthe log-quadratic fit in equation (47). The red dashed line is the log-linear fit log(
M/M ⊕ ) =2 .
06 log(
R/R ⊕ ) from Lissauer et al. (2011b). 32 –of the single-tranet Kepler systems are expected to have additional planets (Figure 5), andmany of these may be detectable by transit timing variations. Ford et al. (2011) estimatethat ∼
6. Summary
We have described a methodology for analyzing the multiplicity function—the fraction ofhost stars containing a given number of planets—in radial-velocity (RV) and transit surveys.Our approach is based on the approximation of separability, that the probability distri-bution of planetary parameters in an n -planet system is the product of identical 1-planetdistributions ( § § §
3) assuming a given form forthe inclination distribution (the Fisher distribution, § ◦ were poor fits to the data. We believe that theseconclusions reflect the restricted, though plausible, range of models for the multiplicityfunction considered by Lissauer et al. (2011b), although their estimated upper limit tothe mean inclination is entirely consistent with our conclusions below based on othermethods.2. Systems with large rms inclinations are only consistent with the Kepler data if at leastsome of them contain a large number of planets. The relation between rms inclinationand maximum number of planets is given by equation (38). 33 –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 5). Thesefractions can be probed observationally using transit timing variations.4. The rms inclination can be constrained by combining estimates of the observable mul-tiplicity function from Kepler and RV surveys, but only after estimating the effectivenumber 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 h sin i i / ≤ . ◦ .5. Although the range of rms inclinations is tightly constrained by this analysis, themultiplicity function is not: the data are well-fit by (presumably) pathological modelscontaining no zero-planet systems, no one-planet systems, etc.This research was supported in part by NASA grant NNX08AH83G, and has made useof the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. Workby SD was performed under contract with the California Institute of Technology (Caltech)funded by NASA through the Sagan Fellowship Program. We acknowledge helpful conver-sations with Dan Fabrycky, Debra Fischer, Matt Holman, Boaz Katz, Darin Ragozzine, andJason Wright. REFERENCES
Abt, H. A. 2010, PASP, 122, 1015Baluev, R. V. 2011, arXiv:1105.4696Bean, J. L., & Seifahrt, A. 2009, A&A, 496, 249Black, D. C. 1997, ApJ, 490, L171Borucki, W. J., et al. 2011, arXiv:1102.0541Correia, A. C. M., et al. 2010, A&A, 511, A21Cumming, A., Butler, R. P., Marcy, G. W., Vogt, S. S., Wright, J. T., & Fischer, D. A. 2008,PASP, 120, 531Ford, E. B., et al. 2011, arXiv:1102.0544Gould, A., Dorsher, S., Gaudi, B. S., & Udalski, A. 2006, Acta Astron., 56, 1 34 –Holman, M. J., et al. 2010, Science, 330, 51Howard, A. W., et al. 2011, arXiv:1103.2541Ida, S., Kokubo, E., & Makino, J. 1993, MNRAS, 263, 875Konacki, M., & Wolszczan, A. 2003, ApJ, 591, L147Latham, D. W., et al. 2011, ApJ, 732, L24Lissauer, J. J., et al. 2011a, Nature, 470, 53Lissauer, J. J., et al. 2011b, arXiv:1102.0543McArthur, B. E., Benedict, G. F., Barnes, R., Martioli, E., Korzennik, S., Nelan, E., &Butler, R. P. 2010, ApJ, 715, 1203Moorhead, A. V., et al. 2011, arXiv:1102.0547Morton, T. D., & Johnson, J. A. 2011, arXiv:1101.5630Papaloizou, J. C. B., & Terquem, C. 2001, MNRAS, 325, 221Ragozzine, D., & Holman, M. J. 2010, arXiv:1006.3727Ribas, I., & Miralda-Escud´e, J. 2007, A&A, 464, 779Strum, J. E. 1972, Two-Year College Math. J., 8, 260Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010, ApJ, 718, L145Wright, J. T., Upadhyay, S., Marcy, G. W., Fischer, D. A., Ford, E. B., & Johnson, J. A.2009, ApJ, 693, 1084Wright, J. T., et al. 2011, PASP, 123, 412Youdin, A. N. 2011, arXiv:1105.1782Zakamska, N. L., Pan, M., & Ford, E. B. 2011, MNRAS, 410, 1895