Are Planetary Systems Filled to Capacity? A Study Based on Kepler Results
aa r X i v : . [ a s t r o - ph . E P ] F e b D RAFT VERSION M ARCH
1, 2013
Preprint typeset using L A TEX style emulateapj v. 11/10/09
ARE PLANETARY SYSTEMS FILLED TO CAPACITY? A STUDY BASED ON
KEPLER
RESULTS J ULIA F ANG AND J EAN -L UC M ARGOT
Draft version March 1, 2013
ABSTRACTWe used a sample of
Kepler candidate planets with orbital periods less than 200 days and radii between1.5 and 30 Earth radii ( R ⊕ ) to determine the typical dynamical spacing of neighboring planets. To derive theintrinsic (i.e., free of observational bias) dynamical spacing of neighboring planets, we generated populations ofplanetary systems following various dynamical spacing distributions, subjected them to synthetic observationsby the Kepler spacecraft, and compared the properties of observed planets in our simulations with actual
Kepler detections. We found that, on average, neighboring planets are spaced 21.7 mutual Hill radii apart with astandard deviation of 9.5. This dynamical spacing distribution is consistent with that of adjacent planets inthe Solar System. To test the packed planetary systems hypothesis, the idea that all planetary systems are dynamically packed or filled to capacity, we determined the fraction of systems that are dynamically packedby performing long-term (10 years) numerical simulations. In each simulation, we integrated a system withplanets spaced according to our best-fit dynamical spacing distribution but containing an additional planet onan intermediate orbit. The fraction of simulations exhibiting signs of instability provides an approximate lowerbound on the fraction of systems that are dynamically packed; we found that ≥ ≥ ≥
45% of2-planet, 3-planet, and 4-planet systems are dynamically packed, respectively. Such sizeable fractions suggestthat many planetary systems are indeed filled to capacity. This feature of planetary systems is another profoundconstraint that formation and evolution models must satisfy.
Subject headings: methods: statistical – planetary systems – planets and satellites: general – planets and satel-lites: detection INTRODUCTION
We examine the question of whether planetary systems gen-erally consist of closely spaced planets in packed configu-rations or whether planets in the same system are generallymore widely spaced apart. Here we adopt the traditional def-inition of dynamical spacing as the separation between adja-cent planets in terms of their mutual Hill radius, and we definea planetary system to be dynamically packed if the system is“filled to capacity”, i.e., it cannot accept an additional planetwithout leading to instability.The degree of packing in planetary systems has impor-tant implications for their origin and evolution. It hasbeen codified in the packed planetary systems (PPS) hy-pothesis (e.g., Barnes & Quinn 2004; Raymond & Barnes2005; Raymond et al. 2006; Barnes & Greenberg 2007),the idea that all planetary systems are dynamicallypacked. Previous works have invoked the PPS hypothe-sis to predict the existence of additional planets in sys-tems with observed planets located far apart with an in-termediate stability zone (e.g. Menou & Tabachnik 2003;Barnes & Quinn 2004; Raymond & Barnes 2005; Ji et al.2005; Rivera & Haghighipour 2007; Raymond et al. 2008;Fang & Margot 2012b), since the PPS hypothesis requiresthat an undetected planet is located in that stability zone. Sys-tems that are observed to have dense configurations could sup-port the PPS hypothesis if they were shown to be dynam-ically packed. Such systems may include Kepler-11, withsix transiting planets within 0.5 AU (Lissauer et al. 2011b),Kepler-36, whose 2 known planets have semi-major axesdiffering by only ∼
10% (Carter et al. 2012), and KOI-500, Department of Physics and Astronomy, University of California, LosAngeles, CA 90095, USA Department of Earth and Space Sciences, University of California, LosAngeles, CA 90095, USA which has 5 planets all within an orbital period of 10 days(Ragozzine et al. 2012).In this study we seek to investigate the underlying distribu-tion of dynamical spacing in planetary systems by fitting tothe observed properties of
Kepler planet candidates. By un-derlying or intrinsic , we mean our best estimate of the truedistribution of dynamical spacing between neighboring plan-ets in multi-planet systems, i.e., free of observational biases.After we derive the underlying distribution of the dynamicalspacing between planets, we create planetary systems whoseplanets have separations that obey this distribution. We thensubject these planetary systems to N-body integrations to ex-amine their stability properties, which allows us to determineif they are dynamically packed or not. By determining thefraction of systems that are packed, we can test the PPS hy-pothesis.In a related study published by Fang & Margot (2012a), weinvestigated the underlying multiplicity and inclination dis-tribution of planetary systems based on the Kepler catalogof planetary candidates from Batalha et al. (2012) in Febru-ary 2012. We created population models of planetary systemsfollowing different multiplicity and inclination distributions,simulated observations of these systems by
Kepler , and com-pared the properties of detected planets in our simulationswith the properties of actual
Kepler planet detections. Weused two types of observables: numbers of transiting systems(i.e., numbers of singly transiting systems, doubly transitingsystems, triply transiting systems, etc.) and normalized tran-sit duration ratios. Within our orbital period and planet radiusregime ( P ≤
200 days , . R ⊕ ≤ R ≤ R ⊕ ), we found thatmost planetary systems had 1 - Kepler mission.This paper is organized as follows. In Section 2.1, we de-fine our stellar and planetary parameter space. We also de-scribe how we created model populations of planetary sys-tems and how we compared them to the properties of
Kepler planetary candidates. In Section 2.2, we present the best-fitmodel representing our best estimate of the intrinsic distribu-tion of dynamical spacing in planetary systems. In Section 3,we compare this distribution of dynamical spacing with that ofthe Solar System. We also make comparisons with two othersystems, Kepler-11 and Kepler-36, to quantify how rare suchsystems are. In Section 4, we test and quantify whether such adistribution of dynamical spacing implies that planetary sys-tems are dynamically packed, by performing an ensemble ofN-body integrations. We briefly describe implications for thePPS hypothesis. Section 5 summarizes the main conclusionsof this study. DYNAMICAL SPACING OF PLANETS
Methods
Our methods for deriving the intrinsic dynamical spacingof planetary systems are as follows. First, we created modelpopulations of planetary systems obeying different underlyingdistributions of dynamical spacing. Second, we performedsynthetic observations of the planetary systems in these pop-ulations by the
Kepler spacecraft. At this stage we identi-fied which simulated planets were detectable by the
Kepler telescope, and which were not. Third, we compared the re-sulting distribution of dynamical spacing of detectable plan-ets from synthetic populations with that of the actual
Ke-pler detections. The actual distribution can be easily ob-tained from
Kepler transit data with an assumed planet radius-mass relationship. Most of these steps are fully explained inFang & Margot (2012a), and we refer the reader to that paperfor details. In the following paragraphs, we summarize themost salient points of our procedure as well as any differenceswith Fang & Margot (2012a).Each model population consists of about 10 planetary sys-tems, and we created various model populations that followeddifferent underlying distributions of multiplicity, inclinations,and dynamical spacing. To generate these populations, weneeded to restrict the range of physical and orbital proper-ties of the stars and planets that we considered in our sim-ulations. We selected ranges that would adequately over-lap those of a Kepler sample that can be considered nearlycomplete (Howard et al. 2012; Youdin 2011). Stellar proper-ties such as radius R ∗ , stellar noise σ ∗ , effective temperature T eff , surface gravity parameter log( g ), and Kepler magnitude K were randomly drawn from the Kepler Input Catalog (seeFang & Margot 2012a). We only considered bright solar-likestars that obeyed the following ranges:4100 K ≤ T eff ≤ , . ≤ log( g [cm s - ]) ≤ . , (1) K ≤
15 mag . Planet radii and orbital periods were drawn from debiased dis-tributions, and we obtained these debiased distributions byconverting the observed sample of Kepler Objects of Inter-est (KOI; based on detections up to Quarter 6 released inFebruary 2012, Batalha et al. 2012) into a debiased sample us-ing calculations of detection efficiencies (see Fang & Margot2012a). We filtered the KOI sample (and correspondingly limited the parameter space of the synthetic populations de-scribed in this paper) to the following orbital period P , planetradius R , and signal-to-noise ratio (SNR) boundaries: P ≤
200 days , . R ⊕ ≤ R ≤ R ⊕ , (2)SNR( → Q8) ≥ . . These limits were imposed in order to choose a sample ofplanets with properties unlikely to be missed by the
Kepler detection pipeline. For SNR, we required an SNR ≥
10 forobservations up to Quarter 6, which corresponds to aboutSNR ≥ √ N , where N is the number of ob-served transits. This scaling is performed because the SNRsof observed KOIs have been reported for observations up toQuarter 8 in Batalha et al. (2012), whereas the actual detec-tions have been reported up to Quarter 6 only. Planet masses M were calculated by converting from planet radii R . Weused a broken log-linear M ( R ) prescription obtained by fittingto masses and radii of transiting planets (see Fang & Margot2012a): log (cid:18) MM Jup (cid:19) = 2 . (cid:18) RR Jup (cid:19) - .
261 (3)for (cid:18) RR Jup (cid:19) < . , log (cid:18) MM Jup (cid:19) = - . (cid:18) RR Jup (cid:19) + .
777 (4)for (cid:18) RR Jup (cid:19) ≥ . . Additionally, we repeated all of the methods described inthis section by using an alternate mass-radius relationship:( M / M ⊕ ) = ( R / R ⊕ ) . (Lissauer et al. 2011b). By adopt-ing this alternate mass-radius equation, our results showedthe same best-fit dynamical spacing distribution as definedin Equation (7) with σ =14.5 (see results presented in Sec-tion 2.2). We note that both of these mass-radius equa-tions were obtained by fitting to the sample of planets withknown masses and radii. Errors in the mass-radius relation-ships can potentially affect our results because our determi-nation of dynamical spacing is a direct function of plane-tary masses. In Figure 1, we investigate how uncertaintiesin the mass-radius relationship map into dynamical spacinguncertainties. Specifically, we plot histograms showing howthe observed dynamical spacing changes if there is a 1 - σ in-crease or decrease in mass for the mass-radius equation. Formasses lower than nominal, adjacent planets appear to be lessclosely spaced and so the distribution shifts to the right. Formasses higher than nominal, adjacent planets appear to bemore closely spaced and so the distribution shifts to the left.The shifts are quantified at the end of Section 2.2.For orbital eccentricities, we adopted circular orbits, as wedid in Fang & Margot (2012a). Eccentricities do not directlyaffect our calculation of dynamical spacing, as we will definebelow in Equation (5).Regarding the multiplicity distribution in our model pop-ulations, we used a bounded uniform distribution with λ =1 . - . λ and is defined as follows: first, draw a value2 M a ss [ M J up ] N o r m a li z ed N u m be r Figure 1. (Top) Observed transiting exoplanets with knownmasses and radii shown in green, and the corresponding mass-radius relationship. (Bottom) Dependence of the observeddynamical spacing distribution on the choice of mass-radiusrelationship. In both panels, the nominal mass-radius rela-tionship as defined in Equations (3) and (4) is shown in black,and the mass-radius relationship shifted by one sigma to lower(higher) masses is shown in brown (red). ∆ represents thenumber of mutual Hill radii between adjacent planets in multi-planet systems and is defined in Equation (5). N max (maximum number of planets) from a Poisson distribu-tion with parameter λ that ignores zero values, and second,draw the number of planets from a discrete uniform distribu-tion with range 1 - N max (Fang & Margot 2012a). Thus, eachplanetary system will have at least one planet. For the inclina-tion distribution of planetary orbits, we used a Rayleigh dis-tribution with σ = 1 , ◦ as well as a Rayleigh of Rayleigh dis-tribution with σ σ = 1 , ◦ . A Rayleigh of Rayleigh distributionhas a single parameter σ σ and is defined as follows: first, drawa value σ from a Rayleigh distribution with parameter σ σ , andsecond, draw a value for inclination from a Rayleigh distribu-tion with parameter σ (Lissauer et al. 2011b). These specificmultiplicity and inclination distributions were chosen becausethey yielded fits consistent with transit numbers and transitduration ratios from Kepler detections (see Fang & Margot2012a). Combinations of these specific multiplicity and in-clination distributions add up to a total of 36 possibilities.The difference between model populations generated in Fang & Margot (2012a) and the model populations generatedin this study is the treatment of planetary separations, sincehere we wish to determine the underlying dynamical spacingof planetary systems. We used a separation criterion ∆ toassess the dynamical spacing between all adjacent planets inmulti-planet systems, where ∆ is defined as (e.g., Gladman1993; Chambers et al. 1996) ∆ = a - a R H , , (5)with R H , = (cid:18) M + M M ∗ (cid:19) / a + a . (6)In these equations, a is the semi-major axis, R H , is the mu-tual Hill radius, and M is the mass. Subscripts ∗ , 1, and 2refer to the star, the inner planet, and the outer planet, respec-tively. For a two-planet system not in resonance, the planetsare required to be spaced with ∆ & .
46 in order to be Hillstable.In our model populations, adjacent planets in multi-planetsystems were spaced according to a prescribed ∆ distribution.We used a shifted Rayleigh distribution, which is the same asa regular Rayleigh distribution except shifted to the right by3.5 (since we require this distribution to provide values of ∆ that meet the minimum Hill stability limit). Such a distribu-tion, as we will show, matches the observed sample well. Themathematical form of a shifted Rayleigh distribution f is f ( ∆ ) = ∆ - . σ e - ( ∆ - . / (2 σ ) , (7)and is described by a single parameter σ . In our model pop-ulations, we explored values of σ =10 -
20 with increments of0.5 for a total of 21 possibilities. We chose this range of σ values based on the location of the observed ∆ distribu-tion (blue histogram in Figure 2) with its approximate peakat about 20 mutual Hill radii. This chosen range of σ val-ues allowed us to explore different distributions of dynamicalspacing that spanned a reasonable range of possible model ∆ distributions. Increments of 0.5 were chosen as a trade-offbetween resolution and computational limitations. As will beseen in Section 2.2, the statistically good match between thedata and the best-fit model demonstrates that our incrementsare sufficiently small and have appropriately sampled the pos-sible range of ∆ distributions.In order to enforce that adjacent planets are spaced ac-cording to the prescribed ∆ distribution, we performed thefollowing steps. For each synthetic planetary system, thefirst planet’s orbital period is drawn from the debiased pe-riod distribution. If the system’s multiplicity is greater thanone, for the second planet we draw its separation ∆ from thefirst planet using the prescribed ∆ distribution and we alsodraw a value from the debiased period distribution. If thatvalue is less/greater than the first planet’s period, then the sec-ond planet will be the inner/outer planet and its exact periodwill be calculated by satisfying the drawn ∆ separation fromthe first planet. This process repeats if the system has addi-tional planets. These steps are different from Fang & Margot(2012a), where in that study all planetary orbital periods werechosen by drawing them from a debiased distribution. Weverified that the periods drawn to match the ∆ distributionprovide a very close match to the debiased period distributionas well.3fter the creation of each model population, we performedsynthetic observations of each population’s planetary systemsby the Kepler spacecraft in order to determine which planetswere transiting and detectable (see Fang & Margot 2012a).The transiting requirement was evaluated by picking a ran-dom line-of-sight (i.e., picking a random point on the celestialsphere) and computing the planet - star distance projected onthe plane of the sky. The minimum of that distance was com-pared to the radius of the host star to determine if the planetin our simulations transited or not. The detection requirementwas assessed by calculating each transiting planet’s SNR, de-fined as SNR = (cid:18) RR ∗ (cid:19) √ N σ ∗ , (8)where the first fraction gives the depth of the transit, N repre-sents the number of transits up to Quarter 6, and σ ∗ representsstellar noise (Combined Differential Photometric Precisionor CDPP; Christiansen et al. 2012). Since CDPP is quarter-to-quarter dependent, we used the median CDPP value overall available quarters. In the calculation of SNR, we tookinto account gaps between Kepler quarters, the fact that notall stars are observed each quarter, and a 95% duty cycle(Fang & Margot 2012a). If the calculated SNR for a transit-ing planet met or exceeded the SNR threshold for detection(SNR=10), then it was considered detectable.Lastly, we determined the goodness-of-fit between eachmodel population’s detected planets and the actual
Kepler de-tections. This was ascertained by comparing the ∆ distribu-tions of adjacent planets in their multi-planet systems. Weperformed a Kolmogorov-Smirnov (K-S) test to assess the fitbetween the ∆ distributions, and this comparison yielded ap-value that we used to evaluate the null hypothesis that thedistributions emanate from the same parent distribution. ThisK-S probability was used to determine how well a particu-lar model matched the observations. We also calculated thegoodness-of-fit for multiplicity (by comparing with observed Kepler numbers of transiting systems using a chi-square test)and for inclination (by comparing with observed, normalizedtransit duration ratios using a K-S test) to check that they wereconsistent with the data (see Fang & Margot 2012a). Whilewe only generated model populations with underlying multi-plicity and inclination distributions that are considered to begood fits to the data based on our previous work, this extrastep allowed us to confirm that any models with acceptable ∆ fits also produced acceptable multiplicity and inclination fitsto the Kepler data. We determined which model populationswere most consistent with the data by combining (multiply-ing) the probabilities associated with each one of the 3 statis-tical tests that probed multiplicity, inclination, and dynamicalspacing. We assumed that these probabilities are independent.Accounting for all combinations of multiplicity, inclination,and dynamical spacing distributions, in total we created 756model populations with about 10 planets each. As describedearlier, each of these model populations underwent syntheticobservations by Kepler as well as statistical tests. The nextsection presents our results.
Results
We report our results on the dynamical spacing (representedby the criterion ∆ ) in multi-planet systems based on Kepler data. Using the methods described in the previous section, wedetermine that our best-fit model for the intrinsic ∆ distribu- tion is a shifted Rayleigh distribution (see Equation (7)) with σ = 14 . ∆ = 21 . ∆ separations larger than 20, and about 90% ofneighboring planet pairs have ∆ separations larger than 10.This best-fit distribution was obtained by considering shiftedRayleigh distributions with increments in σ of 0.5. Themean values for distributions with σ = 14 . σ = 15 . ∆ = 22 . ∆ = 21 .
0, respectively. The combined probabil-ities for the match to the data using distributions with σ = 14 . σ = 15 . σ = 14 . R ⊕ and a maximum orbital pe-riod of 200 days. It is possible that planets are actually evenmore closely spaced than this best-fit distribution if there areplanets located in intermediate locations with radii less than1.5 R ⊕ . Therefore, our findings about the ∆ distribution canbe used to represent the spacing of planetary systems confinedto the scope of our study, or can serve as an upper limit for thespacing of planetary systems that include planets with smallerradii.Figure 3 shows the comparison between the ∆ distributionof simulated detections from this best-fit model’s populationand the observed ∆ from actual Kepler detections. Note thatthis figure does not show the underlying ∆ distribution thatis plotted in Figure 2; instead, this figure shows the distribu-tion of simulated planets that would have been detected , inorder to make an appropriate comparison with the observed ∆ distribution. The K-S test for comparing these two dis-tributions yields a p-value of 56%, indicating that we cannotreject the null hypothesis that these distributions are drawnfrom the same parent distribution. In other words, this modelis consistent with the observations.Comparison between Figures 2 and 3 shows that the under-lying ∆ distribution is similar to the observed ∆ distribution–both distributions have peaks near ∆ ∼
20. This suggests thatthe observed ∆ distribution is not indicative of a significantpopulation of non-transiting and/or low-SNR planets (withinthe planet radius and orbital period limits of our study) lo-cated in-between detected planets; otherwise, the underlying ∆ distribution would have on average lower ∆ values thanthose of the observed distribution. We caution again that ourstudy cannot rule out the existence of a population of planetswith R < . R ⊕ , so the actual underlying ∆ distribution couldbe different from what our results indicate. The similarity be-tween the observed and underlying ∆ distributions is due tothe rarity of cases where a system has ∆ observed > ∆ underlying (i.e., an undetected planet located in-between two detectedplanets). The stringent geometric probability of transit meansthat outer planets are more easily missed than inner planets.As a result, it will be rare to find cases where an intermedi-ate planet is non-transiting and therefore missed, but both theinnermost and outermost planets are transiting and detected.For such a case, which rarely occurs, the observed ∆ wouldbe greater than the underlying ∆ .We discuss how we expect these results to change if an al-ternate mass-radius relationship is used. In particular, Fig-ure 1 shows how the observed dynamical spacing distribution4 P r obab ili t y D en s i t y σ = σ =
20 40 60 80Underlying Δ00.20.40.60.81 C u m u l a t i v e P r obab ili t y σ = σ = Figure 2.
The best-fit model’s underlying ∆ distribution is shown as a magenta-colored curve. This distribution represents our bestestimate of the true or intrinsic (i.e., free of observational bias) distribution of dynamical spacing between all neighboring planetsmeeting our orbital period and planet radius criteria. Recall that ∆ represents the difference in semi-major axes between twoadjacent orbits; it is expressed in units of the mutual Hill radius. The best-fit model is a shifted Rayleigh distribution as definedin Equation (7) with σ = 14 .
5; the top plot depicts the probability density and the bottom plot shows the cumulative probability.The black-colored curves show the range and sampling frequency of model distributions that follow different σ parameter valuesranging from σ = 10 to σ = 20 with increments of 0.5, as examined in this study.changes depending on various choices for the mass-radius re-lationship. For computational expediency, we chose to evalu-ate the effects of the mass-radius relationship on the observed ∆ distribution. We use this as an approximation for the effectson the underlying ∆ distribution, with the justification that theobserved and underlying ∆ distributions appear similar (seeprevious paragraph). For the mass-radius relationship shifteddown by 1 sigma, the best-fitting shifted Rayleigh distributionhas σ = 16.5. For the mass-radius relationship shifted up by1 sigma, the best-fitting shifted Rayleigh distribution has σ = 12.5. As a result, we estimate that our derived dynamicalspacing distribution can span σ =12.5 - COMPARISON TO THE SOLAR SYSTEM, KEPLER-11, ANDKEPLER-36
We can compare our results (i.e., the intrinsic dynamicalspacing or ∆ distribution) with the dynamical spacing distri-bution of the Solar System, if we extrapolate beyond the ra-dius and period limits (Equation (2)) of our study. Figure 4shows our intrinsic ∆ distribution of planetary systems over-plotted with a histogram of the ∆ distribution of the SolarSystem planets. From this figure, it is interesting to note thatthe distributions appear to be relatively similar and that theSolar System planets may be similarly spaced as most exo-planets in general. A K-S test between our cumulative ∆ dis-tribution and the sample of ∆ values between adjacent planetsin the Solar System yields a p-value of 66.2%, indicating thatthe Solar System ∆ distribution is consistent with that of Fig-ure 2.The orbital evolution of the planets in the Solar System isknown to be chaotic and unstable (e.g., Sussman & Wisdom1988; Laskar 1989, 1990; Sussman & Wisdom 1992; Laskar5
20 40 60 80Δ00.020.040.060.080.10.120.140.16 N o r m a li z ed N u m be r ObservedDetected From Best-Fit Model
10 20 30 40 50 60Δ00.20.40.60.81 C u m u l a t i v e D i s t r i bu t i on σ = σ = Figure 3.
The top plot shows the comparison between the ob-served ∆ distribution from actual Kepler detections (shown inblue) and the ∆ distribution from simulated detections in ourbest-fit model population (shown in orange); the K-S proba-bility for this match is 56%. The bottom plot shows this com-parison in cumulative form and its histogram points are con-nected by lines for easy viewing, and also shows the compar-ison with the properties of detected planets from other modelpopulations ( σ =10 to σ =20 with increments of 0.5; shownin black). The 99.5% confidence region of acceptable fits in-cludes model populations with σ ranging from 12.5 to 17.0.1994; Michtchenko & Ferraz-Mello 2001; Lecar et al. 2001).The inner Solar System can be potentially unstable withinthe Sun’s remaining lifetime due to a secular reso-nance (Batygin & Laughlin 2008). Laskar (1994) andLaskar & Gastineau (2009) have shown that inner planetscould be ejected or collide. Numerical simulations of theplanets in the outer Solar System suggest that they arepacked (Barnes & Quinn 2004; Raymond & Barnes 2005;Barnes et al. 2008). All of these results suggest that the So-lar System is dynamically packed. If we consider the SolarSystem to be dynamically packed then it is possible that otherplanetary systems with similar planet multiplicities and ∆ dis-tributions are also dynamically packed. This prompted us toverify whether planetary systems in general are dynamicallypacked (see Section 4).Kepler-11 is a planetary system with six known transit-ing planets in a closely spaced configuration (Lissauer et al.2011a). All six transiting planets have orbits smaller than theorbit of Venus, and five of the six transiting planets have or-bits smaller than the orbit of Mercury. This appears to be
20 40 60 80Δ00.20.40.60.81 C u m u l a t i v e D i s t r i bu t i on Best-fit Model Solar System
Figure 4.
Comparison of the cumulative ∆ distribution be-tween the best-fit model’s ∆ distribution (i.e., shiftedRayleigh distribution with σ = 14 .
5; purple solid line) and theSolar System’s ∆ distribution (i.e., histogram based on its 8planets; dotted green line).a very compact system, and we calculate the ∆ separationsof the innermost five planets to be ∆ b - c =5.7, ∆ c - d =14.6, ∆ d - e =8.0, and ∆ e - f =11.2. We did not calculate the ∆ sep-aration of the f - g planet pair because the mass of planetg is not known and only has an upper limit. Accountingfor the 1 σ uncertainties on mass reported by Lissauer et al.(2011a), the dynamical spacing of these pairs have the follow-ing ranges: ∆ b - c =5.1 - ∆ c - d =13.0 - ∆ d - e =7.2 - ∆ e - f =10.0 - ∆ b - c =5.7is more closely spaced than 98.9% of adjacent planet pairsin multi-planet systems, ∆ c - d =14.6 is more closely spacedthan 74.6%, ∆ d - e =8.0 is more closely spaced than 95.3%,and ∆ e - f =11.2 is more closely spaced than 86.8%. Thesehigh percentages indicate that the planetary separations in theKepler-11 system are much smaller than average separationsin planetary systems, and we conclude that Kepler-11 is un-usual in terms of the density of its configuration.Kepler-36 has two known transiting planets with a largedensity contrast (their densities differ by a factor of ∼
8) yetthey orbit closely to one another (semi-major axes differ by ∼ ∆ =4.7. In comparison toour intrinsic ∆ distribution of dynamical spacing, a separationof ∆ =4.7 is more closely spaced than 99.7% of neighboringplanet pairs of planetary systems in general. DYNAMICAL PACKEDNESS OF PLANETS
Section 2.2 described the best-fit ∆ distribution of planetarysystems based on Kepler data, with a mean value of ∆ = 21 . ∆ implies that planetary systems are dynamicallypacked or not. By dynamically packed , we refer to a plan-etary system that is filled to capacity and cannot include anadditional planet without leading to instability.To investigate whether planetary systems are dynamically6acked, we performed long-term N-body integrations of plan-etary systems generated for our best-fit model population (seeSections 2.1-2.2). For each multiplicity (i.e., 2-planet sys-tems, 3-planet systems, 4-planet systems), we randomly chose1000 planetary systems for which we performed long-term in-tegrations. In total, we performed 3000 integrations. We didnot include single planet systems because they are irrelevantfor studies of dynamical packedness and we did not includesystems with multiplicities higher than 4 planets because theyare relatively rare for our parameter space (Fang & Margot2012a).For each integration, we added an additional planet whentesting for stability; this planet had a mass equal to the low-est mass of all original planets and its initial conditions in-cluded an orbital eccentricity of zero, an inclination drawnfrom a Rayleigh of Rayleigh distribution with σ σ = 1 ◦ (seeSection 2.1), and random values for its argument of pericen-ter, longitude of the ascending node, and mean anomaly. Thisadditional planet was placed in-between the orbits of exist-ing planets, and if there were 3 or 4 original planets, we ran-domly determined which 2 adjacent planets would be receiv-ing a new neighbor. The additional planet’s semi-major axiswas calculated so that it was located with equal mutual Hillradii distances between its neighboring planets. These initialconditions for the additional planet (e.g., low eccentricities,low inclinations, a mass equal to the lowest mass of origi-nal planets, a semi-major axis located at equal ∆ distancesfrom neighboring planets) are very conservative in the sensethat we have chosen initial conditions that are very amenableto stability, as we determine whether a planetary system withthis additional planet can remain stable or not.Our simulations were performed using a hybridsymplectic/Bulirsch-Stoer integrator from the Mercury package (Chambers 1999), and we used a timestep thatcovered 1/25 of the innermost planet’s orbital period. Simula-tions were performed for a length of 10 years; the instabilitytimescales had median values less than 10 years. Possibleoutcomes included either a stable system with no instabilitiesor a system with at least one instability defined as a collisionbetween the star and a planet, a collision between planets,and/or an ejected planet. All of these planetary systems wereverified to be stable for 10 years before adding the additionalplanet.The results of our simulations can be divided into twocamps. The first group is composed of planetary systems thatbecame unstable in our integrations. This suggests that theseplanetary systems are dynamically packed, since the additionof another planet in an intermediate orbit resulted in an unsta-ble planetary system. We found that 31%, 35%, and 45% of2-, 3-, and 4-planet systems were unstable, respectively (Ta-ble 1). The second group is composed of planetary systemsthat did not exhibit any signs of instability. For these systems,although they were stable within the scope of our integrations,they may still be dynamically packed. Possible reasons in-clude: an instability may occur on a longer timescale than ourintegration time, there may be additional planets in the sys-tem beyond the scope of our orbital period range of 200 days,or there may be additional planets in the system smaller than1.5 R ⊕ , which is the minimum radius of our study. All ofthese factors would affect the determination of the stability ofthe system, and so for this second group of systems we areagnostic about their dynamical packedness.Accordingly, we can only confidently provide a lower limit Table 1
Lower Limits on the Percentage of Dynamically Packed Systems
System Multiplicity Percentage of Packed Systems ≥ ≥ ≥ on packed systems by concluding that at least 31 -
45% (de-pending on the system’s multiplicity) of planetary systemswith dynamical spacings consistent with our best-fit ∆ distri-bution (Figure 2) are dynamically packed. These lower limitsare also presented in Table 1. Note that systems with lowermultiplicity are more common (Fang & Margot 2012a).The packed planetary systems (PPS) hypothesis is the ideathat all planetary systems are dynamically packed, and there-fore cannot hold additional planets without becoming unsta-ble. The results of our long-term numerical integrations areconsistent with the PPS hypothesis, as we find a sizeablelower limit of 31 -
45% (depending on the system’s multiplic-ity) of planetary systems to be dynamically packed. CONCLUSIONS
We have generated model populations of planetary systemsand simulated observations of them by the
Kepler spacecraft.By comparing the properties of detected planets in our simula-tions with the actual
Kepler planet detections, we have deter-mined the best-fit distribution of dynamical spacing betweenneighboring planets. This best-fit distribution is our best esti-mate of the underlying (i.e., free of observational bias) dis-tribution of dynamical spacing for our orbital period ( P ≤
200 days) and planet radius (1 . R ⊕ ≤ R ≤ R ⊕ ) parameterregime. Stemming from this distribution, the main results ofthis study are:1. On average, neighboring planets are spaced 21.7 mutualHill radii apart, with a standard deviation of 9.5. This distancerepresents the typical dynamical spacing of neighboring plan-ets for all systems included in the parameter space describedabove.2. Our best-fit distribution of dynamical spacing is consis-tent with the dynamical spacing of neighboring planets in theSolar System, with a K-S p-value of 66.2%. If we consider theSolar System to be dynamically packed, then it is not unrea-sonable to ask whether other planetary systems with similardynamical spacing are also dynamically packed.3. Based on our best-fit distribution of planetary spacing: ≥
31% of 2-planet systems, ≥
35% of 3-planet systems, and ≥
45% of 4-planet systems are dynamically packed. Thismeans that such systems are filled to capacity and cannot holdanother planet in an intermediate orbit without becoming un-stable.4. Our results on the dynamical packedness of planetarysystems are consistent with the packed planetary systems hy-pothesis that all planetary systems are filled to capacity, as wefind sizeable lower limits on the fraction of systems that aredynamically packed.5. Compact systems such as Kepler-11 and Kepler-36 rep-resent extremes in the dynamical spacing distribution. For ex-7mple, the two known planets in Kepler-36 are more closelyspaced than 99.7% of all neighboring planets represented byour orbital period and planet radius regime.We thank the referee for valuable comments. This researchwas performed using resources provided by the Open Sci-ence Grid (OSG), which is supported by the National ScienceFoundation and the U.S. Department of Energy’s Office ofScience.