@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Physics and Astronomy, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Granados Contreras, Águeda Paula"@en ; dcterms:issued "2018-11-09T16:56:03Z"@en, "2018"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description "The discovery of exoplanets on short orbital periods (P < 100 days), including hot-Jupiters and Systems with Tightly-packed Inner Planets (STIPs), defies predictions from classic planet formation theory. Their existence requires either large-scale migration of planets through disks or rethinking fundamental steps in the planet formation process or some combination of both. It is further unclear whether the known STIPs harbor additional, undetected planets at even larger stellarcentric distances, which would have fundamental implications for how the systems formed. Through numerical simulations, we explore: (1) the in-situ formation of hot-Jupiters as an extreme outcome of early metastability of STIPs in the presence of gas; and (2) the dynamical effects of distant gas giants on STIPs using two case studies. In addition, we use synthetic systems to explore whether hot-Jupiters could form in-situ within dynamically unstable STIPs through the consolidation of a critical core M > 10 Earth-masses. We compare the dynamical outcomes of gas-free and gas-embedded planetary systems, in which consolidation of a critical core was only possible in the gas-free simulations. In contrast, STIPs are resistant to instability when gas is present, resulting in coplanar and nearly circular systems. The instability of the configurations after 10 Myr increases if the eccentricity is perturbed to e~0.01. In some cases, the planet-disk interaction produces co-orbiting planets that are stable even when the gas is removed. We explore the transit detectability of these configurations and find that the coorbital transit signature is difficult to identify in current transit detection pipelines due to the system dynamics. To explore STIP evolution in the presence of an outer giant planet, we vary the semi-major axis of the perturber between 1 and 5.2 au. We find that the presence of the outer perturber, in most locations, only alters the STIP precession frequencies but not its evolution or stability. In those locations where the perturber causes secular eccentricity resonances, the STIP becomes unstable. Secular inclination resonances can affect the observed multiplicity of transiting planets by driving the orbits of one or more planets to inclinations about 16 degrees."@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/67781?expand=metadata"@en ; skos:note "Orbital outcomes of STIPs andconsequences for hot-Jupiterformation and planet diversitybyA´gueda Paula Granados ContrerasB.Sc. in Physics, Universidad de Guanajuato, Mexico, 2010M.Sc. (Astronomy), Universidad Nacional Auto´noma de Me´xico, Mexico, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR IN PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Astronomy)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)November 2018c© A´gueda Paula Granados Contreras 2018The following individuals certify that they have read, and recommend to the Faculty of Grad-uate and Postdoctoral studies for acceptance, the dissertation entitled:Orbital outcomes of STIPs and consequences for hot-Jupiter formation and planet diversitysubmitted by Agueda Paula Granados Contreras in partial fulfillment of the requirements forthe degree of Doctor in Philosophyin AstronomyExamining Committee:Dr. Aaron C. Boley, AstronomySupervisorDr. Brett Gladman, AstronomySupervisory Committee MemberDr. Ingrid Stairs, AstronomySupervisory Committee MemberDr. Douglas Scott, AstronomyUniversity ExaminerDr. Neil Balmforth, MathematicsUniversity ExaminerAdditional Supervisory Committee Members:Dr. Jaymie Matthews, AstronomySupervisory Committee MemberSupervisory Committee MemberiiAbstractThe discovery of exoplanets on short orbital periods (P . 100 days), including hot-Jupitersand Systems with Tightly-packed Inner Planets (STIPs), defies predictions from classic planetformation theory. Their existence requires either large-scale migration of planets through disksor rethinking fundamental steps in the planet formation process or some combination of both.It is further unclear whether the known STIPs harbor additional, undetected planets at evenlarger stellarcentric distances, which would have fundamental implications for how the systemsformed. Through numerical simulations, we explore: (1) the in-situ formation of hot-Jupiters asan extreme outcome of early metastability of STIPs in the presence of gas; and (2) the dynamicaleffects of distant gas giants on STIPs using two case studies. In addition, we use syntheticsystems to explore whether hot-Jupiters could form in-situ within dynamically unstable STIPsthrough the consolidation of a critical core M ≥ 10M⊕. We compare the dynamical outcomes ofgas-free and gas-embedded planetary systems, in which consolidation of a critical core was onlypossible in the gas-free simulations. In contrast, STIPs are resistant to instability when gas ispresent, resulting in coplanar and nearly circular systems. The instability of the configurationsafter 10 Myr increases if the eccentricity is perturbed to e ∼ 0.01. In some cases, the planet-disk interaction produces co-orbiting planets that are stable even when the gas is removed.We explore the transit detectability of these configurations and find that the coorbital transitsignature is difficult to identify in current transit detection pipelines due to the system dynamics.To explore STIP evolution in the presence of an outer giant planet, we vary the semi-majoraxis of the perturber between 1 and 5.2 au. We find that the presence of the outer perturber, inmost locations, only alters the STIP precession frequencies but not its evolution or stability. Inthose locations where the perturber causes secular eccentricity resonances, the STIP becomesunstable. Secular inclination resonances can affect the observed multiplicity of transiting planetsby driving the orbits of one or more planets to inclinations about 16◦.iiiLay summaryThe variety of exoplanets so far discovered demonstrates that astronomers do not understandplanet formation as well as previously thought. In particular, the formation mechanism ofplanetary systems that are compact and very close to their host star remains hotly debated.My main focus is on understanding these tightly-packed planetary systems, in which there area large multiplicity of planets, like in the Solar System, with orbits compressed within theequivalent of Earth’s orbital distance to the Sun. It is also uncertain whether these systemshave additional, undetected planets at much larger distances from their host star, makingthem like a more densely populated Solar System. With the aid of numerical simulations, Iexplore in this work an alternative planet-formation mechanism for short-period planets, aswell as investigating how these systems would interact with a distant companion and what theobservational signature might be.ivPrefaceThe research presented here is based on numerical simulations of planetary systems usingthe N-body integrators Mercury6 and the 15th-order Integrator with Adaptive Time-Stepping(IAS15), as described in Chapter 2. I modified the original Mercury6 code (Chambers, 1999,publicly available) to include additional physical interactions with the star, which are describedin Sections 2.1.1 and 4.1. IAS15 is my own implementation of the Rein and Spiegel (2015)algorithm (Section 2.1.4 and 3.1.2), and it is utilized in Section 3.3. Due to the number ofnumerical simulations (more than 500) needed for these studies, we used the Orcinus comput-ing cluster, provided by WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada(www.computecanada.ca).Throughout Chapter 4, the secular theory code resmap is employed to find the secular eigen-frequencies and eigenvectors for a given planetary system, as well as the eccentricity andinclination forcing that the system exerts on a test particle. The resmap code, briefly de-scribed in Section 2.2, was written by my supervisor Dr. Boley and is publicly available inhttps://github.com/norabolig/resmap.Some of the material discussed in Sections 3.1.1 and 3.2 has been published in an article co-authored by me and led by my supervisor (Boley, A. C.; Granados Contreras, A. P.; and Glad-man, B. The In Situ Formation of Giant Planets at Short Orbital Periods. Astrophys. J., Lett.817:L17, 2016). Dr. Boley coordinated the overall project, including writing the manuscript andanalyzing the bulk of the simulation output. Dr. Gladman was also involved in the conceptualdevelopment of in-situ formation of hot-Jupiters and in providing feedback for the manuscript.I set up and ran the 1000 N-body simulations for the study, helped to analyze the simulationoutput, and contributed to the writing of the manuscript, specifically the parts included in thisthesis.A large fraction of the material described in Chapter 4 has been previously published in arefereed journal article on which I was the lead author (Granados Contreras, A. P. and Boley,A. C. The Dynamics of Tightly-packed Planetary Systems in the Presence of an Outer Planet:Case Studies Using Kepler-11 and Kepler-90. Astron. J. 155:139, 2018). I was responsible forrunning the numerical simulations, analyzing the data, and writing the manuscript. Dr. Boleywas involved in the interpretation of the results, as well as the writing and editing of themanuscript.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Exoplanets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.1.1 Characteristics of detected planets. . . . . . . . . . . . . . . . . . . . . . 41.1.2 Hot and warm Jupiters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.1.3 Systems with Tightly-packed Inner Planets (STIPs) . . . . . . . . . . . . 131.1.4 Stability of extrasolar planets . . . . . . . . . . . . . . . . . . . . . . . . 141.1.5 Minimum Mass Solar Nebula (MMSN) vs. Minimum Mass ExtrasolarNebula (MMEN) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151.2 Planet formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161.2.1 Formation of rocky planets . . . . . . . . . . . . . . . . . . . . . . . . . . 161.2.2 Giant planets: core accretion plus gas capture model . . . . . . . . . . . 191.3 This thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Methodology overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.1 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.1.1 Beyond Newtonian Gravity: additional forces . . . . . . . . . . . . . . . 222.1.2 Close encounters and collisions . . . . . . . . . . . . . . . . . . . . . . . . 24vi2.1.3 Mercury6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.1.4 15th-order Integrator with Adaptive Time-Stepping (IAS15) . . . . . . . 292.2 Secular code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.2.1 Secular theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.2.2 Forced eccentricity and inclination . . . . . . . . . . . . . . . . . . . . . . 412.3 Synthetic secular theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423 Close-in giant planets as an extreme outcome of planet-planet scattering. . 453.1 Simulations and prototype STIPs. . . . . . . . . . . . . . . . . . . . . . . . . . . 473.1.1 Gas-free gravitational interactions . . . . . . . . . . . . . . . . . . . . . . 483.1.2 Planet-disk interaction: gaseous tidal damping . . . . . . . . . . . . . . . 493.2 Gas-free outcomes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.2.1 Before moving through clouds: comparing the outcomes of Mercury6 andIAS15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.3 Outcomes with tidal damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.3.1 Stability test after complete gas dissipation. . . . . . . . . . . . . . . . . 573.3.2 Very compact STIPs in the presence of gas. . . . . . . . . . . . . . . . . 583.4 Coorbital planets and their detection . . . . . . . . . . . . . . . . . . . . . . . . 614 STIPs and external planetary perturbers: effects on observability and sta-bility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.1 STIPs studied and simulations setup. . . . . . . . . . . . . . . . . . . . . . . . . 734.2 Unperturbed vs. perturbed STIPs: dynamical outcomes . . . . . . . . . . . . . . 764.3 Secular theory vs. synthetic secular theory . . . . . . . . . . . . . . . . . . . . . 804.3.1 Mean motion resonances MMRs . . . . . . . . . . . . . . . . . . . . . . . 925 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.1 Close-in giant planets as an extreme outcome of planet-planet scattering. . . . . 955.2 STIPs and external planetary perturbers: effects on observability and stability. . 985.3 STIPs and giant planets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1006 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1017 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105AppendicesA Complementary material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117A.1 Laplace-Lagrange literal expansion of the disturbing function. . . . . . . . . . . 117viiA.1.1 Direct term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117A.1.2 Indirect terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118A.2 Laplace-Lagrange solution for more than two secondary masses. . . . . . . . . . 118A.3 Orbital evolution of coorbital realizations in the corotating frame . . . . . . . . 121A.4 IAS15 integration procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122viiiList of Tables1.1 Multiplicity of extrasolar planetary systems. . . . . . . . . . . . . . . . . . . . . . 81.2 Summary of the growth sequence to form rocky planets starting from dust particles. 172.1 Numerical values of the Gauss-Radau spacings used in IAS15 . . . . . . . . . . . 313.1 Nominal orbital elements of Kepler-11 (K11) known planets. . . . . . . . . . . . . 473.2 Summary of randomized orbital elements of Kepler-11 (K11) analogues. . . . . . 493.3 Summary of orbital elements of synthetic STIPs. . . . . . . . . . . . . . . . . . . 493.4 Summary of results comparing outcomes of Mercury6 and IAS15. . . . . . . . . . 533.5 Summary of gaseous tidal evolution results. . . . . . . . . . . . . . . . . . . . . . 573.6 Eccentricity of planets in r-628. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.1 Nominal orbital elements of Kepler-90 (K90) known planets. . . . . . . . . . . . . 744.2 Main precession period and amplitude of each planet in K11, K11+, K90, andK90+. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91ixList of Figures1.1 Planetary sizes and masses as a function of orbital period of confirmed andcandidate extrasolar planets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2 Exoplanet period ratio distribution. . . . . . . . . . . . . . . . . . . . . . . . . . 91.3 Measured eccentricities versus orbital periods of detected exoplanets. . . . . . . . 92.1 Notation used in cylindrical coordinates. . . . . . . . . . . . . . . . . . . . . . . . 232.2 Values of the Gauss-Radau spacings. . . . . . . . . . . . . . . . . . . . . . . . . . 312.3 Comparison of relative energy after 100 orbits between Rein and Spiegel’s IAS15and our implementation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.4 Geometry used for the secular theory development. . . . . . . . . . . . . . . . . . 362.5 Graphical depiction of the eccentricity and inclination vectors. . . . . . . . . . . 393.1 In situ formation pathways of short-orbital period planets. . . . . . . . . . . . . . 463.2 Cumulative total mass distribution in synthetic systems . . . . . . . . . . . . . . 503.3 Fraction of critical cores per interval using K11 as a prototype STIP. . . . . . . . 523.4 Comparison of critical core formation through consolidation between Mercury6and IAS15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.5 Example of the dynamical evolution of a synthetic STIP embedded in a gaseousdisk. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.6 Distribution of period ratios for systems with f = 5.0. . . . . . . . . . . . . . . . 553.7 Distribution of period ratios for systems with f = 9.5. . . . . . . . . . . . . . . . 563.8 Final mutual Hill spacing distribution, f , for systems with initial separationsf = 5.0 and f = 9.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.9 Final distribution of period ratios for systems with initial f = 1.5 and f = 3.0. . 593.10 Final mutual Hill spacing distribution, f , for systems with initial separationsf = 1.5 and f = 3.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.11 Semi-major axes, a, pericenter, q, and apocenter, Q, as a function of time fortwo realizations with planet orbit crossings. . . . . . . . . . . . . . . . . . . . . . 623.12 Offset angle, φ, between the planets near the 1:1 MMR and a detailed radialevolution as a function of time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 633.13 Semi-major axes of coorbital planets as a function of time. . . . . . . . . . . . . . 65x3.14 Normalized deviation from a Keplerian orbit, s/a = (al − as)/a, as a function ofoffset angle, φ, for r-628. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.15 Transit lightcurves of a coorbital pair in a horseshoe orbit at two different epochs,and corresponding period change as a function of time. . . . . . . . . . . . . . . . 683.16 Phase folded transit lightcurves from Figure 3.15. . . . . . . . . . . . . . . . . . . 704.1 Evolution of the orbital parameters of a stable Kepler-11 (K11) analogue andK11+ system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.2 Evolution of the orbital parameters of a stable Kepler-90 (K90) analogue andK90+ system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.3 K90+ system with a highly inclined perturber with an I = 50◦ and with a = 5.2au. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.4 Orbital evolution of K11+ in the presence of a migrating outer perturber. . . . . 784.5 Orbital evolution of K90+ in the presence of a migrating outer perturber. . . . . 794.6 Secular map of the forced inclination and eccentricity of Kepler-11 (K11) in thepresence of a Jupiter-like perturber (excluding K11b and K11c in the calculations). 814.7 Secular map of the forced inclination and eccentricity of Kepler-11 (K11) in thepresence of a Jupiter-like perturber (excluding only K11b in the calculations). . . 824.8 Secular map of the forced inclination and eccentricity of Kepler-90 (K90) in thepresence of a Jupiter-like perturber. . . . . . . . . . . . . . . . . . . . . . . . . . 834.9 Comparison of the theoretical and synthetic apsidal precession frequencies andphases of Kepler-11 (K11) and K11+. . . . . . . . . . . . . . . . . . . . . . . . . 854.10 Comparison of the theoretical and synthetic apsidal precession frequencies andphases of Kepler-90 (K90) and K90+. . . . . . . . . . . . . . . . . . . . . . . . . 864.11 Comparison of the theoretical and synthetic nodal precession frequencies andphases of Kepler-11 (K11) and K11+. . . . . . . . . . . . . . . . . . . . . . . . . 874.12 Comparison of the theoretical and synthetic nodal precession frequencies andphases of Kepler-90 (K90) and K90+. . . . . . . . . . . . . . . . . . . . . . . . . 884.13 Zoom-in of the nodal principal component of planet b in Kepler-90 (K90) andK90+. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 904.14 Librating resonant angles in K90. . . . . . . . . . . . . . . . . . . . . . . . . . . . 93A.1 Phase-on corotational evolution of the realizations in Fig. 3.12 in the frame ofthe larger planet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121A.2 IAS15 integration procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123xiList of AcronymsBS Bulirsch-Stoer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26DFT discrete Fourier transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43EPE Extrasolar Planets Encyclopaedia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4FFT fast Fourier transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43FT Fourier transform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .43GR general relativity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22gSPG gas short-period giant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45HJ hot-Jupiter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1IAS15 15th-order Integrator with Adaptive Time-Stepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21K11 Kepler-11. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21K90 Kepler-90. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21MMEN Minimum Mass Extrasolar Nebula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15MMR mean motion resonance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8MMSN Minimum Mass Solar Nebula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15MVS mixed-variable symplectic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26RV radial velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2SPG short-period giant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46STIP System with Tightly-packed Inner Planets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2TDV transit duration variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94TPS Kepler Transit Planet Search pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67TTV transit timing variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3WJ warm-Jupiter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11xiiList of symbolsPtr transit detection probabilityR∗ stellar radiusRp planetary radiusM∗ stellar massMp planetary massP orbital perioda semi-major axise eccentricityI, i orbital inclination (not to be confused with the subindex i)$ longitude of pericenterΩ longitude of the ascending nodeM mean anomalyλ mean longitude (if not otherwise stated)characteristic disk density decay timescaleK radial velocity (RV) amplitudeG gravitational constantT temperaturer stellarcentric distanceNp number of planets in a system or configurationRmH mutual Hill radiusµ planet-star mass fractionη = f semi-major axis difference between two planets in units of RmHσsolids disk surface density of solidsσgas disk surface density of gasFdisk total disk mass relative to the MMSNZrel disk metallicity relative to the MMSNFD drag forceCD drag coefficientρ densitys particle’s radiusv velocity relative to gasxiiib effective cross sectional radiusvesc escape velocityσv velocity dispersion between objectsFg gravitational focusingΣp surface density of a swarm of planetesimalsM˙ mass growth rateΩp = n orbital frequency or mean motion of a planetMiso isolation massσρ(r) surface density of solids at a stellarcentric distance raGR General Relativity corrected acceleration vectorc speed of lightaTD gaseous tidal damping acceleration vectorv velocity vectorρ radial component of the position vector in cylindrical coordinateskˆ unit vector of the vertical componentvρ magnitude of velocity vector’s radial component (in cylindrical coordinates)vz magnitude of velocity vector’s vertical component (in cylindrical coordinates)τe eccentricity damping timescaleτi inclination damping timescalevc circular Keplerian velocityτwave characteristic wave timecs isothermal sound speedR gas constantτlocal local τwave scaled at 1 auH disk scale heightτ1auo scaled value of τlocal at 1 auσMMEN(ρ) MMEN surface density at a radial distance ρamig migration acceleration vectorτmig migration timescale at 1 aur position vector in Cartesian coordinatesr˙ velocity vector in Cartesian coordinates(q, p) canonical variables, generalized coordinates and momentat time variableτ timestep variableH HamiltonianK(rij) factor to account for close encounters in Mercury6’s hybrid integratorE error term of an approximation (in this case the Gauss-Radau quadrature)Pn(x) Legendre polynomials of order nhn n-spacing of Gauss-Radau quadraturexivgk integration constants of the IAS15bk integration constants of the IAS15Uk central potential felt by mass mkRk disturbing function of mass mkµk product of gravitational constant and mass mkαij ratio of semi-major axis of particle/body i and jα¯ij conditioned semi-major axis ratioRD direct part of the disturbing functionRI indirect part of the disturbing function due to an internal perturberRE indirect part of the disturbing function due to an external perturberb(k)s (αij) Laplace coefficientsAij elements of the secular matrix A associated with e and $Bij elements of the secular matrix B associated with I and Ω(hi, ki) eccentricity vector of particle/body i(pi, qi) inclination vector of particle/body igl eigenvalues of the secular matrix A of mode lβl phase associated with glfl eigenvalues of the secular matrix B of mode lγl phase associated with flefree free eccentricity of a test particleIfree free inclination of a test particleFt[f(t)](ν) Fourier transform (FT)Fk[fkN−1f=0 ](n) discrete Fourier transform (DFT)ν linear frequencyAn amplitude of the FTφ phase of the FToffset angle between two planets (as seen from the central mass)si difference of semi-major axes of two coorbiting planets at time iH(φ) function of the angle between two coorbiting planetsW average arc length of a coorbiting planets180 difference in semi-major axes of two planets when separated by φ = 180◦φmin minimum angle between two coorbiting planetsRl radius of the larger planet in a coorbital configurationRs radius of the smaller planet in a coorbital configurationq,Q pericenter and apocenter values of a planet, respectivelyΩ˙, $˙ nodal and apsidal precession frequenciesPnode, Pap nodal and apsidal precession periodsϕ resonant angle between two planetsxvAcknowledgmentsI acknowledged my supervisory committee for their support and patience after a few changesin the delivery of this thesis. I am grateful to my supervisor, Dr. Aaron Boley, for his support,feedback, time, comments, and patience. Also to Dr. Brett Gladman for his comments andsuggestions on my research and simulations.I am grateful to my parents for their support and for raising me as a strong and independentindividual. I could not have gotten this far without them.This research has made use of the NASA Exoplanet Archive, which is operated by the CaliforniaInstitute of Technology, under contract with the National Aeronautics and Space Administra-tion under the Exoplanet Exploration Program.xviChapter 1IntroductionThroughout history, astronomers and planetary scientists have been intrigued by the formationand evolution of planets. Before the discovery of 51 Peg b by Mayor and Queloz (1995),which is the first confirmed extra-solar planet (exoplanet) orbiting a main-sequence star, planetformation theory focused on explaining the origin of the Solar System (Lissauer, 1993). However,the characteristics of 51 Peg b required a general re-examination of planet formation theory.The planet is similar to Jupiter (Mp = 0.47MJ), but has an orbital period 1000 times shorter(P = 4.27 days). The surface temperature of 51 Peg b is expected to be as high as 1300 K, dueto its short orbital distance of 0.05 au, or 8 times closer to its star than Mercury is to the Sun.The planet became the first of a new class of planets: hot-Jupiters (HJs). At the time, Mayorand Queloz (1995) suggested that 51 Peg b formed at ∼ 5 au and then migrated inward by twoorders of magnitude in semi-major axis to reconcile its location with standard ideas for giantplanet formation.All modern planet formation theories posit that planet building takes place in young circum-stellar disks of gas and dust (a protoplanetary disk) (Lissauer, 1993; Armitage, 2007; Raymond,2010; Perryman, 2011, and references therein). These disks arise as a natural consequence ofthe gravitational collapse of interstellar material, in which high-angular momentum gas anddust cannot fall directly onto the newly-forming star. Although the details are not understood,planetesimals form from small solid condensates in the disk, as well as any dust that survivedthe disk formation process. Depending on a number of different factors, including the efficiencyof planetesimal formation, the growth of planetary embryos, and the availability of gas, a rangeof planetary types is possible. In this scenario, the planet’s birth location with respect to thehost star is of importance (Pollack et al., 1996), because the protoplanetary disk temperaturegenerally decreases with increasing distance from the star. This causes chemical differentiationamong both gaseous and solid disk materials to be a function of stellar distance. Depending onthe details of the disk, condensation thresholds can lead to rapid changes in the availability ofsolids, particularly ices, in some portions of the disk. Thus, some disk regions could favor therapid growth of planetary cores. This framework works well for explaining the Solar System’sarchitecture, which is divided into an inner and outer system, delineated by planetary size andcomposition. Rocky small planets make up the inner Solar System, while only gas and ice giantsare found in the outer system. For context, all Solar System giant planets are between 5 and130 times more distant from the Sun than is Earth. The discovery of massive planets at shortorbital periods, such as 51 Peg b, was therefore a surprise, when considering the ideas built onSolar System formation theory.While short-period giants are fairly rare, short-period planets are not. Several different ex-oplanet surveys have been carried out, using multiple detection techniques, with the mostsuccessful being the transit and radial velocity (RV) methods. As of May 31, 2018, 3735 plan-ets have been discovered and confirmed (NASA Exoplanet Archive, 2018). The orbital perioddistribution of the known exoplanets shows that 80% have a P ≤ 100 days. According to obser-vations, most of the stellar hosts with planets detected have a single confirmed planet (78%),while 614/2825 of stellar host have two or more planets. From these multiplanet systems,around 20% are tightly packed, which means that three or more planets in the system haveorbital periods shorter than 100 days (Schneider et al., 2011). These Systems with Tightly-packed Inner Planets (STIPs) can have a range of planetary masses, and are inherently differentfrom the Solar System’s inner configuration. The fraction of extrasolar systems with planetson short orbital periods is non-negligible,1 and it is now clear that the formation theory basedon the Solar System is either incomplete or inadequate to address the diversity of exoplanets.Several questions have arisen from the discovery of exoplanets, such as: why is the SolarSystem so different from the observed extrasolar planetary systems? What is missing in theformation mechanisms developed to date? How do short-period giant planets get to their currentlocations? Are there more distant (undetected) planets in the systems already discovered? Arethere other formation mechanisms that we should explore?As a consequence of these open questions, the formation scenario of STIPs is hotly debated,with two basic mechanisms in the literature: in-situ assembly/formation and planet migration(see Raymond et al., 2014, for a review). In the planet migration scenario, giant planets format disk distances at which water can condense (i.e., beyond the water ice line) and then migrateinwards due to interactions with the gaseous disk or due to dynamical interactions followedby tidal circularization. Initially, planet migration seemed to be the only viable formationmechanism, particularly for gas giants, due to the inferred lack of materials for planet buildingat short orbital periods (Lin et al., 1996). However, as will be discussed in the following section,the frequency of planets in or near period commensurabilities, a fundamental prediction ofconvergent migration, is lower than expected in STIPs (e.g., Hands et al., 2014; Fabrycky et al.,2014). Recent work has also shown that disks could be massive enough to form planets at shortorbital periods, especially if small solids migrate inward first due to aerodynamic drag (Hansenand Murray, 2012; Boley et al., 2014; Tan et al., 2015).In the following sections, we will provide a summary of the cumulative characteristics of the1In the solar neighborhood the occurrence rate of HJs is a tenth that of STIPs. This point will be discussedfurther throughout Sects. 1.1.1 to 1.1.3.2detected exoplanets (Sect. 1.1), as well as a description of the classical formation mechanismof terrestrial and giant planets (Sect. 1.2). Emphasis will especially be placed on planet con-figurations with short orbital periods and their formation and evolution. We will test thein-situ formation of short-period giant planets through numerical simulations (Chap. 3) basedon stability considerations discussed in Sect. 1.1.4.1.1 ExoplanetsThe number of known and confirmed exoplanets has increased markedly since the discovery ofthe first set of exoplanets orbiting the millisecond pulsar PSR B1257+12 (Wolszczan and Frail,1992). To date, we know about 4000 exoplanets, with different orbital periods, sizes and masses.Of these planets, 2919 were detected using the transit method, 673 using radial velocities (RVs),60 using microlensing, 15 using transit timing variations (TTVs) and 68 with other techniques(NASA Exoplanet Archive, 2018). The most efficient detection method for exoplanets so farhas been the transit method, with the Kepler mission making the largest contribution.In the transit method, photometric monitoring of a star looks for periodic decreases in thestellar flux, which could be due to a planet passing in front of the star. The orbital periodand the planet-to-star size ratio can be measured directly from the observed lightcurves. If thestellar size and mass are known, then the candidate’s orbital semi-major axis and planetaryradius can be derived. The transit detection probability averaged over the orbit orientation isgiven byPtr = R∗ +Rpa (1− e2) ,≈ 0.241− e2(Pday)−2/3(1 +RpR ),where R and Rp are the solar and planetary radii, a is the planetary semi-major axis, and eand P are the orbital eccentricity and period. For a given stellar size, larger planets at shortorbital periods are more likely to be detected with this technique. In addition, the transitmethod is biased towards configurations with low mutual inclinations due to the near edge-ongeometry necessary to observe a transit.The RV approach, which is the second most successful detection method, relies on measuringvariations in the star’s motion along the observer’s line-of-sight (i.e., the radial motion withrespect to the observer). As the planet and the star orbit each other, the stellar motion leads toperiodic Doppler shifts in the measured locations of stellar atmospheric absorption features. Inthe case of a single planet, those variations can be characterized by a period and an amplitude,3where the RV amplitude is given byK =(2piGP)1/3 Mp sin I(M∗ +Mp)2/31√1− e2 ,∝Mp sin I P−1/3,in which Mp and M∗ are the planetary and stellar masses, G the gravitational constant and Iis the angle between the orbital plane’s normal and the line-of-sight. For a given stellar mass,the RV method is most sensitive to massive planets and to planets at short-orbital periods.With the observed amplitude and period, the product Mp sin I can be determined, assumingthe stellar mass is known and that M∗ Mp. The planetary mass itself can only be derived ifthe orbital inclination is measured separately. It would seem that for both detection methods, alarge eccentricity increases the likelihood of detecting a planet, especially for the transit method.Nonetheless, a large eccentricity reduces the transit duration, depending on the observed phaseof the orbit, due to its dependence on 1 − e for a given planetary mass and orbital period. Ineither technique, there are additional variables that affect the planet detection, e.g., the age ofthe star and the presence of more than one planet in the system.There are different online databases that keep track of the orbital parameters of the newly foundexoplanets, e.g., NASA Exoplanet Archive (2018), Extrasolar Planets Encyclopaedia (EPE)2(Schneider et al., 2011) and the Exoplanet Orbit Database3, which are meant to facilitate thestatistical analysis of exoplanets. The total number of planets that each database providesdepends on their definition of a confirmed planet and whether they also include candidates.As an example, according to the NASA Exoplanet Archive (2018), there are 3735 confirmedplanets (updated to May 31, 2018), while the EPE cites 3781. From these nearly 4000 confirmedplanets, to date, there are 973 planets with measured masses, 449 with Mp sin I measured and2953 with measured radii (NASA Exoplanet Archive, 2018). Only a small fraction of planets,around 10% of the total confirmed planets, have both mass and size independently determined(Chen and Kipping, 2017; NASA Exoplanet Archive, 2018).1.1.1 Characteristics of detected planets.In this section, we describe some of the characteristics of the exoplanet population. The datapresented here come from a variety of sources, including published journal articles and exoplanetdatabases (Schneider et al., 2011; Wright et al., 2011; Han et al., 2014; NASA Exoplanet Archive,2018). Special attention is paid to the distribution of orbital periods, planetary sizes and masses.The planetary densities are only known for some exoplanets, as independent mass and radiusmeasurements are needed. For those planets that do have inferred densities, it is possible to2http://www.exoplanet.eu3http://www.exoplanets.org4constrain the planet’s bulk composition.The Kepler mission’s discoveries have demonstrated that planetary systems with multiple plan-ets on short orbital periods are common. Among the Kepler candidates, 23% of the host starsharbor multiple planet candidates; and 46% of all candidates reside in known multiplanetsystems (Burke et al., 2014). For the solar neighborhood, the frequency of planets at shortorbital periods is between 30 and 50% (Howard et al., 2012; Mayor et al., 2011), suggestingthat short-period multiplanet systems may be present around at least 5% of stars in the solarneighborhood. Likewise, distributions for planetary sizes and masses are emerging for theseplanets (Fig. 1.1), and demonstrate that planets with periods P < 100 days are diverse in size,mass, and composition (see Wright et al., 2011).The masses and sizes shown in figures throughout this introduction will be given in Jupiterunits. The conversion factors between Jovian and terrestrial size and mass are R⊕ ≈ 0.09RJand M⊕ ≈ 0.003MJ.Fig. 1.1 shows the orbital periods, planetary masses and radii of the available confirmed planetsand planet candidates in the EPE. This plot includes the combination of observables from thetransit and RV techniques, the strong biases inherent to these detection techniques have notbeen removed from the presented data. The top panel displays the orbital period distributionof: (1) all the exoplanets in the database with measured periods (black solid line); and (2)Jupiter analogues (planets with 0.1 < Mp < 60MJ and Rp > 0.5RJ). Different features areidentified in the period distributions. The majority of the detected planets have a period shorterthan 100 days, regardless of their size or mass. There is also a distinctive peak at 3 days forthe Jupiter analogue population, commonly known as the “3-day pileup” (e.g. Udry et al.,2003; Gaudi et al., 2005; Cumming et al., 2008), which is absent in the complete exoplanetpopulation (Howard et al., 2012; Fressin et al., 2013). This apparent pileup is meaningful forsome migration-driven formation mechanisms of HJs (Rasio and Ford, 1996; Weidenschillingand Marzari, 1996; Fabrycky and Tremaine, 2007; Beauge´ and Nesvorny´, 2012). The countsof both populations drop for 1 < P < 4 d, which seems to be a physical cutoff rather thanobservational bias (Cumming et al., 2008).The planetary size and mass distributions are shown on the right of Fig. 1.1. After almost30 years of exoplanet detection, it is clear that exoplanets come in a variety of sizes andmasses, ranging from a few lunar-radii to radii greater than 1 RJ,4 and from lunar-mass objects4Objects with sizes above 2 RJ are problematic to explain using solely the equation of state that describes gasgiant planets (Zapolsky and Salpeter, 1969). The corresponding mass-radius relationship has a critical mass ofaround 3MJ which correlates to a maximum radius of 1RJ. In general, intrinsic or extrinsic heating mechanismscould explain the slow cooling and contraction of these enlarged planets (for a description of these heatingmechanisms, see Ginzburg and Sari, 2016), but not sizes much larger than about 2RJ.5Figure 1.1: Planetary sizes and masses as a function of orbital period of confirmed and candidateextrasolar planets. The data from exoplanet.eu includes 3781 confirmed and 2723 candidateplanets Schneider et al. (2011). Top panel: orbital period distribution of all objects withmeasured periods (black line) and of Jupiter analogues (0.1 < Mp < 60MJ and Rp > 0.5RJ).Middle panels: planetary size as a function of orbital period (left) along with the planet sizedistribution (right). The colored circles correspond to planets with both mass and radiusmeasured while the open circles have only the planet’s size measured. The colorbar indicatesthe planet mass range in Jovian masses. Red shows Earth-mass planets, yellow Neptune-massplanets, green Saturn-mass planets, cyan Jovian mass planets, and blue brown dwarfs. Noticethat the objects with P > 2000 days have masses larger than 1MJ. On the histogram onthe right, the planet size distribution shows two clear peaks at 0.17RJ and 1.1RJ. Bottompanels: measured planetary masses as a function of the orbital periods. The red circles indicatemeasurements of Mp sin I. The histogram at the right is the observed mass distribution inJupiter masses. There are two peaks again the distribution, this time at 0.025MJ and 1.2MJ.6(Mp ∼ 10−5MJ) up to few tens of MJ5 in mass. Both distributions present two prominentpeaks. The primary peak in the size distribution is located at ∼ 0.17RJ ≈ 1.8R⊕, while asecondary peak is at 1.1RJ ≈ 12R⊕. A recent analysis of the size distribution of the Keplerplanets, as part of the California-Kepler Survey, found that the there is a deficit of planets withsizes between 1.5 < Rp/R⊕ < 2.0 (Fulton et al., 2017).The mass distribution (either Mp or Mp sin I) peaks at 0.025MJ ≈ 8M⊕ and at 1.2MJ,with a higher cumulative fraction of planets within the larger mass range. The mass andsize distributions both have peaks that are separated by a valley at 0.13MJ ≈ 40M⊕ and0.4RJ ≈ 4R⊕, respectively, which seems to distinguish the giant planets from all other types.When the mass, size and orbital period distributions are considered altogether, small planets atshort orbital periods are more common than giants planet at those same periods (Howard et al.,2010; Mayor et al., 2011). Furthermore, Jupiter-mass planets are divided into short-period andlong-period groups, with a gap between 20 and 100 days. The planets with a measured massand/or radius are shown in the size versus orbital period plot in Fig. 1.1. The planets with solelya measured size are indicated in open circles, while the colored points indicate different rangesof mass. Earth-like planets would correspond to red points and sizes . 0.1RJ. The density ofbigger and more massive planets is easier to determine than smaller planets in the same periodrange. Nonetheless, a combination of observations and planet interior models suggests thatrocky planets have sizes Rp < 1.6R⊕ or < 0.14RJ (Rogers, 2015).Exoplanet surveys have shown that about 40% of the planets reside in systems with multiplemembers (Np ≥ 2), which, as previously noted, means that 20% of host stars have two ormore planets. Table 1.1 presents a summary of the number of stellar hosts with a given planetmultiplicity. The different values, corresponding to two different exoplanet databases, dependon the planet’s status and also whether all candidates suspected to orbit a single star have beenconfirmed. Low eccentricities (Wright et al., 2009) and planets with small sizes (Latham et al.,2011) characterize a significant number of multiplanet systems, or multis. Giant planets havealso been observed in multis, but at a lower occurrence rate. More recently, Weiss et al. (2018)found that the planet sizes and spacings in multiplanet systems are correlated, meaning thatneighboring planets are likely to have similar sizes and have regular spacings.Although there are similarities among multis, there are also clear differences, particularly inthe orbital distributions. Fig. 1.2 shows the period ratio distribution for systems6 in the EPE.All possible planet combinations within systems are included. The period ratio distributionprovides a way to test formation theories because disk migration models predict that planet5The upper limit in mass used by the EPE for planets is 60MJ, based on a revised density-mass relationshipfor planets, brown dwarfs and stars. Gas giant planets and brown dwarfs seem to follow the same relation, whichhas a clear cutoff at 60MJ (Hatzes and Rauer, 2015). Brown dwarfs are typically thought to have a lower masslimit of 13MJ based on the deuterium fusion limit (Burrows et al., 2001).6Confirmed planets and planet candidates.7No. of stellar hosts with given planet multiplicityMultiplicity Confirmed Confirmed + candidateNASA EPE EPE1 2174 2193 42912 399 419 5503 136 139 1524 51 46 525 20 16 166 6 8 97 1 1 18 1 1 1Total 2788 2823 5071Table 1.1: Multiplicity of extrasolar planetary systems. The second column shows the numberof systems with a given multiplicity, according to the NASA Exoplanet Archive (2018), accessedin May 2018. The third and fourth columns correspond to the EPE (Schneider et al., 2011),accessed in May 2018.pairs should tend to be trapped in or near mean motion resonances (MMRs), particularly first-order resonances (Goldreich and Tremaine, 1980; Lee and Peale, 2002). Dynamically importantcommensurabilities are indicated in Fig. 1.2 with dotted lines. There are indeed weak featuresfor period ratios in or near first-order MMRs (Lissauer 2011), i.e., commensurabilities of theform k + 1:k, which means that the inner planet completes k + 1 orbits when the outer planetcompletes k orbits. Moreover, Fig. 1.2 shows that a pattern seems to surface near the 2:1and 3:2 MMRs, where at slightly shorter period ratios there is a count deficit, followed by anoverabundance at slightly larger commensurabilities. Some recent studies have been able toreproduce this pattern using either the disk migration theory, but accounting for planetaryeccentricity (Goldreich and Schlichting, 2014), or accounting for interactions with the planetes-imal disk (Chatterjee and Ford, 2015). Nevertheless, the majority of planet pairs are not in aresonance. Altogether, this suggests that disk migration may be occurring, but may not be theonly (or even the dominant) mechanism for the existence of planets on short orbital periods.In the sample of confirmed and candidate planets, there are six pairs with period ratios closeto 1 ± 0.06, which suggest that either these pairs are strongly coupled or that they share thesame orbit. Three of the pairs have been confirmed: Kepler-132 b and c; Kepler-271 d and b(Lissauer et al., 2012; Rowe et al., 2014; Morton et al., 2016); and Kepler-1625 b and Kepler-1625 b I (Teachey et al., 2018). The Kepler-132 and Kepler-271 stars might have a bound stellarcompanion, and therefore the pairs might not orbit the same star. In the case of Kepler-1625b I, the transit signature is thought to be the first detection of an exomoon (extrasolar moon)(Teachey et al., 2018; Teachey and Kipping, 2018; Heller, 2018). The status of the other threepairs (K06242.02 and K06242.03, K00521.01 and K00521.02, K02248.04 and K02248.01) remainunknown at this time. In Sect. 3.4, we will discuss 1:1 MMRs further.8Figure 1.2: Exoplanet period ratio distribution of all possible combinations among the planetsin their respective system (black line) and with the nearest neighbor (in blue). The labelPout/Pin indicates that the period ratio shown is that of an outer planet to that of an innerplanet, and therefore it is always larger than 1. Confirmed and candidate planets are includedin this plot. Planets seem to avoid exact first-order MMRs, but do concentrate at slightly longerperiod ratios.Figure 1.3: Measured eccentricities versus orbital periods of detected exoplanets (data obtainedfrom Schneider et al., 2011). The different colors indicate the measured mass of each planet inJupiter masses, as shown in the colorbar. There is a significant fraction (24%) of planets withnearly zero eccentricity, as shown on the right histogram.Eccentricity and inclination are two orbital elements that can also be constrained from observa-tions, although they require more detailed observations and modeling than for determining theperiod. The majority of the measured eccentricities were obtained with the RV technique. How-9ever, it is also possible to estimate the eccentricity of transiting planets from transit durationsand TTVs, if the period, stellar size and impact parameter are known (Winn and Fabrycky,2015). Fig. 1.3, shows the exoplanet eccentricity distribution as a function of orbital period.The different colors indicate the planetary mass in Jovian units. A significant fraction of exo-planets are consistent with circular orbits (right panel), with smaller planets tending to havelower eccentricities (Wright et al., 2009; Mayor et al., 2011). Giant planets (green to dark blue)have a wider range of eccentricities, especially those with long orbital periods. The eccentricitydistribution spans from 0 to 1 and the typical eccentricity decreases with decreasing period.This trend may result from a combination of physical phenomena, along with RV detectionbias against e & 0.6, which is due to poor data sampling at periastron when eccentricities arelarge (Cumming, 2004). In Fig. 1.3 there is a sharp edge at about 1 day. For shorter peri-ods, all eccentricities are zero,7 regardless of their mass, which is usually attributed to tidalcircularization (Ogilvie, 2014). Fabrycky and Tremaine (2007) suggested a Kozai cycle8 withtidal friction during phases of high eccentricity to explain the shape of the upper boundaryin the distribution. In this mechanism, planets with an inclined companion are driven to ahigh mutual eccentricity by Kozai cycles, and suffer periodic close approaches with the star.During the close approach, the tidal forces cause the planet to heat, which is then dissipatedgradually. The energy dissipation decreases the eccentricity and semi-major axis (and thereforethe period) over time.The determination of the orbital inclination of exoplanets is of importance, mainly for RVdetection and mass estimation. Furthermore, the mutual inclination of planets in multis isof interest to planet formation theory. The low mutual inclinations in the Solar System areattributed to its formation from a rotating disk. Neither the transit nor the RV techniquesare particularly sensitive to mutual inclinations. Fabrycky et al. (2014) found that the Keplermultiplanet systems have mutual inclinations smaller than a few degrees, and are consistent witha Rayleigh distribution9 with σI ≈ 1.8◦. Several observational and dynamical studies comparingthe mutual inclinations of the Kepler multiplanet and single-planet systems observed differencesamong both distributions, in which the best fit for the multis was inconsistent with the best fitfor the singles (Lissauer et al., 2011b; Johansen et al., 2012; Hansen and Murray, 2013; Ballardand Johnson, 2016). As result, the singles are thought to belong to two different populations:some part of a flat system in which only the innermost planet transits, while others are part7Except for PSR J1719-1438 b, with a mass similar to Jupiter, an orbital period P = 0.090706293 d and aneccentricity e < 0.06 (Bailes et al., 2011). PSR J1719-1438 b is thought to be a white dwarf rather than a planet,due to its density.8The Kozai cycles occur when the orbit of a binary system is perturbed by a third distant object, which causesthe libration of the argument of pericenter. This libration is observed as a cyclic exchange of the eccentricityand inclination.9The Rayleigh distribution has the functional form:f(x;σ) =xσ2e−x2/2σ2being σ the mode of this distribution.10of a low multiplicity system or a system with much higher mutual inclinations (Lissauer et al.,2011a).The orbital inclination with respect to the stellar spin axis has been determined for a fewmultiplanet systems using the Rossiter-McLaughlin effect.10 Most of the measured stellar spinaxes are aligned with the normals to the planetary orbital planes (Fabrycky and Winn, 2009;Southworth, 2011). Yet for several systems (e.g. HAT-P-07, KELT-17, Kepler-63, WASP-08,WASP-33), there is a significant spin-orbit misalignment, which might be the outcome of awide range of processes, such as early misalignment of the protoplanetary disk with respect tothe stellar equator (Batygin, 2012; Crida and Batygin, 2014; Lai, 2014; Bate et al., 2010), puredynamical interactions between the planets (Kaib et al., 2011; Batygin et al., 2011; Boue´ andFabrycky, 2014a,b) or intrinsic stellar processes (Rogers et al., 2012, 2013).1.1.2 Hot and warm JupitersThe giant planet population, i.e., planets with 0.1 < Mp/MJ < 60 and Rp > 0.5RJ, is typicallysub-divided into period ranges, which are then related to a temperature adjective (cold, warmand hot). This temperature association comes from the expected equilibrium surface temper-ature of the planet, T , which is a function of stellarcentric distance, r, with a proportionalityT ∝ r−1/2 or T ∝ P−1/3, in terms of the orbital period. A hot-Jupiter (HJ) is a giant planetwith an orbital period P < 10 d, with an equilibrium surface temperature of about 1300 K.Following this reasoning, warm-Jupiters (WJs) are expected to have lower temperatures thanHJs because their orbital periods are longer, 10 < P < 200 days, and hence they are furtheraway from their host star. The rest of the giant planet population is referred to as cold-Jupiters.HJs are uncommon. Their occurrence rate in FGK stars in the solar neighborhood is ≈ 1%Wright et al. (2012), while in the Kepler data it is 0.5% (Howard et al., 2012).Different studies using RVs and TTVs have found that HJs, in contrast to WJs, seem tolack close companions (Wright et al., 2009; Steffen and Agol, 2005; Latham et al., 2011; Steffenet al., 2012; Huang et al., 2016). A key word here is “close”, which refers to smaller planets withsimilar orbital periods to either the WJs or HJs. The only known HJ with a close companionis WASP-47b (Becker et al., 2015), which is flanked by two mini-Neptunes. The most recentstudy to determine the frequency for which WJs and HJs have companions (Huang et al., 2016)investigated systems with planets that transit interior to the Kepler sample of gas giants withorbital periods within 0.5 and 200 days. A planet was deemed to be a gas giant if its radiuswas 8 < Rp/R⊕ < 20. Huang et al. found that the companion fraction of WJ and HJ aremutually exclusive and confirmed that HJs typically do not have close inner companions. In10The Rossiter-McLaughlin effect is the change to a rotationally broadened stellar spectral line due to theplanet blocking different portions of the star at different times during the transit. The shape of the spectral linesprovides information on the direction of the planet with respect to the rotation of the star.11comparison, half of the WJ have close companions. This led the authors to suggest that WJand HJ are formed by two different mechanisms, in which the latter might represent in-situformation. It must nonetheless be noted that planetary frequency does drop off rapidly withina few days, so other effects might be at play to shape the distributions. When companions at5 < a < 20 au are analyzed (Bryan et al., 2016), the fraction of outer companions of WJs with1 < Mp/MJ < 20 is lower than that of HJs. Both observations, low fraction of close companionsand high fraction of distant companions, point to a disruptive formation or evolution mechanismfor HJs.Studies of the stellar host metallicity and giant planet occurrence have found that the detectedgiant planets often orbit metal-rich hosts (Gonzalez, 1997; Santos et al., 2004; Fischer andValenti, 2005). The stellar metallicity is usually used as proxy for the metallicity of the pro-toplanetary disk, and therefore, the correlation between giant planets and stellar metallicity isinterpreted as an indication that giant planets are formed in protoplanetary disks with a highsolid-to-gas ratio. A correlation between stellar metallicity and orbital eccentricity has alsobeen found. Dawson and Murray-Clay (2013) determined that giant planets around metal-richstars have higher eccentricities than those around metal-poor stars. They suggest that metal-rich stars had protoplanetary disks that were rich in solid material and consequently formedmore planets, which could have then engaged in planet-planet interactions.In addition, the metallicity and planet occurrence rate do not correlate for all planet sizesand orbital periods (Petigura et al., 2018). Warm super-Earths (10 < P < 100 days and1.0 < Rp/R⊕ < 1.7) have a near-constant distribution over metallicities, from 0.4 to 2.5 timessolar values. Warm sub-Neptunes, on the other hand, double their occurrence rate over thesame metallicity range.11 Metallicity seems to correlate with the occurrence rate of planetswith P < 10 days and with planets that have sizes R > 1.7R⊕ (Petigura et al., 2018). Theinterpretation of these correlations remains unclear.The orbital period distribution of short-period giant planets is still a matter of debate, partic-ularly the existence and significance of the 3-day pileup. This feature was first observed in theRV data (Wright et al., 2009), but later found to be absent when taking into account the Ke-pler population (Howard et al., 2012). However, a recent study of the RV signal among Keplergiant planets, with false-positives removed from the sample, found that the occurrence rate ofthese planets in the orbital period range 1 < P < 400 days decays by an other of magnitude,with the highest occurrence at 400 days (Santerne et al., 2016). This occurrence distributionhas an important deficit between 10 to 40 d. This deficit could be interpreted as a pileup atabout 4 days, which is a prediction of HJ formation through eccentricity excitation plus tidalcircularization.11Meaning that at a metallicity 2.5 solar the sub-Neptune occurrence rate is twice than that at 0.4 solar (seeFig. 10 of Petigura et al., 2018).121.1.3 Systems with Tightly-packed Inner Planets (STIPs)The Systems with Tightly-packed Inner Planets (STIPs) are the main focus of this thesis. Theseare multiplanet systems with short orbital periods, ranging from 1 to 100 days. The orbitalperiod and multiplicity of these systems is not well constrained. The inner orbital periodcutoff is real, while the outer one is still debatable and dependent on observational limits andchallenges. Here, “tightly-packed” means that the planet period ratios are within 1 and 3(Ford, 2014). These small spacings are likely to affect the stability of the systems (discussed inmore detail in the next subsection). The minimum multiplicity for a system to be considereda STIP is not well-established and hence one needs to pick a working definition. We requirethat the multiplicity of a STIP be Np ≥ 3, because their dynamical evolution is richer than atwo-planet system, while their stability state cannot be predicted analytically. The formationand dynamics of STIPs are a challenge that increases in complexity with planet multiplicity.The fraction of confirmed planetary systems that are STIPs based on the above definitions isabout 5 to 6%. This fraction does not change significantly if the outer period cutoff is extendedto larger periods, as it is dominated by shorter orbital periods. In literature, the period cutoffis at 100 days for STIPs, which will be used throughout this work to avoid confusion with otherstudies. We further require that there are at least three planets with orbital period shorterthan 100 days. These two requirements select those systems that are difficult to reconcile withplanet formation theories based on our Solar System (described in Sect. 1.2).While STIPs are not limited to planets with a particular range of sizes and masses, the majorityof planets in known STIPs are either super-Earths or mini-Neptunes, i.e., planets with masseswithin 1.0 to 10 M⊕. Hence, the composition of most of these planets is likely dominatedby rock, ice and/or water, but not gas, even for low density planets. Moreover, formationmechanisms should be able to reproduce the mass-radius relationship and orbital distributionof STIPs.As mentioned earlier, two main formation mechanisms for STIPs have been suggested: in-situ assembly/formation and planet migration. Planet migration is the favored theory but ithas a few apparent inconsistencies. For example, we cannot reliably explain the direction ormagnitude of planet migration (Goldreich and Tremaine, 1980; Lin and Papaloizou, 1986; Ward,1997; Paardekooper and Mellema, 2006; Paardekooper et al., 2010, 2011), which depends ona large number of disk parameters. Neither can we explain why migration stops at particulardistances (Ida and Lin, 2004a,b; Alibert et al., 2005; Masset et al., 2006). Furthermore, wecannot easily discern from observations whether the inner regions of disks have a differentmetallicity from their host star, e.g., due to the local enhancement of solids, which woulddirectly impact the ability of a disk to form planets at short orbital periods. Planet migrationmight explain some STIPs and short-period giant planets, but not all. Alternative formationscenarios need to be explored, particularly those that address planetary diversity using a single13framework.1.1.4 Stability of extrasolar planetsDetermining whether a planetary system is likely to be stable for long timescales is of broad in-terest, both for exoplanets and the Solar System. Several studies have looked into this problem,but so far there is not a straightforward answer.A pair of planets in a two-planet system with initially circular-orbits will always be stable ifthey are separated by a dimensionless distance ∆ > 2.4 (µ0 + µ1)1/3 (Gladman, 1993), where µ0and µ1 are the mass ratios of the inner and outer planet relative to the star, respectively. Thesecond planet’s semi-major axis is thus a1 = a0(1 + ∆), where a0 is the semi-major axis of theinner planet. Planet separations can also be described in units of mutual Hill radii (Smith andLissauer, 2009; Lissauer et al., 2014a; Marzari, 2014; Pu and Wu, 2015), where the mutual Hillradius RmH = 0.5 (a0 +a1) [(µ0 + µ1)/3]1/3. This gives η ≡ (a1−a0)/RmH, for which two-planetstability requires η & 3.46. For example, a two Jupiter-mass planet system around a solar-massstar must have planet separations of 3.5RmH or greater.For higher multiplicities, there is no analytical separation limit. Usually, only the typicaltimescale for planets to become orbit crossing in a given system can be determined statistically(e.g., Obertas et al., 2017). This timescale depends on the initial η and the number of adjacentplanets with that η (e.g. Chatterjee et al., 2008). For Earth-mass planets with Np & 3, a mutualHill radius separation of 10 allows long-term stability (Smith and Lissauer, 2009). Planetmultiplicities (Table 1.1) and spacings are relevant to STIPs. Few STIPs have planets withη < 10, suggesting that the known planets should be stable over the lifetime of most stars(Lissauer et al., 2014a; Obertas et al., 2017). In a similar study, Pu and Wu (2015) found throughN-body simulations that the dynamical instability spacing threshold for planetary systems withNp ≥ 4 is 10 mutual Hill spacings for circular orbits and 12 for those with e ' 0.02. Theycompare their findings with the observed Kepler exoplanets’ mutual Hill radii, and find thatthe typical separation is 12, consistent with the limit determined from simulations. However,secular interactions can lead to the disruption of a system, even if the planetary spacing alonesuggests long-term metastability. A classic example is the evolution of Mercury in the SolarSystem, which has a small but non-negligible probability of being driven to orbit crossing withVenus on 5 Gyr timescales (Laskar, 1994).Statistical stability studies of known STIP analogues with different multiplicity has shown thatthese systems are metastable, or prone to decay (Volk and Gladman, 2015), with equal fractionsof systems going unstable in each decade of integration. This metastability implies that, overtime, the multiplicity of an initially more populated system decreases as a result of instabilityand collisions. Pu and Wu (2015) reach a similar conclusion, and suggest that STIPs were more14tightly-packed and evolved to their current configurations through instability. This might meanthat Kepler planets formed in a more dissipative environment than the terrestrial planets inthe Solar System.Because many STIPs have high multiplicity, we must ask whether a large fraction of the systemswith only two or three planets are decay products themselves.12 There are also observabilityconsiderations; in particular, the presence of outer perturbers can affect the observed planetmultiplicity of transiting systems, which can have a bearing on how we interpret planet-starmisalignment (e.g. Winn et al., 2005; Kaib et al., 2011; Boue´ and Fabrycky, 2014a,b), at leastin part.1.1.5 Minimum Mass Solar Nebula (MMSN) vs. Minimum MassExtrasolar Nebula (MMEN)The determination of protoplanetary disk densities, compositions and profiles is of great im-portance to planet formation models. The Minimum Mass Solar Nebula (MMSN) is a diskmodel that was derived by spreading the masses of the planets into annuli centered on theircurrent orbits and then augmenting those masses to restore the composition to solar values(Weidenschilling, 1977b; Hayashi, 1981). A commonly used version of this model isσsolids = 3.7× 102Fdisk Zrel( a0.2 au)−1.5g cm−2,where Fdisk and Zrel are the total disk mass and metallicity relative to the MMSN (Chiang andYoudin, 2010). This construction is usually used to establish order of magnitude calculationsdue to its intrinsic uncertainties, which neglect the possibility that planets did not form in theircurrent locations.With the discovery of exoplanets, this approach could be compared with other systems. Chiangand Laughlin (2013) estimated a Minimum Mass Extrasolar Nebula (MMEN), which is the solar-metallicity disk of gas and solids out of which the super-Earths uncovered by Kepler could haveformed, if planet formation were 100% efficient and orbital migration were negligible. Theauthors employed 1925 planet candidates with radii R < 5R⊕ and P < 100 d from Batalhaet al. (2013). In this model, each planet is assigned a surface densityσsolids,i ≡ Mi2piai∆ai=Mi2pia2i,which is the surface density required to form the known planets in-situ. The mass-radiusrelationship Mi = (Ri/R⊕)2.06M⊕ (Lissauer et al., 2011b) and ai = (P/yr)2/3 were used to12In planet building through planetesimal accumulation and then giant impacts, this is always true to someextent. Here we are referring specifically to otherwise fully built planetary systems that achieve instability inless than approximately 1 Gyr.15obtain the surface density of solids:σsolids = 6.2× 102Fdisk( a0.2 au)−1.6g cm−2,with a corresponding gas surface density about 200 times higher:σgas = 1.3× 105Fdisk( a0.2 au)−1.6g cm−2.These surface densities are related by Zsolar×Zrel = 0.0049, where Zrel is the fraction of metalsthat could have condensed as solids in the hot inner disk (Zsolar = 0.015). We point out thatthe MMSN given here has been adapted from Chiang and Youdin (2010) to match the scalingof the MMEN.Again, these surface densities have intrinsic uncertainties and therefore, their use is mainly toprovide order of magnitude estimates for formation models. For example, using the MMEN,the typical coagulation times for close-in rocky planets can be as short as 104 to 107 yrs forRp = 1− 5R⊕ within a = 0.1− 0.5 au (see Eq. 14 from Chiang and Laughlin, 2013).The MMEN depends strongly on the mass-radius relationship, the chosen planet sample andthe temperature profile. It does not account for different compositions throughout the disk orpossible density profile changes. Nonetheless, we will use the MMEN in Chap. 3 as a baselinefor gas surface densities.1.2 Planet formationA number of formation mechanisms have been proposed and developed to explain the originof both rocky planets and gas/ice giants. While there are many parts of the planet formationprocess that are believed to be understood, there are also many outstanding questions, partic-ularly regarding planetesimal formation. In the case of giant planet formation, there are twomain accepted ideas: the core accretion model and gravitational instability. In this section,we will only briefly describe the formation of rocky planets and the core accretion model forforming giant planets. Giant planet formation by disk instability is thought to be most relevantfor wide (50-100 au) orbit planets (Boley, 2009). Because we are principally concerned withSTIPs and planets on moderate orbits, we do not discuss this mechanism further and refer thereader to Helled et al. (2014) for a review.1.2.1 Formation of rocky planetsBriefly, the formation of a rocky planet starts with the settling of interstellar dust particlesand disk condensates to the mid-plane of the protostellar (or circumstellar) disk, collisionally16growing from micron to planetary scales (13 order of magnitude in size) (Armitage, 2007). Ra-diogenic dating of the Solar System materials and observations of protoplanetary disks aroundstars of different ages constrain the bulk of planet formation to occur in a few Myr (Wadhwaet al., 2007), followed by additional evolution over 107 to 108 yr. Select stages of planetarygrowth and potentially dominant processes during that stage are highlighted in Table 1.2.Stage of growth Physical process Diameter range Time scale(m) (yr)Dust to pebbles Settling and radial migration 10−6 − 10−2 103?Pebbles to rocks Uncertain mechanism 10−2 − 103103 − 106?Pebbles to rocks andplanetesimalsPairwise collisional growth 102 − 104Large planetesimals Runaway growth (σv . vesc) 104 − 105 103 − 104Large planetesimalsto planet embryosOligarchic growth (σv & vesc) 105 − 106 105Planet embryos toterrestrial planetsChaotic growth 106 − 107 106 − 108Table 1.2: Summary of the growth sequence to form rocky planets starting from dust particles.The physical processes and typical growth timescales are also included. The question marksindicate timescales that are not well defined. Adapted from Perryman (2011).First, we introduce two important physical phenomena affecting the interactions of solid objectswith the gas and other particles: aerodynamical drag and gravitational focusing. Each acts ondifferent size/mass regimes. Aerodynamic drag (gas drag) is due to the relative velocity of asolid body moving through gas. The gas in a protoplanetary disk has partial radial supportthroughout most of the disk due to an outward pressure gradient. This causes the gas totypically orbit more slowly than is necessary for solids, the latter of which do not directly feelthe same pressure gradient. The drag force is given byFD = −12CD(pis2)ρvv,where ρ is the gas density, and s and v are the particle’s radius and velocity relative to the gas,respectively. The drag coefficient, CD, depends on the gas mean free path13 λ, the particle’sshape (usually assumed to be spherical) and its size s (see Armitage, 2007, for more details).Depending on the local disk conditions, solids with sizes around the millimeter to decameterscales can experience rapid inward radial drift. For example, bodies with R ' 1 m migratefrom 1 au all the way to the star in about 100 yr, for a MMSN. Gas drag places an apparentbarrier to the growth beyond 1 meter scales (or smaller) due to the expected inward migrationtimescales, commonly known as the “meter-barrier” (Weidenschilling, 1977a).13For a MMSN, λ ∼ 0.01 m (Perryman, 2011).17While gas drag seemingly frustrates the growth of planetesimals, once formed, the collisioncross-section of an object can be enhanced by its gravity. The effective cross sectional radius,b, is given byb = R(1 +v2escσ2v)= RFg.Here vesc =√2GM R−1 is the escape velocity of the massive object with mass M , and radiusR, and σv is the velocity dispersion between objects. Fg is known as the gravitational focusingfactor and it is important for objects with R & 10 km.14 A massive body, such as a planetaryembryo, that is immersed in a swarm of planetesimals with surface density Σp will have a massgrowth rate given byM˙ =dMdt=12Σp Ωp piR2Fg,where Ωp is the orbital frequency (Armitage, 2007). The growth of an embryo is then a functionof vesc/σv; if the planetesimals have a smaller velocity dispersion than the embryo’s escapevelocity, then M˙ ∝ R4 and the embryo grows fast, or “runs away”. On the contrary, if vesc < σvthen Fg → 1 and M˙ ∝ R2. This means that the embryo would grow at a slower pace thanin the runaway stage but still faster than smaller objects. The latter is known as “oligarchicgrowth” because it gives opportunity to less massive planets to catch up.As Table 1.2 indicates, the dust condenses from the protostellar disk as it cools down, followedby growth while it rapidly settles to the mid-plane of the disk. Once in the mid-plane, theparticles drift inward due to gas drag until they reach a terminal velocity and in the meantime,grow through electrostatic and mechanical sticking. When the silicate dust grains reach amillimeter in size, the collisions do not continue to form aggregates because they bounce offeach other. This is known as the “bouncing barrier” (Zsom and Dullemond, 2008; Gu¨ttler et al.,2009). The growth mechanism from cm-sized particles to full planetesimals is unknown, due toa combination of radial drift and particle bouncing. Different processes have been suggested toexplain the rapid growth to planetesimal scales; we refer the reader to Armitage (2007).Once objects reach sizes around 1 km, the planetesimals become decoupled from the gas onshort timescales, allowing the dominant growth mechanism to be pairwise collisions. Initially,the presence of the gas keeps the velocity dispersion of the planetesimals low, and thereforethose with a bigger size run away. When the bigger aggregates of planetesimals are about 100km, they stir up the orbits of the smaller objects and the velocity dispersion increases, reducingthe gravitational focusing and entering the oligarchic growth phase.The growth of a planet embryo is limited by the availability of planetesimals at a disk distance14This threshold depends on the assumed σv, which roughly scales as the product of the eccentricity, e andthe circular velocity, vc, of the planetesimals at a given stellarcentric distance, such that: σv ∼ e vc ∝ e r−1/2,implying that at larger stellarcentric distances the gravitational focusing becomes important for smaller objectsthan closer to the star.18r, within an annulus 2∆r. The maximum mass of an embryo, or isolation mass, is thenMiso = 4pi r∆r σρ(r) = 16pi r2(Miso3M∗)1/3σρ(r),=√[16pi r2 σρ(r)]33M∗, (1.1)where M∗ is the stellar mass, and σρ(r) is the surface density of solids at r. These isolated em-bryos can still affect the orbits of the neighboring planetesimals and other embryos, particularlyif the mass of the embryo is comparable to the total mass of planetesimal disk. The cumulativeeffect is a chaotic environment that could lead to collisions and merging of the embryos or toprotracted accretion of planetesimals. The final assembly of rocky planets continues with theembryo mergers, until the orbital spacing results in a quasi-stable configuration, which takesa few Myr to potentially over 100 Myr under standard planet formation models (Chambersand Wetherill, 1998; Goldreich et al., 2004; Moorhead and Adams, 2005; Kokubo et al., 2006;Kenyon and Bromley, 2006).1.2.2 Giant planets: core accretion plus gas capture modelThe core accretion model, also known as “core nucleated instability”, posits that giant planetsform in three distinct stages. First, there is an initial core formation of rocky planets, asdescribed above. Second there is a protracted period of slow gas and planetesimal growth,which transitions to a runaway gas capture phase when a critical mass threshold is reached.A rocky planet core keeps accreting both disk solids and gas as available within its gravitationalreach. When there are no more solids in the vicinity, the planet opens a gap in the solid surfacedistribution, but not necessarily in the gaseous disk. The planet steadily accretes gas, but thecapture onto an envelope depends on the cooling and contraction of the accreted gas as muchas the planetary mass. Under the right conditions, the planetary mass keeps increasing untilit reaches a critical mass at which gas runaway accretion starts. The critical mass definitiondepends of the time dependence of the model used. In static models, where there is a steadyaccretion of planetesimals, the critical core mass is that above which the gaseous envelope cannotsupport hydrostatic equilibrium (Mizuno, 1980; Stevenson, 1982; Rafikov, 2006, 2011). On theother hand, if the model is time-dependent, a critical core mass is defined when the envelopegrows superlinearly. This occurs when the self-gravity of the envelope becomes significant, i.e.,when the envelope has as much mass as the core, and the thermal disequilibrium increases(runaway cooling, Pollack et al., 1996; Ikoma et al., 2000). The growth stops either when a gapis opened in the gaseous disk or when the disk dissipates, starting the isolation stage of growth.The core accretion model is thought to be more efficient in regions of the disk with a higherlocal solid-to-gas ratio, particularly beyond the ice line. The temperature gradient within a19disk is ∝ (r/1au)−p, where p is a positive exponent, and therefore, the temperature decays withstellarcentric distance, eventually allowing additional solids to condense from the gas.In principle, if rocky planets are formed efficiently and fast within the lifetime of the gaseousdisk, regardless of their stellarcentric distance, it facilitates the subsequent gas acquisition. Thisopens the door to the possibility of in-situ formation of gas giant planets. Simplistically, theamount of gas accreted onto a planet depends on: (1) availability of gas in the disk; (2) thegravitational well of the planet; and (3) cooling rate of the planet to hold the envelope. Theavailability of gas is related to the disk lifetime and gap opening, while the gravitational welldepends on the mass of the planet. The cooling rate depends on the heat transport of theenvelope, which can be affected by dust. We then could ask: can rocky planets form efficientlyin the innermost region of the disk before the gas disperses? If so, can they reach a criticalmass to start runaway gas accretion? Is orbital planetary metastability present when the gasis around? How many rocky planets can be formed?1.3 This thesisThroughout this work we seek to answer two main questions. Firstly, assuming that rockyplanets form efficiently in the inner region of the disk before the gas disperses, is metastabilitypresent when the gas is still around? If so, can critical cores form from this initial instability?Secondly, what dynamical effects originate from the presence of a planetary perturber outsideknown STIPs?. This thesis is divided into six additional chapters. In Chap. 2, we provide anoverview of the methodology used in the following two chapters, including a description of thenumerical methods, i.e., N-body integrators and secular code. The experimental developmentand analysis of the two main science questions are developed in Chap. 3 and Chap. 4. InChap. 3, we evaluate whether critical cores can be formed through the metastability of pri-mordial high-multiplicity STIPs in a gas-free and gas-rich environment; while, in Chap. 4, wedetermine dynamical effects that an outer planetary perturber has on STIPs, and how theyaffect the stability and observability of the inner system. We discuss the results of the numeri-cal simulations and analysis in Chap. 5, including limitations of each experiment with possibleways to improvement as future work, as well as the importance and implications of the results.In Chap. 6, we summarize the main take away points of each experiment. Finally, a list ofcomplementary and improved experiments as possible future work is provided in Chap. 7.20Chapter 2Methodology overviewWe designed a series of N-body simulations: 1) to explore the in-situ formation scenario of hot-Jupiters (HJs), in which giant planet formation is envisaged to be the outcome of planet-planetscattering followed by collisions and consolidation; and 2) to investigate the effects of the long-term gravitational perturbations, i.e., secular effects, of an outer giant planet on known STIPs.In both studies, numerical experiments are carried out using a large number of planetary systemrealizations that have initial conditions based on STIP prototypes with a planet multiplicityNp ≥ 6.As part of these studies, we evaluate system stability for any given STIP prototype by perturbingthe initial conditions for each realization, keeping initial planetary masses and semi-major axesfixed. The specifics of the initial conditions will be described in Sect. 3.1. The first experimentis divided into two stages, each with increasing dynamical complexity. Initially, planet-diskinteractions are ignored so that we can estimate the fraction of systems in which at least oneplanet with Mp ≥ 10M⊕ forms through consolidation. We then include a gaseous tidal dampingprescription in the N-body integrator to account for planet-disk interactions, as described inSect. 2.1.1. We run the realizations without planet-disk interactions with Mercury6, while15th-order Integrator with Adaptive Time-Stepping (IAS15) is employed whenever the gaseoustidal damping is used, on account of the non-conservative forces. I should emphasize that theIAS15 code is my personal implementation, in C language, of the algorithm developed by Reinand Spiegel (2015), which will be described in Sect. 2.1.4.It is uncertain whether the known STIPs have any outer companions. Should such companionsexist, they could affect the dynamics of STIPs through gravitational perturbations. To explorethis situation, we run N-body simulations of known STIPs, but include a Jupiter-like perturber.We use Kepler-11 (K11) and Kepler-90 (K90) as prototypes (see Sect. 4.1). The perturber’sdistance to the host star is varied from 5.2 au to 1 au to explore a range of secular dynamics. Forcertain system configurations, we calculate the theoretical secular precession eigenfrequenciesand determine whether there is strong inclination or eccentricity forcing. In addition, we obtainthe synthetic secular precession frequencies and compare them to the theory. Discrepanciesbetween the theoretical and synthetic spectra can reveal strong gravitational interactions, suchas mean motion resonances or non-linear secular resonances.212.1 Numerical simulationsIn this section, we describe the main integration codes and algorithms used in this work:Mercury6 (Chambers, 1999) and IAS15 (Rein and Spiegel, 2015). The physical forces in-cluded in the simulations, along with explanations regarding their relevance, are discussedin Sect. 2.1.1. The handling of close encounters and collisions is detailed in Sect. 2.1.2, whichis crucial to addressing planetary consolidation. Sect. 2.1.3 provides an overview of Mercury6’shybrid-symplectic integration algorithm, and finally, Sect. 2.1.4 details the IAS15 integrator,comparing my implementation to that of Rein and Spiegel (2015).2.1.1 Beyond Newtonian Gravity: additional forcesThe effects of general relativity (GR) become important for planets with short orbital periods,as it changes the precession of the planet’s longitude of pericenter, $. In all simulations, weinclude a correction for GR as defined by (Nobili and Roxburgh, 1986), which is modeled as adipole-like potential and only depends on the position of the body, i.e.,aGR = −6G2M2c2r4r. (2.1)With this approximation, GR can be modeled in both Mercury6 and IAS15. Some simulationsalso include gaseous tidal damping, i.e., damping due to the gravitational interaction betweenthe planet and the spiral waves launched by the planet. This interaction is prescribed here as:aTD = − 1τe(v · ρ)‖ρ‖2 ρ−(v · kˆ)τikˆ= −vρτeρˆ− vzτikˆ. (2.2)Here, we use cylindrical coordinates (as shown in Fig. 2.1), where ρ = x iˆ + y jˆ, and vρ and vzare the magnitude of the radial and vertical components of the velocity. This prescription suc-cessfully damps the eccentricity and inclination of planets in a way that is in overall agreementwith the results of Cresswell et al. (2007). Furthermore, this damping is similar to that usedby Papaloizou and Larwood (2000), but different by a factor of 2. This difference is ultimatelyabsorbed by the eccentricity and inclination damping timescales τe and τi, respectively. Thesetimescales control the magnitude of the damping acceleration, as shown in Eq. 2.2, and aregiven by1τe=KeKeτwave +∣∣∣vρvc ∣∣∣3 and1τi=KiKiτwave +∣∣∣vzvc ∣∣∣3 . (2.3)22Figure 2.1: Notation used in cylindrical coordinates.where vc =√GM∗/ρ is the circular velocity at ρ. Both timescales are related to the characteris-tic wave time, τwave, which is determined through tidal torque theory of planet-disk interactionsin isothermal disks (see Tanaka and Ward, 2004, Eq. 49). This wave time describes the meantimescale for the orbital elements of a planet to evolve due to the forcing by density waveswithin a gaseous disk; namelyτwave =(MpM∗)−1(σp a2M∗)−1(csaΩp)4Ω−1p , (2.4)where Mp and M∗ are the masses of the planet and star, respectively. The wave time alsodepends on the planet’s semi-major axis, a, and mean motion, Ω2p = GM∗a−3, often presentedas n in astrodynamics. The disk properties are included in the surface density, σp(ρ), and theisothermal sound speed at radius ρ, i.e.,cs(ρ) =√RTµwith molecular weight of the gas µ = 2.3 (in g/mol), which corresponds to the value for theinterstellar gas that was thought to lead to the formation of the Solar System. In IAS15,Ketwave = 0.00253 and Kitwave = 0.00075 are treated as constants (Cresswell et al., 2007).Consequently, as τwave changes according to the local disk conditions, the values of Ke and Kialso change.In practice, we do not use τwave directly. In IAS15, we introduce the timescale τlocal, which isin part a scaled version of Eq. 2.4, at 1 au. This reduces the number of calculations in the codeand allows us to include the vertical distribution of gas in the disk and the density decay withtime. The timescale τlocal depends on the local disk conditions and the planet’s position and23mass in a similar way to Eq. 2.4. We account for the vertical distribution of the gas, z, on thetorque by assuming a Gaussian profile for the density, with scale height H = cs/Ωp. We alsoassume that the disk density decays exponentially with time using a characteristic timescale λ,such thatτlocal(ρ, z, λ,Mp; t) = τ1auo(MpM∗)−1 ( ρ1au)q−2p+3/2ez2/2H2 et/λ, (2.5)where p and q are given by the chosen thermal and surface density profiles, respectively. Thefactor τ1auo is a scaled wave time at ρ = 1 au, given byτ1auo ≡(σo ρ2oM∗)−1(RTo ρoµGM∗)2(GM∗ρ3o)−1/2∣∣∣∣∣ρo=1auτ1auo =(ρ3oG5M3∗)1/2(RToµ)2σ−1o∣∣∣∣∣ρo=1au(2.6)τ1auo ≈ 0.078 days.To obtain the numerical value, we assume M∗ = 1M , ρ = 1 au, σo = 9.899 × 103 g cm−2and To = 300 K. For reference, at t = 0 the local damping timescale of an Earth mass planetlocated at 1 au and in the disk mid-plane (z = 0) is τlocal ≈ 71 yr. The gaseous surface densityprofile, from which σo proceeds, is the MMEN (Chiang and Laughlin, 2013)σp = σMMEN(ρ) = σo( ρ1 au)−q, (2.7)where q = 1.6 and σo is a scaled surface density at 1 au with value σg = 9.899× 103 g cm−2 forthe gas or σs = 47.21 g cm−2 for the solids. The thermal profile is set according toT (ρ) = 300( ρ1 au)−pKwith p = 0.6.The use of τlocal provides a convenient way of setting the damping timescale (through τo) forexperiment design. The convenience of Eqs. 2.5 and 2.6 is that they can be easily adapted fordifferent assumptions of disk density and temperature profiles.2.1.2 Close encounters and collisionsIn most N-body integrators, the interacting bodies are taken as point masses as long as thereare no close encounters or collisions. Close encounters occur whenever two bodies are withina threshold distance of each other. A collision, in it simplest form, occurs when the distancebetween two bodies is smaller than the sum of their “physical” radii (rij ≤ Ri +Rj), a strategy24that we also adopt for this work. Commonly, in simulations of planetary systems the closeencounter threshold is set to be an integer multiple of the planets’ mutual Hill radiusRmH =ai + aj2(mi +mj3M∗)1/3, (2.8)where a is the semi-major axis, and m and M∗ are the planetary and stellar masses, respectively.In practice, the multiple of RmH ranges between 1 and 3 and is a trade-off between resolving theencounter properly and minimizing computational time, particularly for the hybrid-integratorin Mercury6. In IAS15, a close encounter is flagged whenever the distance between two bodiesis rij ≤ 3RmH,15 which is evaluated at each substep.The planetary size is indispensable in determining whether a collision has occurred within asimulation. Mercury6 and IAS15 calculate the size of the planets with the provided bulk densityand mass. For most planets discovered by the transit method, only the size is known. On theother hand, for planets detected by the RV method, only the product of the mass times sin I isknown. About 300 planets have independent measurements of both radius and mass (Chen andKipping, 2017). As a consequence, the bulk density of the majority of planets is unknown. Forplanets with unknown masses or radii, several power-law mass-radius relationships have beensuggested, such asRp = R⊕(MpM⊕)ν, (2.9)where ν takes different values according to the data used to fit the function (Valencia et al.,2006; Lissauer et al., 2011b; Wolfgang et al., 2016; Chen and Kipping, 2017). The databaseexoplanets.org takes the value ν = 0.485 (Lissauer et al., 2011b), using a fit to the massesand radii of the planets in the Solar System from Earth to Saturn. Alternatively, Valenciaet al. (2006) who included interior models of Earth-like planets with different compositions andmantle-core ratios, derived ν = 0.267− 0.272 for masses 1 < Mp/M⊕ < 10.It is often desirable to infer the bulk density of a planet, which can also be estimated from themass-radius relationship:ρ = ρ⊕(MM⊕)1−3ν, (2.10)ρ = ρ⊕(3.33× 105)1−3ν ( MM )1−3ν. (2.11)Each expression is scaled either to the Earth’s (Eq. 2.10) or the Sun’s mass (Eq. 2.11), andρ⊕ = 5514 kg m−3 is Earth’s bulk density.15As discussed in Sect. 1.1.4, two planets need to be separated by a mutual Hill radii of η & 3.46 to bestable (Gladman, 1993). Therefore, tagging a close encounter with smaller separation increases the probabilityof identifying a pair of planets that might be involved in a collision later in the simulation. By choosing thetagging distance at 3RmH, we are being conservative in the early identification of close encounters.25Collisions in both codes are assumed to be inelastic. The new mass is determined as the sum ofthe colliding bodies, which is then used to determine the bulk density using Eq. 2.11. The newbody’s radius is then calculated from the mass and density. The position and velocity vectorsof the consolidated body are then estimated from the center of mass of the parent bodies. Eachtime that a collision occurs, the size of the body array16 remains the same but the order ofthe elements is changed. The lost planet is moved to the end of the array and assigned a zeromass, while the other planets are shifted accordingly one place lower in the array. As a finalstep, the total number of bodies is reduced by one, which shortens the loop for evaluating thenumber of planets, reducing the integration time. Other possible dynamical outcomes, besidescollisions, are the loss of planets due to scattering outside the system or due to collisions withthe star. In these two cases, the total number of interacting masses is reduced each time, aswell as the total mass available for planet consolidation. When a planet is lost into the star,the center of mass is affected by the redistribution of mass. On the other hand, if a planetreaches a threshold distance from the star that is considered unbound to the planetary system,the code stops following the orbit.2.1.3 Mercury6Mercury6 (Chambers, 1999) is primarily used to run the simulations that explore the feasibilityof critical core consolidation in the absence of planet-disk interactions. Mercury6 possessesfive integration algorithms: a second-order mixed-variable symplectic (MVS) integrator (Wis-dom and Holman, 1991; Wisdom et al., 1996); a general Bulirsch-Stoer (BS) (Stoer and Bu-lirsch, 1980; Press et al., 1992a); a conservative BS; a non-symplectic fifteenth-order integratorwith Gauss-Radau spacings (RADAU, Everhart, 1985); and a hybrid symplectic/BS integrator(Chambers, 1999). From the five different integrators within Mercury6, the hybrid integratorwas utilized due to its ability to accurately evolve systems for protracted periods of time andcapture close encounters.The hybrid integrator uses mixed-center variables conformed by heliocentric positions andbarycentric velocities (mixed-variable symplectic, MVS) for most of the evolution, but it changesfrom the MVS to the BS integrator to resolve/integrate close encounters. Should close encoun-ters not occur, then the system is integrated solely with the MVS algorithm.A symplectic algorithm conserves the phase-space (q, p)17, and is a numerical method for solvingHamilton’s equationsp˙ = −∂H∂qand q˙ =∂H∂p, (2.12)where H(p, q) = T (p) + V (q) is the Hamiltonian, and q and p are the canonical position and16This means that within the code, in C language, the bodies are defined in an array or vector, where eachbody is identified with an index.17Not to be confused with the power-law indices q and p given above26momentum coordinates, respectively. The total variation with time of a quantity z = z(qi, pi, t)for a system of N particles isdzdt=N∑i=1(∂z∂qidqidt+∂z∂pidpidt)=N∑i=1(∂z∂qi∂H∂pi− ∂z∂pi∂H∂qi)= Fz,where F is an operator defined asF ≡N∑i=1(∂H∂pi∂∂qi− ∂H∂qi∂∂pi). (2.13)The general solution of Hamilton’s equations (Eqs. 2.12) is given by the integration ofz(t)∫z(t−τ)dzz=t∫t−τFdt ⇒ z(t) = eτF z(t− τ), (2.14)where z(t) is either the canonical position or momentum at time t and τ is the timestep.In practice, it is common to divide the Hamiltonian into two or more components that can besolved analytically in the form of H = HA + HB + · · · , where HB is a perturbation to HA.When the Hamiltonian is separable, the respective operator is the sum F = A + B and thesolution changes toz(t) = eτ(A+B) z(t− τ), (2.15)= eτAeτB z(t− τ). (2.16)Numerically, the exponential in Eqs. 2.14 and 2.16 are expanded in Taylor series of τ :eτF = 1 + τF +τ2F 22+ · · · . (2.17)This approximation is exact to first or second-order depending on the expansion of eτ(A+B).Applying each operator sequentially gives a first-order integrator (Eq. 2.16). The integratorcan be made second-order by splitting the B operator into a half time-step and applying itbefore and after the A operator (Eq. 2.18). Because the B operator represents a perturbation,27this can be viewed as a form of a kick-drift-kick leapfrog algorithm,18z(t) = eτB/2eτAeτB/2 z(t− τ). (2.18)The expanded exponential in Eq. 2.15 is the same to first-order in Eq. 2.16 and to second-orderin Eq. 2.18. To be explicit,eτ(A+B) = 1 + τ (A+B) +τ22(A+B)2 + · · ·eτAeτB = 1 + τ (A+B) +τ22(A2 + 2AB +B2) + · · ·eτB/2eτAeτB/2 = 1 + τ (A+B) +τ22(A+B)2+τ34[(A+B)(A2 + 2AB +B2) + (B2 +BA−A2)A] + · · ·Recall that the operators A and B are not assumed to commute, i.e., AB 6= BA and thatA2 +AB +BA+B2 6= A2 + 2AB +B2.For a second-order integrator, the error between eτ(A+B) and eτB/2eτAeτB/2 is O(τ3). If τis small, then the difference between the analytic and numerical solutions is also small. Theenergy error does not build up because the method is symplectic.The Hamiltonian in the MVS integrator (Chambers, 1999) for an N-body system is H = H0 +H1 +H2, whereH0 =N∑i=0(p2i2mi− Gmim∗ri∗)−GN∑i=1N∑j=i+1mimjrij[1−K(rij)] , (2.19)H1 = −GN∑i=1N∑j=i+1mimjrijK(rij), (2.20)H2 = 12m∗(N∑i=1~pi)2. (2.21)Here, H0 is the unperturbed motion of the secondary bodies around the primary (e.g., planets18This means that the velocity, vi, and position, xi, are updated from the time ti to ti+1 by a “leapfrog”method. This involves giving a kick with half timestep to the velocity, drifting the position with this kick by afull step, and then using this drift to update the second derivative of x, i.e., ai, to ai+1, to finally update thevelocity to vi+1 using a second half timestep kick. This scheme is thus written as:vi+1/2 = vi + ai∆t2(kick),xi+1 = xi + vi+1/2 ∆t (drift),vi+1 = vi+1/2 + ai+1∆t2(kick),where the subindex i+ 1/2 indicates a half timestep variation.28around a star) if there are no close encounters. Perturbations due planet-planet interactionsare given by H1. Finally, because barycentric velocities are used, an additional Hamiltonian,H2, is needed to account for the kinetic energy due to the star. In Mercury6, H1 and H2, arealways less than H0 due to the dependence of H0 and H1 on K(rij):K(rij) =0 y < 0y22y2−2y+1 0 < y < 11 y > 1(2.22)andy =(rij − 0.1 rcrit0.9 rcrit), (2.23)which accounts for possible close encounters (see Chambers, 1999, for more details). The thirdterm in H0 (Eq. 2.19) equals zero if rij > rcrit19. When rij ≤ 0.1 rcrit, H1 → 0 and the mutualplanetary interaction is included in H0. At the same time, the integrator changes to the BS sothat it can resolve the close encounter and a possible collision.The second-order integrator for the MVS isz(t) = eτB/2eτC/2eτAeτC/2eτB/2 z(t− τ). (2.24)This means that the integration scheme follows these steps: (1) accelerates each body by planet-planet perturbations weighted by the factor K(rij) (the coordinates are fixed) with a step τ/2(Eq. 2.20); (2) the positions are then modified by a contribution from the momenta (which isfixed at this point) (Eq. 2.21); (3) the bodies are next accelerated by the star for τ and weightedby the term 1−K(rij), which is equal to 0 for a non-encounter or increasing to 1 as rij → rcrit,i.e., when a close encounter is about to happen, as expressed in Eq. 2.19. The scheme finalizeswith the repetition of step (2) and then (1).2.1.4 15th-order Integrator with Adaptive Time-Stepping (IAS15)Unlike Mercury6, IAS15 is a 15th-order implicit predictor-corrector integrator with adaptivetime-stepping (Rein and Spiegel, 2015). It is based on the 15th-order Runge-Kutta algorithm20of Everhart (1985), which uses Gauss-Radau spacings (Abramowitz and Stegun, 1964) to inte-grate a second-order derivative that depends on position, velocity and time: y′′ = F (y′, y; t).The algorithm adjusts the timestep to efficiently resolve the physical scale of the given systemand achieve convergence to machine precision.19rcrit is given in terms of Hill radii, i.e., RH = ap(Mp/3M∗)1/3, where Mp and ap are the planetary mass andsemi-major axis, and M∗ the mass of the star.20The Runge-Kutta algorithm better known is that of 4th-order, that is described in several books of numericalmethods (e.g. Abramowitz and Stegun, 1964; Hildebrand, 1987a; Press et al., 1992b).29Gauss-Radau quadrature is a scheme for integrating a function f(y, t) by approximating theintegral as a weighted sum of the function evaluated at predetermined spacings hn:1∫0f(x)dx =2n2f(0) +n−1∑k=1Wkf(hk) + E, (2.25)whereE =22n−1n[(n− 1)!]4[(2n− 1)!]3 f(2n−1)(ξ) (2.26)is the error of the Gauss-Radau quadrature of order 2n−1 and f (2n−1)(ξ) is the (2n−1)th deriva-tive of f evaluated in the range 0 ≤ ξ ≤ 1 (Hildebrand, 1987b). The number of Gauss-Radauspacings and their values depend on the order of the polynomial that is used to approximatethe function in the integral in Eq. 2.25. The method requires n points to fit polynomials ofdegree 2n− 1 or lower, and the number of free abscissas is n− 1, where the first evaluation isat h0 = 0 for the Radau-like quadratures.There are different ways to calculate the Gauss-Radau spacings for degree 2n− 1 polynomials.Abramowitz and Stegun (1964) present an analytic equation to find the optimal spacing hn bytreating them as roots of a combination of Legendre polynomials, which are the polynomialsused in Gauss-Radau quadrature:Pn (2h− 1) + Pn−1 (2h− 1)2h.In IAS15, the degree is 2n − 1 = 15, and therefore n = 8. As a consequence, the values of hnare the roots ofP8 (2h− 1) + P7 (2h− 1)2h= 6435h7 − 24024h6 + 36036h5 − 27720h4+ 11550h3 − 2520h2 + 252h− 8 = 0. (2.27)Table 2.1 shows the Gauss-Radau spacings to seven decimal places, and Fig. 2.2 shows thegraphical representation of Eq. 2.27.IAS15 uses a predictor-corrector scheme to integrate a second-order derivative, y′′ = F (y′, y; t).It predicts the position, y, and velocity, y′, at the next step, t, by first fitting a polynomial toy′′ at several computed spacings, hn, and then integrating the fit. With the predicted y and y′,the acceleration is determined at t to correct the values of y and y′. The process is iterated untilthe values of both y and y′ converge to machine precision. Under this scheme, the predictor isgiven by the polynomialy′′(h) ≈ y′′0 + g0 h+ g1 h (h− h1) + g2 h (h− h1)(h− h2)+ · · ·+ g7 h (h− h1) · · · (h− h7), (2.28)30Figure 2.2: Values of the Gauss-Radau spacings (in red) obtained from the roots of Eq. 2.27,which are used to numerically integrate a polynomial of degree 8.n hn0 0.01 0.05626262 0.18024073 0.35262474 0.54715365 0.73421026 0.88532097 0.9775206Table 2.1: Numerical values of the Gauss-Radau spacings to 7 decimal places. These values arethe roots of the Radau polynomial in Eq. 2.27 and represented graphically in Fig. 2.2.where y′′0 is the acceleration at the start of the step21, gk are constants to be determined withunits of acceleration, and the hns take the values of the Gauss-Radau spacings (Table 2.1). Forthis approximation, it is assumed that the derivative is smooth and well behaved, and as aconsequence, the second-order derivative can also be expressed as a polynomial of t:y′′(h) ≈ y′′0 + b0h+ b1h2 + · · ·+ b6h7; (2.29)y′′(t) ≈ y′′0 + a0t+ a1t2 + · · ·+ a6t7.Here, y′′(h) = y′′(t) for timestep dt, h ≡ t/dt, and bk ≡ ak dtk+1. The velocity, y′(h), andposition, y(h), can be found by integrating Eq. 2.29 once and then twice, with expressionsy′(h) ≈ y′0 + h dt(y′′0 +h2(b0 +3h2(b1 + · · · ))), (2.30)y(h) ≈ y0 + y′0 h dt+h2 dt22(y′′0 +h3(b0 +h2(b1 + · · · ))), (2.31)where y′0 and y0 are the velocity and position at the beginning of the step. To predict the21When h = h0 = 0.31position and velocity at t, the values of the bk constants are needed. The constants bk canbe calculated from gk by equating Eq. 2.28 and 2.29. On the other hand, the estimation of gkrequires that the acceleration, in this case y′′ = GM/r2 (for only point-like gravitational forces),is evaluated at each substep hn. Using the value y′′n = y′′(hn) in conjunction with Eq. 2.28 onefinds thath = h1 gives g0 =y′′1 − y′′0h1,h = h2 gives g1 =y′′2 − y′′0 − g0h2h2(h2 − h1) ,h = h3 gives g2 =y′′3 − y′′0 − g0h3 − g1h3(h3 − h1)h3(h3 − h1)(h3 − h2)and so forth. The notation y′′(hn) means that Eqs. 2.31 and 2.30 have been evaluated at thegiven substep hn. The values of bk are either taken to be bk = 0, such as at the start of thesimulation, or have been predicted from the previous step. In general,gk =y′′k+1 − y′′0 −k−1∑i=0gii∏n=0(hk+1 − hn)k∏n=0(hk+1 − hn)where k = 0, 1, ..., 6. (2.32)It is clear, that the process to find the values gk is iterative, with a loop over the spacings hk.After estimating gk, those values are used to update the bks according toc00 g0 + c10 g1 + c20 g2 + c30 g3 + c40 g4 + c50 g5 + c60 g6 = b0c11 g1 + c21 g2 + ... + c61 g6 = b1c22 g2 + ... + c62 g6 b2. . .......c66 g6 = b6(2.33)where the coefficients cjk are constant through all executions of the algorithm and, therefore,calculated just once using the relationcjk =1 for j = k,−hj cj−1,0 for j > 0 and k = 0,cj−1,k−1 − hj cj−1,k for k < j.(2.34)The position and velocity can be corrected by substituting the estimated values of bk intoEqs. 2.31 and 2.30, and using h = 1. The procedure from Eq. 2.28 to 2.33 is looped until the32convergence criterion is reached. Rein and Spiegel (2015) define this criterion to beδ˜b6 ≡ maxi |δb6,i|maxi |y′′i |< 10−16. (2.35)Here, δb6,i = bnew6,i − bold6,i is the difference between the value of b6 in the previous iteration andthat of the present one. The subindex i runs over the three vector components of all particlesfor both δb6 and the acceleration y′′. Eq. 2.35 takes advantage of the fact that the values ofbk vary slowly and that b6 is the smallest bk coefficient. The acceleration is usually dominatedby the b0 coefficient and for this reason, when machine precision is achieved, the value of b6 isinsignificant compared with y′′.The number of iterations for each step can be reduced by predicting the coefficient bk from theprevious values. This is done explicitly usingB0 = Q ( 7 b6 + 6 b5 + 5 b4 + 4 b3 + 3 b2 + 2 b1 + b0 ),B1 = Q2 ( 21 b6 + 15 b5 + 10 b4 + 6 b3 + 3 b2 + b1 ),B2 = Q3 ( 35 b6 + 20 b5 + 10 b4 + 4 b3 + b2 ),B3 = Q4 ( 35 b6 + 15 b5 + 5 b4 + b3 ),B4 = Q5 ( 21 b6 + 6 b5 + b4 ),B5 = Q6 ( 7 b6 + b5 ),B6 = Q7 b6.(2.36)In this case, Q = dtnew/dtold and Bk denotes the predicted values. This procedure reducesthe number of iterations, usually to just two, in the following timestep and accelerates theconvergence. With these Bk → bk, the gk are recalculated using the reversal of Eqs. 2.33 and2.34:d00 b0 + d10 b1 + d20 b2 + d30 b3 + d40 b4 + d50 b5 + d60 b6 = g0d11 b1 + d21 b2 + ... + d61 b6 = g1d22 b2 + ... + d62 b6 = g2. . .......d66 b6 = g6(2.37)with coefficients given bydjk =1 for j = k,h1 dj−1,0 = hj−11 for j > 0 and k = 0,dj−1,k−1 + hk+1 dj−1,k for k < j.(2.38)The time loop, Eqs. 2.29 to 2.35, continues until the integration time is fulfilled. The timeadvances by dt when δ˜b6 converges.33A considerable advantage of IAS15 is the adaptive timestep of the algorithm. After y and y′are converged, the current timestep is evaluated and compared withdtneed = dtused(10−9b˜6)1/7where b˜6 =maxi |b6,i|maxi |y′′i |. (2.39)The subindex i runs in the same way as in Eq. 2.35. If dtneed < dtused, then the step is rejectedand restarted with dt = dtneed; otherwise, the step, with the estimated values of y and y′, areaccepted, the time advanced by dtused, and the timestep changed to dt → dtneed in the nextstep. Rein and Spiegel (2015) chose to update the timestep using Eq. 2.39 to make the stepindependent of the physical scale of the problem, and therefore the timescales are capturedproperly without the intervention of the user. The value 10−9 was determined to address thepossible error sources, mainly the error due to the truncated approximation (Eq. 2.29), whichbuilds-up as dt15 if the timestep is constant. In Figure 2.3, we compare the behavior of therelative energy after 100 orbits of a simulation of the outer Solar System as run by Rein andSpiegel’s IAS15. The error, and therefore, the relative energy loss in the Rein and Spiegelimplementation and in ours starts to build as dt15 when the timestep is approximately 500days, and overall, both implementations agree well. When adaptive time-stepping is turned on,most of the additional error sources are random and with values under machine precision. Theother strong source of error is the precision lost due to truncation in sums between large andsmall numbers. Compensated or Kahan summation (Kahan, 1965) is used to track and reducethe build up of random errors, which grow ∝ t1/2. A complementary graphical summary of thesteps in IAS15 is shown in Fig. A.2.Figure 2.3: Comparison of relative energy after 100 orbits between Rein and Spiegel’s IAS15 andour implementation. The timestep was constant, i.e., the adaptive stepping was deactivated.The fractional energy of IAS15 must behave as dt15 when the error term is dominated by theGauss-Radau quadrature algorithm instead of random round-off error. Our implementationagrees with the previously published implementation.342.2 Secular codeThe equations of motion for a gravitational system with more than two bodies (N > 2) arenon-integrable. To understand the dynamics of an N-body system, it is necessary to integratethe equations of motion numerically, but this does not elucidate the functional form of thegravitational interaction or highlight how the system evolution depends on different orbitalelements. Secular theory is an analytical approximation to the long-term evolution of an N-body system. The approximation takes into account mutual gravitational perturbations fromplanets or other sources on a given Keplerian orbit and then averages those perturbations overtime. These perturbations are modeled as a power series of either Cartesian coordinates or theorbital elements.We calculate the long-term, or secular, eigenstructure of prototype STIPs by using the Laplace-Lagrange solution as laid out by Murray and Dermott (2000)22. The code is available athttps://github.com/norabolig/resmap and will be described briefly in the following subsec-tions. The code first calculates the characteristic frequencies (eigenfrequencies) of the longitudesof ascending node and pericenter of the secondary masses (e.g., the planets) as described inSect. 2.2.1. Knowing these eigenfrequencies allows us to evaluate how the system would perturba test particle. In some cases, an individual, low-mass planet can be treated as a test particle,and if so, we can evaluate the secular forcing on its eccentricity and inclination. However, thismust be done with caution, as described below (Sect. 2.2.2). The radial position of the testparticle can be varied to explore a range of parameter space and it helps to identify possibleeccentricity and inclination resonances.2.2.1 Secular theoryA secondary body in a hierarchical gravitational system with N > 2 masses feels the gravita-tional effects of the central mass and the other N−2 secondary bodies. As an example, considera three-body system as shown in Fig. 2.4. The equation of motion for the secondaries i and jassuming that the mass i is closer to the star than j (i.e., ri < rj) is given byr¨i = −G(mc +mi) rir3i+ Gmj(rj − ri|rj − ri|3 −rjr3j), (2.40)r¨j = −G(mc +mj)rjr3j+ Gmi(ri − rj|ri − rj |3 −rir3i). (2.41)22The original development of the Laplace-Lagrange solution was developed by Brouwer and Clemence (1961),but we followed the notation of Murray and Dermott (2000).35Figure 2.4: Geometry used for the secular theory development. The central body (the mostmassive) is denoted as mc, and the two companions by mi and mj , where mi is closer to thecentral body than mj (i.e., ri < rj). The dotted arcs represent the orbits of the masses mi andmj . Finally, O is the origin of the reference frame.The right-hand side can be represented by the gradient of two scalar potentials:r¨i =∇i (Ui +Ri) =∇i[G(mc +mi)ri+ Gmj(1|rj − ri| −rj · rir3j)], (2.42)r¨j =∇j (Uj +Rj) =∇j[G(mc +mj)rj+ Gmi(1|ri − rj | −ri · rjr3i)], (2.43)where Uk = G(mc + mk)/rk is the central potential felt by mass mk (k = i, j). The potentialRk, given by the second and third terms in Eqs. 2.42 and 2.43, is the disturbing function, whichis due to the additional secondary masses in the system. This can be explicitly written asRi = µj(1|rj − ri| −rj · rir3j)for the inner secondary and (2.44)Rj = µi(1|ri − rj | −ri · rjr3i)for the outer secondary, (2.45)where µk = Gmk. For simplicity, the disturbing function is divided into direct and indirectterms. The first term in Eqs. 2.44 and 2.45 is the direct part, RD, while the indirect part,given by the second term in both Eqs. 2.44 and 2.45, is represented by RE and RI. RE is theindirect contribution due to an external perturber (in Eq. 2.44), and RI is due to an internalperturber (Eq. 2.45). Using the ratio of semi-major axes αij = ai/aj < 1 and the dot productri · rj = rirj cosψ, the disturbing function can be expressed asRi = µjaj(RD + αijRE) (2.46)36Rj = µiai(αijRD + RIαij). (2.47)Furthermore,RD ≡ aj|rj − ri| , (2.48)RE ≡ −(riai)(ajrj)2cosψ, (2.49)RI ≡ −(rjaj)(airi)2cosψ. (2.50)Here, ψ is the angle between ri and rj .The direct and indirect parts of the disturbing function can be represented in a series expansionof the orbital elements that has a general formRi = µj∑S(ai, aj , ei, ej , Ii, Ij) cosϕ, (2.51)where S is a function of the semi-major axes (a), eccentricities (e) and inclinations (I) ofmasses i and j. The argument, ϕ, is a linear combination of the mean (λ), apsidal ($) andnodal longitudes (Ω), such thatϕ = k1λj + k2λi + k3$j + k4$i + k5Ωj + k6Ωi. (2.52)The coefficients ki are integers that sum to zero. Whether the disturbing function representssecular and/or resonant interactions depends on the product ki with the corresponding longitude(λ, $, or Ω). The number of relevant terms can be reduced if some of the cosine terms can beignored. This is motivated by the averaging principle, which assumes that all of the unimportantterms of the disturbing function will have short periods and their effects average to zero overlong timescales. In a two-body system, the mean longitude of the secondary mass increaseslinearly with time because λ˙ = n, but the apsidal and nodal longitudes are constant because$˙ = Ω˙ = 0. In comparison, when there are two or more secondaries in the system, λk variesrapidly, but $k and Ωk slowly. This means that the mean longitude has a shorter period thanthe apsidal and nodal longitudes.The long-period behavior of the argument in Eq. 2.52 depends as well on the possible commen-surabilities between k1λi and k2λj . If k1λi ≈ k2λj , then k1ni − k2nj ≈ 0 and the timescale ofthe interaction is long. In this case, terms associated with k1 = q and k2 = p − q, with p andq being integers, must also be considered in both RD and RE,I. The value of p provides theorder of the resonant terms and the strength of the resonance is proportional to the eccentricityS ∝ e|p| (with e < 1). As a consequence, the first-order resonances are typically the strongest(Murray and Dermott, 2000).37To proceed, the disturbing function is often expanded in orbital elements to the desired orderand a solution is then found for the simplified equations. In the Laplace-Lagrange solution,only the terms that are second order in e and I and first order in mass are retained. As aconsequence of the averaging principle, the secular terms of the disturbing function are thosethat have long periods, and therefore, do not depend on the mean longitudes, λ. This in turnrequires k1 = k2 = 0 in Eq. 2.52. Murray and Dermott (2000) show that the literal expansion ofthe direct part of the disturbing function to second order in e and I, assuming that I is small,is:R(sec)D =18αij b(1)3/2(e2i + e2j)− 18αij b(1)3/2(I2i + I2j)− 14αij b(2)3/2 eiej cos($i −$j) +14αij b(1)3/2 IiIj cos(Ωi − Ωj). (2.53)Here, I represents the inclination and the subindex i denotes a particle, and b(k)s (αij) are theLaplace coefficients defined as12b(k)s (αij) =12pi∫ 2pi0cos kψ(1− 2αij cosψ + α2ij)s dψ. (2.54)Then, from Eq. 2.46 and 2.47, the secular disturbing function for the inner and outer planetsareRi = GmjajR(sec)D (2.55)Rj = Gmiaiαij R(sec)D =GmiajR(sec)D (2.56)given that all terms in RE and RI depend of the mean longitude (see Eqs. 6.110 and 6.111Murray and Dermott, 2000). Expression 2.53 is only valid whenever mean motion commensura-bilities are absent. Otherwise, additional terms should be considered in R(sec)D and the indirectpart should be included.Substituting Eq. 2.53 into Eqs. 2.55 and 2.56, and using G = n2i a3i /(mc +mi) = nj2aj3/(mc +mj), where n is the mean motion, we can represent both Ri and Rj asRi = nia2i[12Aii e2i +Aij eiej cos($i −$j)+12Bii I2i +Bij IiIj cos(Ωi − Ωj)], (2.57)for i = 1, 2 and j = 2, 1, andAii = +14mjmc +mini αij α¯ij b(1)3/2(αij) (2.58)38Aij = −14mjmc +mini αij α¯ij b(2)3/2(αij) (2.59)Bii = −14mjmc +mini αij α¯ij b(1)3/2(αij) (2.60)Bij = +14mjmc +mini αij α¯ij b(1)3/2(αij) (2.61)whereα¯ij =αij external perturber, if i < j1 internal perturber, if i > j. (2.62)In this way, Eq. 2.57 can be represented in matrix form after recognizing that the terms Aijand Bij are the elements of the matricesA =(A11 A12A21 A22)and B =(B11 B12B21 B22).The Laplace-Lagrange secular solution makes use of the eccentricity and inclination vectors,as shown in Fig. 2.5, and defined ashi = ei sin$i, ki = ei cos$i, (2.63)pi = Ii sin Ωi, qi = Ii cos Ωi, (2.64)to substitute the corresponding terms in Eq. 2.57. Using trigonometric identities, the resultcan be expressed asRi = nia2i[12Aii (h2i + k2i ) +Aij (hihj + kikj)+12Bii (p2i + q2i ) +Bij (pipj + qiqj)]. (2.65)Now, we seek to use Eq. 2.65 to find the time evolution of the system, namely h˙i, k˙i, p˙i, andFigure 2.5: Graphical depiction of the eccentricity (left) and inclination (right) vectors describedin Eqs. 2.63 and 2.64, and used throughout Sect. 2.2.1.39q˙i. The corresponding relationships between the disturbing function and the derivatives of theinclination and eccentricity vectors are (Murray and Dermott, 2000, p. 278)h˙i =1nia2i∂Ri∂ki, k˙i = − 1nia2i∂Ri∂ki, (2.66)p˙i =1nia2i∂Ri∂qi, q˙i = − 1nia2i∂Ri∂pi. (2.67)For two secondaries (i = 1, 2), this results in a system of equationsh˙1 = A11k1 +A12k2, h˙2 = A21k1 +A22k2, (2.68)p˙1 = B11q1 +B12q2, p˙2 = B21q1 +B22q2, (2.69)for which the corresponding solutions have the form, recalling the expressions for the eccentricityand inclination vectors (Eqs. 2.63 and 2.64):hi = ei sin(g t+ β), ki = ei cos(g t+ β), (2.70)pi = Ii sin(f t+ γ), qi = Ii cos(f t+ γ). (2.71)By substituting Eq. 2.70 in Eq. 2.68, the problem is reduced to finding the eigenvalues of thematrix A and B, denoted by g and f , respectively. For example, if the system of equations inEq. 2.68 is expressed in matrix form after taking the derivative of hi (Eq. 2.70), then(h˙1h˙2)= g(e1e2)=(A11 A12A21 A22)(e1e2),which is an eigenvalue problem. The same is applicable for B in terms of p and q (Eq. 2.71). Thenumber of eigenvalues depends on the number of secondaries (N) considered in the calculation.Thus, the coefficients ei and Ii are the eigenvectors23 of the matrices A and B, whose valuesneed to be determined together with the phases βi and γi. The eigenvectors can be foundby replacing each eigenvalue in the corresponding set of equations and solving for the vector,but the resulting eigenvector is unscaled. The boundary conditions are necessary to scale theeigenvectors and to calculate the phases (see Sect 7.3 of Murray and Dermott, 2000, for moredetails).Once the eigenvalues and the eigenvectors are known, the time dependence of the componentsof the eccentricity and inclination vector are expressed ashi =N∑leil sin(gl t+ βl), ki =N∑leil cos(gl t+ βl),23Here, ei → eil and Ii → Iil, where l denotes the mode number associated with the eigenvalue gl or fl.40pi =N∑lIil sin(fl t+ γl), qi =N∑lIil cos(fl t+ γl).The eccentricity and inclination value at any time is given by the magnitude of the analogousvector.Physically, the eigenvalues found in the Laplace-Lagrange secular solution represent the the-oretical characteristic precession rates of the longitudes of pericenter and ascending node ofthe system. Throughout this work, we will refer to the eigenvalues as eigenfrequencies of agiven system due to their association with the precession rates. The gi eigenfrequencies areassociated with the precession rate of $i, while fi are associated with Ωi, where i indicates themode number and not a mass. The comparison of the theoretical secular eigenfrequencies andsynthetic precession frequencies (described in detail Sect. 2.3) provide insight into the dynam-ics of a system, particularly if the number or values of the fundamental frequencies (or both)mismatch. The discrepancy can be due to neglected resonant interactions in the theoreticalcalculation. For a case study presented in Ch. 4, we compare the eigenfrequencies from seculartheory with the numerically-derived principal frequencies for two different systems that havehigh planet multiplicity. The same principles laid out here can be applied to a system withmore than two secondaries, in which Eqs. 2.57, 2.65, 2.68, 2.69, 2.70 and 2.71 change to thecorresponding equations in Appendix A.2.2.2.2 Forced eccentricity and inclinationAfter the eigenfrequencies have been estimated, a test particle can be included in the seculartheory to map the gravitational influence of the secondary masses on the particle’s location.The number of eigenfrequencies in this new system increases by one without modifying thevalue of the eigenfrequencies found in Sect. 2.2.1 because the extra body has zero mass. Oursecular code exploits this feature to map the system’s secular behavior by placing many testparticles (one a time) over a range of semi-major axes. Thus, for a test particle with orbitalelements n, a, e, I, $ and Ω in a system with N secondaries, the particle’s secular disturbingfunction (analogous to Eq. 2.65) isR = na2[12A (h2 + k2) +12B (p2 + q2)+N∑j=1Aj (hhj + k kj) +N∑j=1Bj (p pj + q qj)]. (2.72)The solutions to this equation have the form of a forced harmonic oscillator:h = efree sin(A t+ β) + h0(t), k = efree cos(A t+ β) + k0(t), (2.73)41p = Ifree sin(B t+ γ) + p0(t), q = Ifree cos(B t+ γ) + q0(t), (2.74)whereh0(t) = −N∑iN∑lAieilA− gl sin(gl t+ βl),k0(t) = −N∑iN∑lAieilA− gl cos(gl t+ βl), (2.75)p0(t) = −N∑iN∑lBiIilB − fl sin(fl t+ γl),q0(t) = −N∑iN∑lBiIilB − fl cos(fl t+ γl), (2.76)andA =14nmcN∑j=1mj αj α¯j b(1)3/2(αj) (2.77)B = −14nmcN∑j=1mj αj α¯j b(1)3/2(αj). (2.78)In Eqs. 2.73 and 2.74, the first terms are the free (proper) eccentricity and inclination vectorsof the test particle, while the second terms are the forced eccentricity and inclination vectors.The forced eccentricity and inclination values are obtained from the normalization of vectors(k0, h0) (Eq. 2.75) and (p0, q0) (Eq. 2.76), respectively.Secular resonances originate when the denominators A−gl and B−fl in Eqs. 2.75 and 2.76 areclose to zero. This occurs when the particle’s semi-major axis generates an A or B value thatmatches one of the system’s eigenfrequencies. This situation, for either a secular eccentricityor inclination resonance, is of particular interest for the STIPs case study.2.3 Synthetic secular theoryThe synthetic secular theory comes from analyzing the frequency components of the orbitalprecession rates among planets in a given system. A system is evolved in an N-body integratorfor a set period of time. Then, the osculating orbital elements are calculated and transformed toFourier space to obtain the principal frequency components of the eccentricity and inclinationvectorse(t) = (h, k)(t) = (e sin$, e cos$) (2.79)42I(t) = (p, q)(t) = (sin I sin Ω, sin I cos Ω), (2.80)where e is the eccentricity, I the inclination, $ the longitude of pericenter and Ω the longitudeof ascending node. Finally, the value of the principal frequencies of the synthetic secular theoryare compared to the eigenfrequencies of secular theory. We emphasize that both synthetic andtheoretical secular frequencies represent the precession of the eccentricity or inclination vector,$˙ or Ω˙, respectively. When a system with N masses is characterized well by the Laplace-Lagrange theory, the synthetic and theoretical precession rates should match; otherwise, thedisagreement indicates that strong-gravitational or long-term interactions were neglected in thesecular theory calculation.The Fourier transform (FT) converts a time series into the frequency domain. The transfor-mation is implemented by using either the orthogonality of sin(2piνt) and cos(2piνt), or byusing exp(−j2piνt) functions. For the same reason, different normalizations can be found in theliterature. The FT used throughout this work isF (ν) = Ft[f(t)](ν) =∞∫−∞f(t) e−j2piνt dt, (2.81)where ν is the linear frequency, f(t) is a time-dependent function and F (ν) is the FT of f(t),which is a function only on ν. Because we will examine discrete times and frequencies, thecorresponding discrete Fourier transform (DFT) is used:Fn = Fk[{fk}N−1k=0 ](n) ≡N−1∑k=0fk e−j2pink/N . (2.82)The kernels in Eqs. 2.81 and 2.82 are related through the k-time t → tk = k∆t, and the n-frequency ν → νn = nN∆t , where N is the total number of time samplings and ∆t is the timesampling interval. In Eq. 2.82, fk is the discrete time series to be transformed to Fourier space,which is e(tk) and I(tk) in our case. We use the Python predefined function np.fft.rfft toperform the DFT, which implements a fast Fourier transform (FFT) for real fk data.The outcome of the DFT is the complex precession spectrum of e(jνn) and I(jνn). The am-plitude and phase of the spectra are of interest because the former indicates which frequenciesdominate the data and the latter shows how the sinusoids are aligned. The amplitude is givenbyAn =|Fn|/N for n = 0,2|Fn|/N for n > 0, (2.83)where |Fn| is the modulus of the complex spectrum, which is the magnitude of the correspondingeccentricity and inclination vectors as functions of the precession frequency νn. The frequenciesthat correspond to the highest amplitude values are the principal components of the spectrum.43The phase is given by the complex argumentφ = tan−1(Im(Fn)Re(Fn)). (2.84)Here, Im() and Re() denote the imaginary and real parts of the complex number, respectively.In the following two chapters, we will use the methods and codes described throughout thischapter. We point out that Mercury6 was just slightly modified to include the GR correctionusing Eq. 2.1. On the other hand, the IAS15 implementation was completely written by mefollowing Rein and Spiegel (2015) algorithm. In IAS15, we included, besides the planet-starand planet-planet gravitational forces, the GR correction used in Mercury6, and a prescriptionfor gaseous tidal damping of eccentricity and inclination (see Sect. 2.1.1, particularly Eq. 2.2).44Chapter 3Close-in giant planets as an extremeoutcome of planet-planet scattering.Volk and Gladman (2015) showed that the observed extrasolar planetary systems are metastable,i.e., in statistical realizations with slightly perturbed initial conditions, these systems tend to-ward instability in various timescales, especially those at short orbital periods. As a conse-quence, the multiplicity of STIPs tends to decay with time due to planet-planet scattering,collisions and consolidation. Thus, the high-multiplicity STIPs observed today are the longest-lived variants and could be representative of typical formation configurations. We suggest thatthe intrinsic metastability of STIPs leads to a variety of planetary types that form in situ atshort orbital periods. The properties of the planets produced in this scenario depend of thedetails of the consolidation, including the relative timing of instability compared with gaseouslifetimes. The formation of a HJ could be an extreme dynamical outcome of a STIP, should alarge core form through consolidation and enable runaway gas accretion. Here, we consider agas giant to be any planet that has at least half of its mass in captured nebular gas.The critical mass needed to initiate runaway accretion depends on different parameters, e.g.,the distance from the star, gas density profile, planetesimal accretion rate, etc. (Lee et al.,2014; Piso et al., 2015). In this study, we use a simple parametrization and define a criticalmass to be mc = 10M⊕, independent of the distance from the star. We also assume that theformation of short orbital period planets, with masses between 1 and 4 M⊕, occurs within aMyr of the disk lifetime (Dawson et al., 2015) and therefore, the gaseous disk is still present.Observations of circumstellar disks around young stars in clusters show that the fraction ofstars with detectable disks decays exponentially with time, ∝ exp(−t/λ), with an averagetimescale λ = 2.5 Myr (Mamajek, 2009). If we consider this to represent, on average, the decaytimescale for gas in a disk, then the disk’s density would decrease to 50% of the initial value byt = 1.7 Myr. Hence, any considerable accretion of gas by a planetary embryo should take placewithin the first 2 Myr. While the zero point for the time is ill-defined, we take t to be the timefrom STIP formation, which is envisaged to occur essentially along with disk formation.Fig. 3.1 shows different planet formation outcomes that could occur in the in situ formationframework. The successful in situ formation of a gas short-period giant (gSPG), with M >10M⊕ and 50% of its mass in gas requires: (1) the ubiquitous formation of a high-multiplicity45UbiquitousSTIP formationSTIP formed in GAS-RICH environment,NO critical-mass planet.One or more gSPGSTIP with a gas-poor SPGSTIP formed in GAS-FREE environment,NO critical-mass planet.YES Instability and planetconsolidation?Critical core viamerger BEFORE gasdisperses?Critical core viamerger AFTER gasdisperses?STIP forms BEFOREgas disperses?Critical core viamerger AFTER gasdisperses?NOYES YESYESYESNO NONONOFigure 3.1: In situ formation pathways of short-orbital period planets (adapted from Boleyet al., 2016). Timing between consolidation and gas dissipation/availability is key in thisparadigm, with the formation of a gSPG representing an extreme outcome. Dashed-lines rep-resent evolution pathways that we conceive uncommon.STIP within 1 Myr, and thus, before the gas disperses from the disk; (2) instability of theprimordial STIP in the presence of a gaseous disk, which leads to the collision and consolidationof planetary cores; (3) the subsequent consolidation of a critical core to start runaway accretionof gas, and; (4) a consolidation occurring within the first 2 Myr of the disk lifetime.Among the confirmed exoplanets, there is a large population with sizes between 1R⊕ and 4R⊕(Howard et al., 2010; Batalha et al., 2013; Petigura et al., 2013; Dong and Zhu, 2013; Fressinet al., 2013; Rowe et al., 2014), also known as super-Earths, with measured masses within2 ≤ Mp/M⊕ ≤ 20 (Wu and Lithwick, 2013; Weiss and Marcy, 2014). Close-in super-Earthsare puzzling due to their upper mass range and current separation from the host star, forwhich the upper masses could have underwent runaway accretion. Nonetheless, the estimatedgas content by mass of these objects is . 10% (Lopez and Fortney, 2014; Rogers and Seager,2010b,a), which suggests forestalled runaway gas accretion. In the paradigm suggested in thiswork, massive super-Earths originate from a different pathway than gSPGs. Super-Earths arecritical cores consolidated after 2 Myr, when the gas has significantly dissipated. Dependingon the accreted gas mass, the planet resembles either a mini-Neptune or a super-Earth. UnlikegSPGs, short-period giants (SPGs) are planets at or above the critical mass but that failed toaccrete a significant gaseous envelope. Examples of such a planet include K2-66 b (Sinukoffet al., 2017) and Kepler-62 d (Borucki et al., 2013). As shown in Fig. 3.1, other pathways arealso possible, e.g., rocky planets can form through accretion of solid material after the gaseousdisk disperses.Using numerical simulations, we determine whether it is plausible for critical cores to consolidatefrom metastable STIPs, and if so, whether such consolidation occurs at times during which gas46may still be present.3.1 Simulations and prototype STIPs.We evaluate the outcomes of STIP metastability and planet consolidation by running two sets ofN-body simulations. The first set of simulations includes, for simplicity, only planet-planet andplanet-star gravitational interactions, i.e., we neglect interactions between the planets and thegaseous disk. This part of the study serves as a proof of concept, exploring whether metastabilitycan drive the consolidation of critical cores within a few million years (Sect. 3.1.1).For the second part of the study (Sect. 3.1.2), we explore how planet-disk interactions can alterthe dynamics of STIPs. While this could be studied directly using hydrodynamics simulations,such simulations are very expensive, even for short integration timescales. An alternative is toinclude prescriptions for planet-disk interactions in an N-body integrator. This is the approachwe take, using a gaseous tidal damping algorithm based on the results of hydrodynamics simu-lations with embedded planets (Cresswell et al., 2007; Tanaka and Ward, 2004). Gaseous tidaldamping is an outcome of the interaction between the spiral density waves within a disk andthe planets that excite those waves. These spiral waves, in turn, reduce the orbital eccentricityand inclination of planets within a disk and can drastically affect planet orbital stability.Due to the different physics required for each of these parts, we use different N-body integratorsfor each study. In the simulations without gas drag, we ran 1000 realizations of K11 (see Table3.1) in Mercury6. For the simulations with gas damping, we used our implementation of IAS15to include non-conservative forces. We also generated synthetic six-planet STIPs, as we foundPlanet P Rp Mp a e I ω(days) (R⊕) (M⊕) (au) (◦) (◦)b 10.3039+0.0006−0.0010 1.80+0.03−0.05 1.9+1.4−1.0 0.09 0.045+0.068−0.042 0.12 45.0c 13.0241+0.0013−0.0008 2.87+0.05−0.06 2.9+2.9−1.6 0.11 0.026+0.063−0.013 0.07 51.3d 22.6845+0.0009−0.0009 3.12+0.06−0.07 7.3+0.8−1.5 0.15 0.004+0.007−0.002 0.15 146.3e 31.9996+0.0008−0.0012 4.19+0.07−0.09 8.0+1.5−2.1 0.19 0.012+0.006−0.006 0.63 90.0f 46.6888+0.0027−0.0032 2.49+0.04−0.07 2.0+0.8−0.9 0.25 0.013+0.011−0.009 0.05 90.0g 118.3807+0.0010−0.0006 3.33+0.06−0.08 < 25.0 0.47 < 0.15 0.35 90.0M∗ = 0.961± 0.025M R∗ = 1.065+0.017−0.022R Table 3.1: Nominal orbital elements of the known planets of Kepler-11 (K11) (Lissauer et al.,2011a, 2013). These values are a product of observational measurements and dynamical mod-eling. The nominal values provided by (Lissauer et al., 2013) for K11 are within 1σ confidenceintervals, except for the mass and eccentricity of K11g with 2σ upper bounds.47the K11 configuration to show no signs of instability in the presence of gas. All the simulationsuites with more than 300 members were run on the Orcinus computing cluster, provided byWestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca).3.1.1 Gas-free gravitational interactionsThe masses and spatial distributions of STIPs shortly after formation is unknown. However,due to their metastability, we argue that the high multiplicity STIPs observed today are amongthe longest-lived configurations of “initial” conditions after the main phases of planet building.Following this line of thought, we concentrated on using a well-characterized STIP, with morethan five detected planets, as a template for our first study.The number of systems detected to date with more than five planets is small. Nine of the2816 confirmed systems harbor more than five planets, and only Kepler-90 (K90) (with eightplanets now known) resembles the multiplicity of the Solar System (Schneider et al., 2011). Atthe start of this work, before the detection of TRAPPIST-1 (Gillon et al., 2016), the STIPswith the highest multiplicities included Kepler-11 (K11) (Lissauer et al., 2011a, 2013) andKepler-90 (K90) (Cabrera et al., 2014; Lissauer et al., 2014b; Schmitt et al., 2014), with six andseven planets24, respectively. K11 was chosen to serve as a prototype because the planetarymasses and semi-major axes are well characterized. In addition, constraints on the orbitaleccentricities, inclinations and arguments of pericenter (Table 3.1) were also available due toa combination of transit observations, TTVs, and dynamical modeling (Lissauer et al., 2013).The exception is K11g, the outermost planet, which only has an upper mass and eccentricitylimit of M < 25.0M⊕ and e < 0.15 provided by stability tests (Lissauer et al., 2013). Detailsof the K11 system are given in Table 3.1.We performed a suite of 1000 numerical simulations based on the nominal masses and semi-major axes of the K11 planets. We set the mass of K11g to be 8M⊕, which is well underthe upper limit constrained by stability analyses and reduces the possibility that this planetdominates the dynamics of the STIP. Thus, the total planetary mass of the K11 analogues isMT = 30.1M⊕. The orbital eccentricity, e, longitude of the node, Ω, argument of pericenter,ω, and mean anomaly, M, are uniformly randomized within the ranges shown in Table 3.2.The inclination, I, is also randomized but with a Rayleigh distribution with mode σI = 1.8◦,consistent with the observed inclination distribution of exoplanets (Fabrycky et al., 2014). Eachrealization is integrated for 20 Myr with a timestep of 0.1 days. Collisions are allowed to occurin the simulations.24In December 2017, an eighth planet was discovered in Kepler-90, Kepler-90 i, using deep machine learning(Shallue and Vanderburg, 2018)48Parameter Distribution Value range σI Rayleigh 0◦ − 6◦ 1.8◦ω, Ω, M Uniform 0◦ − 360◦ · · ·e Uniform 0− 0.05 · · ·Table 3.2: Summary of randomized orbital elements of K11 analogues. Here, σ is the charac-teristic parameter of the corresponding probability distribution.3.1.2 Planet-disk interaction: gaseous tidal dampingTo explore the effects of eccentricity and inclination damping on the stability of STIPs, weproduce 1000 random six-planet configurations. The systems were designed to be compact,more so than the observed typical mutual separation of neighboring exoplanets, i.e., 18−30RmH(Lissauer et al., 2014b), so that instability could be driven even in the gaseous disk, as completecircularization of planets in a system greatly increases the timescales for metastability. Forexample, many of the known multiplanet systems can exhibit instability on short timescaleswhen their orbital elements are given small perturbations, but are stable on very long timescalesif all the planets are given zero eccentricities (Fabrycky et al., 2014). The compact configurationsused in this study seek to explore potential initial states of planetary systems during theirformation in the presence of gas.The masses, eccentricities, and inclinations for the planets in the synthetic systems are randomlydrawn form a Rayleigh distribution with mode σm = 2M⊕ for the mass, σI = 1.8◦ for theinclination and σe = 0.025 for the eccentricity (as shown in Table 3.3). The initial planetarymasses are constrained to be in the range 0.1 ≤ Mp/M⊕ ≤ 8.0; if a mass is either smaller orbigger than this range, the generating script redraws the mass until all six planets have beenassigned a value. Fig. 3.2, shows the cumulative distribution of the resulting total planetarymasses for the systems. If instabilities occur during the dynamical evolution of these systemsin the disk, 90% of the systems have enough mass to form one critical core under the most idealsituations. However, in practice a much smaller fraction of systems could actually form criticalcores through collisions because typically only two to a few planets are consolidated. Using moreParameter Distribution Value range σMp Rayleigh 0.1− 8.0M⊕ 2.0 M⊕I Rayleigh 0◦ − 6◦ 1.8◦ω, Ω, M Uniform 0◦ − 360◦ · · ·e Rayleigh 0− 0.1 0.025M∗ = 1M R∗ = 1R Table 3.3: Summary of orbital elements of synthetic STIPs. Here, σ is the characteristicparameter of the corresponding probability distribution.49Figure 3.2: Cumulative total mass distribution in synthetic systems. From the synthetic STIPsgenerated, 90% have a total mass of > 10M⊕. The average total mass is 15±3M⊕, and a smallfraction of systems (∼ 10%) having less than 10M⊕ or more than 20 M⊕ in total. From thisplot, at least 90% of the systems have enough mass to form one critical core under the mostideal situations.massive systems is an area of future study. It should be noted that the K11 realizations havetwice the mass of a typical synthetic system (Fig. 3.2). This was due in part because we wereoverly cautious about setting up near-critical cores in the initial conditions. As a consequence,we unfortunately cannot directly compare the critical core formation rate in K11 with thispart of the study, which will instead be the focus of future work. The initial eccentricity andinclination distributions are based on results from the Kepler sample (Fabrycky et al., 2014).Once again, ω, Ω and M are uniformly randomized between 0◦ to 360◦. The mass of the staris set to M∗ = 1M .Following the discussion in Sect. 2.1.2, the planet densities are estimated from their mass usingEq. 2.11, which uses the Valencia et al. (2006) power law (ν = 0.272).The compactness of the synthetic systems is determined by the free parameter f = ∆a/RmH,which is the number of mutual Hill radii between the planets. The mutual Hill radius of aplanet j and its inner neighbor, j − 1, is given byR(j)mH =aj−1 + aj2(mj−1 +mj3M∗)1/3=∆aj + 2aj−12Kj , (3.1)whereKj ≡(mj +mj−13M∗)1/3. (3.2)Here, M∗ and m are the stellar and planetary masses, respectively, and ∆aj = aj − aj−1 is the50difference of planetary semi-major axes. The semi-major axes are iteratively calculated fromthe allocated masses and a constant value of f , using the expressionaj = aj−1(2 + f Kj2− f Kj)for j > 1. (3.3)The innermost planet semi-major axis is always set to a0 = 0.09 au.In addition to the initial conditions of the planets, the gaseous disk must be parametrized.In Sect. 2.1.1, we described how the gaseous tidal damping is included in IAS15 and thatthe strength of the damping depends of the local damping timescale (Eq. 2.5). We use thatprescription here, such thatτlocal(ρ, z, λ,Mp; t) = 0.078 days(MpM∗)−1 ( ρ1 au)1.9ez2/2H2 et/λ, (3.4)which is evaluated at each timestep, t, for a planet-to-star mass ratio Mp/M∗, and verticaland radial components, z and ρ. The only free parameter in Eq. 3.4 is the surface densitydecay timescale λ (recall the disk density and temperature profiles are fixed power laws for thisstudy). In preliminary simulations, we used a λ = 3 Myr and a mutual Hill radius parameterf = 9.5. For these parameters, the gaseous disk can still affect the dynamics after 20 Myr.The corresponding simulations required ∼ 30 days of computing time, in part because thetidal damping decreased the semi-major axis of the innermost planet, and in response, IAS15shortened the timestep to resolve the new orbit. To reduce the integration wall-clock time, weinstead used λ = 1 Myr. With this timescale, the system is practically gas-free at 10 Myr.Qualitatively, the change of decay timescale does not affect the dynamics of the STIP. Thepreliminary simulations also showed that with a mutual Hill separation f = 9.5, the numberof collisions was small. Therefore, we used f = 5 in the final suite of simulations to explorestability under more compact conditions. For reference, the mutual Hill spacings f = 5 and 9.5are, respectively, about 1.5 and 3.0 times the separation threshold for the stability of planetpairs (Sect. 1.1.4), i.e., f ≈ 3.46 (Gladman, 1993).3.2 Gas-free outcomes.In contrast to the nominal K11 system, which is stable for billions of years, the 1000 realizationsbecome unstable after short timescales. Within 20 Myr, the resulting fraction of unstablerealizations was 66.2%, producing consolidated masses between 5 and 17 M⊕. Most of themergers involved the two innermost planets, K11b and K11c, consolidating into a mass of5 M⊕.51Figure 3.3: Fraction of critical cores per interval using K11 as a prototype STIP. The differentcolors indicate the number of critical cores produced during that time interval. In black, weshow a proxy for gSPG formation after critical core consolidation. This is approximated byweighting each blue bin by exp(−t/2.5 Myr), which is intended to capture the decay of thegaseous disk. Adapted from Boley et al. (2016).In Fig. 3.3, we show the fraction of realizations that produced critical cores per time interval.In total, 20% of the realizations resulted in the consolidation of at least one critical core. Thedifferent histograms represent the number of critical cores in the indicated time interval. Twocritical cores form from consolidation in some cases, but we concentrate here on the fractionof systems that produce at least one critical core (blue histogram in Fig. 3.3). The highestfraction per interval corresponds to the first million years and it decreases with increasing time.In order to estimate whether the formation of a gSPG is possible after critical core consolidation,we introduce a gas disk mass decay timescale. With this in mind, we weight the core productionper time interval with an exponential decayexp(− t¯2.5 Myr),where t¯ is the bin’s midpoint. The proxy for gSPG formation is shown in black in Fig. 3.3.Cumulatively, 6.3% of gSPGs are generated within 10 Myr, with a minimal increase from 10 to20 Myr.3.2.1 Before moving through clouds: comparing the outcomes of Mercury6and IAS15The next stage of this study includes the implementation of gaseous tidal damping in IAS15.Before moving forward, we confirmed that Mercury6 and IAS15 give consistent results when52Figure 3.4: Comparison of critical core formation through consolidation between Mercury6(dashed lines) and IAS15 (solid lines).the simulations with identical initial conditions are run in the different codes without tidaldamping. To this end, we ran a suite of 1000 simulations in Mercury6 and IAS15. The initialconditions are similar to those described in Sect. 3.1.1, but with an eccentricity distributiondrawn from a Rayleigh distribution with σe = 0.025.Each simulation was run for 10 Myr. As before, the fraction of unstable K11 analogues is high.Using Mercury6, we found 80.2% of the systems to be unstable during the integration period,compared with 79.8% using IAS15. These unstable systems generated a range of planetarymasses through consolidation. Most of these unstable realizations originated from the consol-idation of planetary masses Mp < 10M⊕. Furthermore, of the total population, 22.4% and23.5% produced at least one critical core in Mercury6 and IAS15, respectively. The results ofthis comparison are shown in Fig. 3.4 and Table 3.4.The dynamical outcomes of both N-body integrators are consistent. The differences could bedue to numerical chaos or to the timestep management and collision detection algorithms. Forexample, identical simulations can diverge if the timestep, which would normally sample theorbit of a single planet well, is too large to properly sample the dynamical interactions betweenmultiple planets (Dones et al., 2018; Lissauer, 2018). We note that the IAS15 integrator iteratesMercury6 IAS15(%) (%)Unstable systems 80.2 79.8At least one critical core produced 22.4 23.5Two critical cores produced 4.8 5.3Table 3.4: Summary of results comparing outcomes of Mercury6 and IAS15.53on the solution for each step, and may better capture some of these dynamics. Regardless ofthe small differences, we found there to be reasonable agreement between Mercury6 and IAS15.3.3 Outcomes with tidal dampingGaseous tidal damping stabilized our synthetic systems, even though they were compact. Theresulting unstable fraction is small and depends on the parameter values for the initial condi-tions. A summary is shown in Table 3.5. No critical cores were produced and the consolidatedmasses ranged from 0.64 to 7 M⊕, well below the critical core threshold. For the decay timescaleλ = 3.0 Myr and Hill mutual spacing f = 9.5, we only found 2 unstable realizations out of 465.A limited number (31) of realization were run with λ = 1 Myr and f = 9.5, of which none wentunstable within 10 Myr. Only a few tens of realizations were run because we noticed that thechange of λ does not modify the overall behavior of the systems; it only changes the evolutiontimescale. This effect is due to the tidal damping timescale being much shorter than the gasdissipation timescale, as planets are damping within a few thousand years, depending on theFigure 3.5: Example of the dynamical evolution of a synthetic STIP embedded in a gaseousdisk. The semi-major axis, a, pericenter, q, and apocenter, Q, vs. time are shown for eachplanet (delineated by color). Because the planets typically have very low eccentricities, the a,q and Q curves appear as one thick curve for most times. In general, the systems spread insemi-major axis, with the inner planets moving to lower a and the outer planets moving tohigher a. The jumps in semi-major axis correspond to the location of first-order MMRs. In thisplot, f = 5.0 and one of the planets is not shown for reasons that will be discussed in Sect. 3.4.54Figure 3.6: Distribution of all possible combinations of period ratios for systems with f = 5.0.The vertical dotted-lines indicate the location of MMRs. The initial (top) and final (bottom)distribution are displayed in black. Top: the colored filled histograms represent different initialcloseness ratios, i.e., the closest neighbor is denoted as Pk+1/Pk, while the third closest asPk+3/Pk, et cetera. The width of these peaks is due to the initial planetary mass distributionchosen. Bottom: the period ratio of the closest planets (in blue) at the end of simulation showthat the planets cluster between MMRs. This clustering is different from that seen in the initialconditions, the planets spread significantly from their initial configurations.actual gas density. However, the gas damping timescale increases with decreasing gas density,allowing planet-planet interactions to have stronger perturbations. This increases the rate thatthe planetary spacing diffuses due to excitation followed by gaseous damping. In the last set,with λ = 1 Myr and f = 5.0, the number of unstable systems was the highest with 45/1000.We attribute this to the reduced Hill separation.The general dynamical behavior of the synthetic systems is highlighted in Fig. 3.525. The plotshows the semi-major axis, a, pericenter, q, and apocenter, Q, for the planets in a syntheticsystem as a function of time. The planetary semi-major axes slowly spread out, as a consequenceof having their eccentricities damped. Jumps in evolution also occur over ∆t < 1 × 104 yrwhenever planets cross first-order MMRs.In Fig. 3.6 and 3.7, we compare the initial and final period ratio distributions of the synthetic25We highlight that, even though only five planets are shown, the system is stable. For reasons that willbecome evident in Sect. 3.4, we neglected a planet in this plot. Regardless, the figure illustrates the generalbehavior of the systems well.55Figure 3.7: Distribution of all possible combinations of period ratios for systems with f = 9.5.The same notation and colors used in Fig. 3.6 are used here. Both, initial and final period ratiodistributions are similar to those in Fig. 3.6, except the initial distribution extends to largerperiod ratios due to the value of f .systems. All possible unique planet combinations are considered. The planet period ratio iscolor coded by the proximity of the pairs. For example, adjacent pairs are denoted k+ 1, pairsseparated by one planet are denoted k+2, etc. The chosen Hill radius spacing causes the periodsto be grouped into distinct peaks, with widths ultimately reflecting the mass distribution. Forf = 5.0, the period ratio of the closest neighbor peaks near the 8:7 MMR in the initial conditions,while, for f = 9.5, this occurs near the 5:4 MMR. This initial period ratio distribution close toan MMR is a fortuitous result of the chosen initial mutual spacing and mass distribution. Byconstruction the initial period ratio distributions are different between systems with f = 5.0and f = 9.5, the latter being spread to higher values.The final period ratio distribution indicates that the planets avoid the first-order MMRs. Thisis due to an eccentricity kick followed by gaseous tidal damping whenever the planet is near oneof these commensurabilities. The combined effect is to cause a small change in the semi-majoraxis, and subsequently remove the planet from the MMR. The excursions into first-order MMRsare short-lived and the width of the excursion is wider for the 2:1 MMR and narrower for the8:7. This behavior is present in the final distributions of f = 5.0 and f = 9.5 systems. Thesimilarity between the final period distributions for the f = 5.0 and f = 9.5 simulations leads usto suggest that in presence of a gaseous disk, the initial planet distribution is unimportant. This56λ f t¯uns Nsim Nuns M > 10M⊕ Effective integration time(Myr) (RmH) (yrs) (Myr)3 9.5 1600, 30000 465 2 No 0.8-2.21 9.5 · · · 31 0 No 10.01 5.0 200 1000 45 No 2.37-6.21Table 3.5: Summary of gaseous tidal evolution results. The mean instability time, t¯uns, is on the3rd column. The 4th and 5th columns are the total number of run and unstable realizations,respectively, in a given set. Ideally, each suite would have the same number of realizations andintegrated the same number of e-folding time relative to λ. The first and third set were run inthe Orcinus computing cluster (provided by WestGrid and Compute Canada Calcul Canada).mechanism might be a byproduct of the assumptions made for the disk, although the combinedeffects of dynamical excitation followed by gaseous dissipation should be general. Interestingly,the 5:2 MMR appears to have a local maximum (although this is not a distinct peak).From Table 3.5, it can be seen that the number of realizations and the effective integrationtime is different between the simulations suites (e.g., those with different λs and fs). Ideally,each suite would have the same number of realizations and would be integrated for the samenumber of e-folding times (relative to λ). However, due to the computing resources involved,we focused on λ = 1 Myr and f = 5.0 suite because it evolved faster and captured the essentialbehavior.3.3.1 Stability test after complete gas dissipation.We finally explore the stability of the final configurations of the synthetic systems in the sim-ulations with gaseous tidal damping. First, we use the orbital elements at the end of eachsimulation for the f = 9.5, λ = 1 Myr suite. We used this suite for three reasons: (1) systemswith f = 5.0 and f = 9.5 have similar end states, with a two-peak mutual spacing distributionat f = 14 and 17.2, the latter being the highest peak (see Fig. 3.8); (2) all the realizations inthis suite were evolved for the same amount of time (10 Myr) and none became unstable; and(3) this suite has a more tractable number of simulations. We integrate the final configurationsfor 100 Myr using Mercury6. Due to the small eccentricities and inclinations, none of the re-alizations went unstable within 100 Myr. This means that the gas removal does not affect theeccentricity and inclination of the planets.Next, we repeated the experiment, but assigned each planet in each system a random eccen-tricity and inclination drawn from a Rayleigh distribution with mode σe = 0.01 and σI = 1.8◦,respectively. After 100 Myr, 12 of the 31 realizations became unstable, from which 50% of theremaining planets have an eccentricity e > 0.018, with a maximum of e = 0.1.The results of these tests seem to suggest that: (1) the final configurations are stable if both57Figure 3.8: Final mutual Hill spacing distribution, f , for systems with initial separations f = 5.0(top panel) and f = 9.5 (bottom panel). The black line shows all possible combinations ofmutual Hill spacings within a system, while the blue highlight shows only the separationsbetween adjacent planets. The end states of both sets of simulations peak at mutual Hillseparations of f = 14.0 and 17.2, particularly if only the nearest neighbor is considered. Theseparation between planets within a system, as shown in black, can span up to 100 mutual Hillspacing.inclination and eccentricity are close to zero; and (2) the instabilities can still occur if there is amechanism that can perturb e and i among the planets, either within the gaseous disk or afterits dispersal.We should point out that only gaseous damping is considered in this work. The physics of aplanet interacting with the surrounding disk is rich, and there are mechanisms that competewith the gaseous tidal damping that can excite the planetary eccentricity and inclination.Ben´ıtez-Llambay et al. (2015) and Eklund and Masset (2017) showed that as a planet heats updue to planetesimal accretion, it re-radiates the thermal energy and changes the local conditionsof the disk, which prevents the inward migration of the planet. This mechanism can also excitethe eccentricity and inclination of the planet.3.3.2 Very compact STIPs in the presence of gas.Two additional sets, with 200 realizations each, were run to further test the instability of verycompact STIPs in the presence of gas. We used the same disk properties, planetary mass58Figure 3.9: Final distribution of period ratios for systems with initial f = 1.5 (top panel) andf = 3.0 (bottom panel). Similarly to Fig. 3.8, the black line indicates all possible combinationsof period ratios within a system, while in blue only the nearest neighbors were considered.The period ratio distribution from these simulations also shows an avoidance of the first-orderMMRs (as in Fig. 3.6 and 3.7) and the main peak (between the 3:2 and 5:3 MMRs) is stillpresent but with a broader width towards the 2:1 MMR. The period ratios shorter than 1.5are not as populated as in the f = 5.0 and f = 9.5. Also, there is a populated region betweenthe 2:1 and 5:2 MMRs, absent in Fig. 3.6 and 3.7.distribution and integration time as that used in the f = 5.0 experiment, but with initialmutual spacings of f = 1.5 and 3.0, respectively. Due to the initial compactness of thesesystems, the instability and subsequent consolidations were more efficient than in the othermutual separation experiments. The instability occurred early in the evolution of the systems,usually within a couple of years and occasionally after a few thousand years. The fraction ofunstable systems in the f = 1.5 and f = 3.0 sets were 96% and 86%, respectively. Theseextremely unstable and compact systems, unlike the more widely spaced realizations, producedconsolidated masses Mp ≥ 10M⊕. In the f = 1.5 set, 13% of the realizations produced at leastone critical core, and even a smaller fraction of systems (2.5%) consolidated two critical cores.On the other hand, in the f = 3.0 set only 1.5% of the systems formed only one critical core.While the setup is itself ad hoc, it does allow us to further explore instabilities in the presenceof gas and how such instabilities might shape, for example, the period distribution.The final period ratio distribution of these sets have a similar distribution to Fig. 3.6 and3.7. As shown in Fig. 3.9, the remaining planets avoid the first-order MMRs and have a very59Figure 3.10: Final mutual Hill spacing distribution, f , for systems with initial separationsf = 1.5 (top panel) and f = 3.0 (bottom panel). The same notation and colors used in Fig. 3.8are used here. Compared to Fig. 3.8, the peak at f = 14 is not significant in the distribution.Instead, a feature at f = 25 has emerged almost as prominent as the peak at f = 17.2.significant peak between the 3:2 and 5:3 MMRs. The main difference is that due to the initialmutual spacing, which causes strong instabilities, the peaks between first-order MMRs shorterthan 3:2 are not as populated as in the f = 5 or f = 9.5. This indicates that the planetsthat consolidated were preferentially those with period ratios shorter than 1.5. The survivingplanets then evolved like the more relaxed configurations.The final mutual spacing distributions of both sets, shown in Fig. 3.10, look very similar amongthem and also have similarities with those correspondent to f = 5.0 and f = 9.5 sets. Theprincipal peak is still located at f = 17.2, but the peak at f = 14.0 has almost disappeared.Instead, a peak at f ≈ 25 has emerged in both f = 1.5 and f = 3.0 initial mutual spacingdistributions. This new peak is almost as significant as the f = 17.2 peak in the f = 1.5 finaldistribution. The f = 25 peak seems to indicate the high instability occurrence in these setsand the subsequent clearing of some of the closest neighboring planets in the early evolution ofthe systems. In many of the simulations, the planets that collide are those with intermediatesemi-major axes, which results in the survival of the planets of the more distant planets. As theplanets are tidally damped by the gas, their mutual spacing also increases. This new feature isa consequence of both instability and tidal damping.603.4 Coorbital planets and their detectionWhile we did not observe the formation of critical cores through the consolidation of planetsin the λ = 1 Myr and f = 5.0 simulation suite, we did find 23 realizations in which planetsevolved into crossing orbits but did not collide. In most of those cases, only one orbit wascrossed, but in three realizations, a planet crossed the orbit of another two planets. In oneparticular case, the second innermost became the outermost planet within 10 years (about 200orbits), crossing the orbits of four planets. Although such sudden instability at the start ofthe simulation reflects very unstable initial conditions, the resulting dynamics are of potentialinterest, as highlighted here.In two of the simulations, neighboring planets apparently became locked in a 1:1 resonanceduring orbit crossing. The evolution of the semi-major axes of the planets in two realizations,noted as r-054 and r-628, are shown in Fig. 3.11. This behavior suggest that either the leastmassive planet was captured as a satellite/binary planet or the planets did indeed enter a 1:1coorbital resonance.We looked at the offset angle between the planets near the 1:1 MMR (Fig. 3.12), as well as theirorbits in a frame corotating with the most massive planet of the pair (Fig. A.1)26, to determinethe nature of the interaction. The offset angle of a satellite orbiting a planet has a maximumangle that is equal to the ratio of the satellite-planet and the planet-star distances27 (' 0.15◦in Earth-Moon case). In contrast, coorbital planets are identified by their libration near theLagrange points, and can be either in a tadpole or horseshoe configuration. The offset anglebetween coorbital planets is associated with three of the equilibrium points of the restrictedthree body problem, i.e., φ = ±60◦ (marked in green in Fig. 3.12 and A.1) and φ = 180◦ forthe Lagrange points L4, L5 and L3, respectively. If in a tadpole orbit, the offset angle libratesclose to either L4 or L5. On the other hand, if a planet is in a horseshoe orbit, the planetencloses L4, L5 and L3, which means that the libration is centered at 180◦ and the amplitudeof libration ∆φ < 180◦.As shown in Fig. 3.12 (bottom panels), in both of the simulations, the smaller planet becomestrapped in the 1:1 MMR after about 75 years. The offset angle between the planets firstlibrates as a tadpole orbit, followed by changes in and out of a horseshoe orbit as the planetsmigrate inward (Fig. 3.11(a)) or outward (Fig. 3.11(b)). For r-054, this only last 3% of the totalintegration time, after which the 0.25 M⊕ planet breaks out of the resonance and migrates to asmaller orbital distance. In the case of r-628, the coorbital configuration changes from a tadpole26Fig. 3.12 succinctly shows the dynamics of the coorbital realizations. The face-on plot of the corotating framerepresents the spatial evolution of these realizations, and for completeness, the corresponding plot is included inAppendix A.3.27The offset angle in this case is calculated from tanφmax = φmax ≈ am/ap, where am is the semi-major axisof the satellite in the planet’s reference frame and ap is the planetary semi-major axis respect to the star. Here,we used the small-angle approximation.61(a) r-054(b) r-628Figure 3.11: Semi-major axes, a, pericenter, q, and apocenter, Q, as a function of time for tworealizations with planet orbit crossings. Each system is identified by its run number. In bothrealizations, there is a planet that within the first million years crosses the orbits of another twoplanets, and is then captured into a 1:1 MMR with either (a) the innermost or (b) outermostplanet.62(a) r-054 (b) r-628Figure 3.12: Offset angle, φ, between the planets near the 1:1 MMR (top) and a detailed radial evolution as a function of time (bottom).The first column represents r-054 and the second r-628. The planetary masses are displayed in the legend. The offset of the L4 (φ = 60◦)and L5 Lagrange (φ = 300◦) points are indicated with a green horizontal dotted-line. (a) In this realization the smaller planet is trappedin a coorbital configuration for about 1000 years, after which it breaks free from the resonance due to the gaseous tidal damping. Forthe system in panel (b), in contrast, the smaller planet is captured into the 1:1 MMR by the first hundred years, starting in a tadpoleorbit which then evolves to a horseshoe orbit until the end of the simulation.63to a horseshoe orbit at 1 Myr and remains in this configuration until the end of the integration.Because the integration with gas is limited in time, we further test the stability of r-628 bytaking the final orbital elements and integrating them forward for 10 Myr using Mercury6,assuming the gaseous disk has dissipated. We run two different versions of this final state. Inthe first one, all the orbital elements are used at face value, which implies co-planarity and nearcircular orbits. For the other, the eccentricity and inclination of each planet in the system areslightly perturbed using a Rayleigh distribution with a mode of 0.01 and 1.8◦, respectively.Only the excited system becomes unstable. Within the first 250 kyr, the coorbital planetsbreak the 1:1 MMR and the least massive planet collides with one of the inner planets. Theplanetary eccentricities for the unperturbed and perturbed simulations are shown in Table 3.6.The planets e and g, those in the 1:1 MMR, have an increased eccentricity of about 2 ordersof magnitude in the perturbed case. These results suggest that if a coorbital configurationforms from a dissipative interaction with the gaseous disk, it can be long-lived should otheradditional perturbations be absent. Long timescale stability studies for equal-mass planets inthe 1:1 MMR, with e1 = e2 = 0.01, found that the mass threshold for horseshoe configurationsis µ = 2mp/(mp + m∗) < 4 × 10−4 for long-lived pairs (Laughlin and Chambers, 2002). Thissuggests that the coorbital pair, with µ ≈ 5×10−6 in r-628, will be stable for longer than 10 Myr.The eccentricities of the coorbital planets in our perturbed case are a factor of 2 higher thanLaughlin and Chambers (2002) study, but they only considered the two coorbital planets intheir simulations. It is likely that the presence of the other four planets in our simulations havea big impact on the stability of the perturbed realization.Fig. 3.13 shows the semi-major axis of each coorbital planet as a function of time for 100 yearsof integration of simulation r-628. In the horseshoe orbit, energy is exchanged periodicallybetween the two planets, which causes the semi-major axes and the orbital periods to varyPlanet a eu ep(au)b 0.04732 7.3× 10−5 1.7× 10−2c 0.06713 7.6× 10−5 6.5× 10−4d 0.11410 8.9× 10−5 1.0× 10−2e 0.20910 2.2× 10−4 1.6× 10−2f 0.16577 1.5× 10−4 7.2× 10−3g 0.20825 2.4× 10−4 2.4× 10−2Table 3.6: Eccentricity of planets in r-628. Here, etext, represents the eccentricity values of theunperturbed system, where the initial conditions where taken directly from the final state ofthe tidally damped simulations. ep, on the other hand, are the values for the perturbed system,drawn from a Rayleigh distribution with σe = 0.01.64Figure 3.13: Semi-major axes of coorbital planets as a function of time. This realization wasrun for 100 years, and it shows that the semi-major axis of each planet varies periodically in awell-constrained range. An additional horseshoe period of 40.3 years is present (indicated witharrows).with time. This orbital period has two components: the azimuthal period28 and a radial, Pr,horseshoe period. In this case, the horseshoe period is Pr = 40.3 years, 421 times longer thanthe azimuthal component. The mean azimuthal period for the planets is P = 34.733 days,about which the more massive planet varies ∆Pl = ±0.028 days and the less massive planetwith ∆Ps = ±0.198 days. The functional form for the radial period of tadpole planets is known(Laughlin and Chambers, 2002), which is a function of the planets-to-star mass ratio and theKeplerian orbital period. That relation, however, does not strictly apply for horseshoe orbits.Our simulations suggest that the radial period of a horseshoe orbit is ≈ 2.5 times longer thanthat of a tadpole.The 1:1 corotational MMR can provide additional information about the planets, other thanjust the orbital periods. If the difference between the semi-major axes, s = a2−a1, is measuredat two different times, i, j, and the corresponding offset angle, φ, can be inferred, then theplanetary masses can be determined using the following relations (Murray and Dermott, 2000):(sia)2 − (sja)2= −43m1 +m2mc[H(φi)−H(φj)] , (3.5)m1W1 = m2W2, (3.6)andH(φ) ≡(sinφ2)−1− 2 cosφ− 2. (3.7)Here, a is the semi-major axis of the Keplerian orbit; mc, m1, and m2 are the masses of thestar and planets, respectively, such that 0 < m2 < m1 and mc  m1,m2; and W is the average28This is the period that either of the two planets would have on a circular orbit at a = 0.2083 au.65Figure 3.14: Normalized deviation from a Keplerian orbit, s/a = (al − as)/a, as a function ofoffset angle, φ, for r-628. The shape of this plot is typical for horseshoe orbits, which highlightsthe orbital evolution (100 years) of both coorbital planets and provides information such as s180and φmin (magenta and green lines, respectively) needed to infer orbital properties. The meansemi-major axis of the planets is a = 0.21 au, as indicated.length of the arc, i.e., the angular spread in the corotating mean frame of each planet’s orbitwhen observed for a prolonged time. Eq. 3.5 can be simplified to(s180a)2=43m1 +m2mc[H(φmin)− 1] , (3.8)where s180 is the difference in semi-major axes of the two planets when separated by φ = 180◦,and φmin is the minimum phase difference between the planets. For r-628, we show in Fig. 3.14the semi-major axis deviation from a Keplerian orbit for the two coorbital planets as a functionof their offset angle. Using Eqs. 3.8 and 3.7, and the values s180/a = (1.852 ± 0.133) × 10−3and H(φmin = 21.46◦) = 1.509, we estimate a total planetary mass m1 + m2 = (1.005 ±0.141)×1025 kg = 1.68±0.24M⊕, which agrees with the known total mass, i.e., 1.617M⊕. Theuncertainty in the total mass calculation comes from the uncertainty in the measurement ofs180 within 170◦ ≤ φ ≤ 190◦. Should coorbital planets exist and be detected, their densitiescould be estimated with a single observation method, in the same spirit as using TTVs (Agolet al., 2005).Similar coorbital configurations have been found in studies of planetary systems embedded indissipative gaseous disks (Beauge´ et al., 2007; Cresswell and Nelson, 2006, 2008). Currently,it is uncertain whether coorbital planets exist, but their detection can provide a constrainton primordial formation mechanisms. For this reason, several detection methods have beenproposed based on TTVs (Schwarz et al., 2015; Ford and Holman, 2007; Vokrouhlicky´ andNesvorny´, 2014; Haghighipour et al., 2013), transits (Janson, 2013; Hippke and Angerhausen,2015) and RV (Ford and Gaudi, 2006; Leleu et al., 2015). Using some of these techniques,66several archival surveys have been carried out (Lillo-Box et al., 2018; Hippke and Angerhausen,2015; Janson, 2013; Ford and Holman, 2007; Madhusudhan and Winn, 2009), with no coorbitalplanets detected.Upon reviewing detection methods for coorbital planets, we realized there is not an explicitrepresentation of a horseshoe transit lightcurve, only ones for tadpole orbits (Hippke and Anger-hausen, 2015). For this reason, and to further efforts to detect trojan planets, we tested thefeasibility of detecting the coorbital pair in a horseshoe configuration with the transit method.We produce observational transit lightcurves from N-body simulations and then put those curvesthrough the Kepler Transit Planet Search pipeline (TPS) to determine whether coorbital planetscould be detected with current transit detection methods.For producing the transit lightcurves, we ran a short N-body simulation that mimicked Kepler’slong cadence observations by integrating for 4 years and using an output of 30 minutes. Theinitial conditions were those used in the stability test of the coplanar, circular r-628 simula-tion (see above). Synthetic lightcurves were then produced from the position vectors of bothcoorbital planets by choosing an arbitrary line of sight aligned with the mean orbital planeand reducing the normalized flux each time one of the planets crossed the projected surfaceof the star. The radii of the planets, Rl = 1.11R⊕ and Rs = 0.6R⊕, were obtained from themass-radius relationship in Eq. 2.9 and ν = 0.272, while the stellar size was set to 1R . In thisapproach, stellar limb darkening was neglected while producing the transit shape. To make thelightcurve more realistic, we included a Gaussian random noise of σ = 5 × 10−6. One impor-tant variant to consider is the initial time of integration, which observationally, corresponds toobserving epochs, as a particular lightcurve will be produced for different starting times due tothe changing relative positions of the planets.We synthesized three different transit lightcurves (Fig. 3.15), each showing distinct behaviordepending on the relative location of the planets when the integration begins. On the leftside of Fig. 3.15, we show the orbital period of each planet (top) and the relative period asa function of time. The colored regions represent the different times during which lightcurvesare produced. In the blue region, the planets just had their closest approach and the relativeperiod has a change in direction, which yields a lightcurve with transits for the two planets thatfirst decreases in relative phase and then increases. We call this signature non-linear behavior.In orange, the relative period changes linearly, which is shown in the middle lightcurve. Onthe third lightcurve, highlighted in green on the left panels, the relative timing of the transitsseems to be constant with time, which occurs when the planets are 180◦ apart.67Figure 3.15: Transit lightcurves of a coorbital pair in a horseshoe orbit at two different epochs (left), and corresponding period changeas a function of time (right). The color highlights in the left panels identify each of the lightcurves on the right. Here, Pl and Ps arethe orbital period of the large and small planet, respectively. We show that there is a linear and non-linear change of periodicity of thesmaller planet, which correspond to the orange and blue highlighting on the left, respectively. The third lightcurve (green highlight)corresponds to an offset angle between the two planets of about 180◦.68The synthetic linear and non-linear transit lightcurves were then given to collaborator MichelleKunimoto to run through the TPS. In both cases, the pipeline easily recovered the largerplanet, in part because the flux drop Rp/R∗ = 0.01 is well above Kepler’s detection limit of30 ppm for solar type stars (Gilliland et al., 2011). For this planet, the pipeline recovered aperiod Pl = 34.8 days in both lightcurves. Using the standard pipeline, the smaller planet(Rp/R∗ = 5×10−3) was not recovered from the lightcurves. However, the planet’s period couldbe recovered upon forcing the pipeline to look for a signal that is within 10% of the largerplanet’s period. In this case, the TPS could retrieve the second planet, with an observed periodof Ps = 34.6 days, but only for the linear lightcurve. The smaller planet evaded detection inthe non-linear case.A common way to determine the transit profiles, amplitudes and period variations is to phase-fold the lightcurve. Folding a time series, like a lightcurve, involves identifying a periodic signaland the corresponding period, P0; then the phase is calculated by dividing the time series byP0, such as ψ = (t− t0)/P0, where t0 is the center of the first transit. At this point the phasevaries from 0 to a given integer number of the period. The folding is finalized when the phaseis normalized to a range −0.5 ≤ ψ ≤ 0.5. With this phase-folding, the first transit in the datais centered at a phase 0, and the following transits should align on top if P0 represents the realperiod of the data. In the case of small deviations from exact periodicity, the feature at phase 0looks smeared or shifted from the center. Although the period of both coorbital planets changeswith time, the relative period variation of the larger planet is shorter than that of the smallerplanet. Therefore the transit lightcurve can be folded using the period of the larger planet(Pl = 34.8 days). By phase-folding the lightcurves with Pl, the transit of the smaller planet willappear smeared, showing an arc (see Fig. 3.16). In the folded lightcurve, the larger planet is thereference and the phase difference behaves much like observing the corotating reference framein phase values rather than relative positions. A phase difference provides a direct estimationof the angular separation between the two planets, related by φ = 360◦ ψ. This is why, thateven if the second planet cannot be detected, the transit signatures of the larger planet in thefolded lightcurve could still provide additional information about the system. For example, thedispersion in the transit signature could constrain W and φ (Eqs. 3.6 and 3.7).Fig. 3.16 shows the folded lightcurve for the three epochs, corresponding to the non-linear,linear and maximum separation regimes. In the non-linear case, when the planets go from alarger to a smaller semi-major axis, or vice versa, the transit of the larger planet is broadenedand the extent of the arc of the smaller planet is relatively small. Roughly, from this lightcurve,we can estimate ml/ms ≈ 8.5 and φ ≈ 19◦ by measuring the lengths of the arcs in the phaseplot and by taking the difference of the left edge of the deeper transit to the right edge of theshallower dip. In the linear case, at phase of zero, there is no smearing and the length of thesecond transit is larger than in the previous case. The mass fraction could also be obtained fromthis lightcurve. On the last lightcurve, when the planets are separated by 180◦, the important69Figure 3.16: The non-linear, linear and maximum angular separation lightcurves from Fig. 3.15,but phase folded to the period of the larger planet. The smaller planet produces the shallowtransit arcs. The position of the arcs depend on the relative angular separation between thetwo planets in the corotating frame. The length of the arcs, and particularly the ratio betweenthem, provides information on the mass ratio of the two planets (Eq. 3.7). From the non-linearlightcurve (at the top), in principle, we can also obtain the minimum angular separation of theplanets. The locations of L4 and L5 for the given period are shown for reference.parameter to estimate would be s180/a. It may be possible to relate this parameter to theobserved period, but this connection has not yet been explored here.These preliminary tests show that the transit method is promising for detecting coorbital planetsand should be explored further in future work. This could include, for example, using systemsthat produce a more realistic signal-to-noise value, provided the synthetic systems are within thestability limit µ = (m1 +m2)/mc < 4× 10−4 (Laughlin and Chambers, 2002). In principle, wecould use planets 100 times more massive for the coorbital configuration than those exploredhere. Updating the present transit lightcurve pipeline to search within 10% of a just found70period will improve the likelihood of recovering a coorbiting planet. Furthermore, if two planetsin a coorbital configuration are large enough to produce a transit above the survey observationallimit, and these are coplanar as well, they will produce a transit lightcurve with characteristicssimilar to those in Fig. 3.15. The phase folded lightcurve could be used to find m2/m1 and µ,but a relationship between s180 and P is still needed.71Chapter 4STIPs and external planetaryperturbers: effects on observabilityand stability.Due to the selection effects of the different planet detection techniques, it is uncertain whetherouter companions exist in the observed STIPs, particularly between semi-major axes of 5 and10 au (or 1000 < P < 1 × 104 d). In other words, STIPs could be the inner region of a morecomplex planetary system, with an outer system perhaps similar to the outer Solar System.Because the dynamics of a planet is not only affected by the host star, but also by the gravita-tional interactions with other planets in the system, deviations from expected orbits of knownplanets could reveal the presence of additional unknown planets. As an example, Neptune’sexistence was predicted mathematically by John Adams and Urbain Le Verrier due to the dis-crepancy between the predicted and the observed position of Uranus. Observations went on toconfirm the predictions. Nowadays, due to the increase in computing power, we can make useof numerical simulations to explore a wide range of possibilities and to investigate the dynamicsof planetary systems more generally.In this chapter, we investigate the observability and stability of STIPs in the presence of aplanetary outer perturber by using direct numerical integration along with secular theory. Wepresent case studies of two systems that have high multiplicity (Np > 5), namely Kepler-11 (K11) and Kepler-90 (K90)29 (Lissauer et al., 2014b; Cabrera et al., 2014; Schmitt et al.,2014). For each system, we first examine the planetary dynamical behavior using just the knownplanets, and then that of the known planets with an additional perturbing Jupiter-like planet.When the perturber is included, it is placed exterior to the last known planet in each system,initially at 5.2 au. The location of the Jupiter analogue respect to the STIP is expected tomodify the secular frequencies of the system, to explore a broader range of secular interactionswe “migrate” the planet inwards in our numerical simulations. Looking at the dynamicalevolution of the STIPs with a migrating perturber, we identify the perturber location thatmight be affecting the inner system through secular resonances, reflected in instability of theinner system or selective inclination changes on the inner planets. Both instability and selective29K90 is also known as KIC 11442793 and KOI-351.72inclination changes affect the observability of STIPs, the former reduces the total multiplicityand the latter can strongly affect the number of detectable transiting planets in a system(depending on the inclination change), both explainable by the existence of unknown planetcompanions. We verify the secular interaction, by following up each tentative configuration,i.e., STIP plus perturber location, with analytical and synthetic secular theory, as laid out inSect. 2.2 and 2.3.The work presented in this chapter has been published as Granados Contreras and Boley (2018).The original article has been adapted to the structure of this thesis, with a few modificationsand additions. In particular, the methodology is described in Sections 2.2 throughout 2.3. Inaddition, we include the phases of the synthetic secular theory instead of the sole precessionspectra.4.1 STIPs studied and simulations setup.The values used for K11 are those used in Sect. 3.1.1 (Table 3.1). The only difference is thatwe use Mp = 20.0M⊕ for K11g, instead of 8.0M⊕, because this value is closer to the uppermass limit found dynamically by Lissauer et al. (2013) and in this way the system resemblesK90 by having the outermost planet as the most massive in the STIP. The other system,K90 (Lissauer et al., 2014b; Cabrera et al., 2014; Schmitt et al., 2014), has seven planetsorbiting an M∗ = 1.2 ± 0.1M star. K90 (along with TRAPPIST-1) has the largest numberof confirmed planets thus far, making it one of the closest in planet multiplicity to the SolarSystem. However, unlike the Solar System, all of K90’s known planets are confined within 1au. The masses of K90’s planets have not been measured directly, and as a result, the massesused here are estimated from the mass-radius relation described by Wright et al. (2011), whichis based on the size distribution reported by Lissauer et al. (2011b). Table 4.1 summarizes themeasured or estimated properties for K90. The values in Tables 3.1 and 4.1 are adopted forthe present study, unless otherwise noted.Published eccentricities and inclinations are used when available; for the K90 planets the orbitaleccentricities are set to zero (e = 0). The inclination values given in Table 3.1 and 4.1 are relativeto a reference plane that was determined by averaging the published orbital inclinations relativeto a perpendicular to the system’s line of sight, which is 〈iK11〉 = 89.52◦ and 〈iK90〉 = 89.68◦.We create 100 realizations of each system, in which the longitudes and anomalies (Ω, $ andM) for each planet are drawn from a uniform random distribution between 0 and 2pi. Theserealizations are first run without the external perturber to establish the stability of the systemsover 10 Myr. We further test the stability of K11 and K90 in the presence of a Jupiter analogueand run the same initial conditions with an additional planet, also for 10 Myrs. The results ofthese simulations are discussed in detail in Sect. 4.2. Unless otherwise noted, the perturber has73a mass Mp = 1MJ and initial orbital elements ap = 5.2 au, ep = 0.05, ip = 1.3◦, ωp = 273.8◦,Ωp = 100.5◦ and Mp = 93.8◦, where the subindex P denotes parameters for the perturber.The realizations that include the perturber will be denoted as K11+ and K90+ throughout thischapter.The stability of the STIP is expected to be sensitive to the semi-major axis of the perturber(due to the secular frequencies), and as such, we systematically explore the semi-major axisparameter space of the perturber by forcing it to migrate inwards. The simulations that includethe inward migration of the perturber are performed using the initial conditions that resultedin stable configurations for K11+ and K90+ after 10 Myrs. For K11+ we used the full set ofstable realizations, but chose only one of the K90+ stable systems. This choice was based onthe analysis of the K11+ simulations, which all showed consistent behavior.N-body integrations were run using a modified version of the Mercury6 code (Chambers, 1999),which includes a general relativity correction as described in Sect. 2.1.1, which is appropriatefor low-eccentricity orbits (Nobili and Roxburgh, 1986). Because we want to resolve closeencounters should they occur, we use the hybrid integrator of Mercury6. When in MVS mode,the time step is set to 10−2 of the initial period of the innermost planet. Again, the totalintegration time, if not otherwise stated, is 10 Myrs.The inward migration of the Jupiter analogue is achieved by applying an acceleration term toPlanet P Rp Mp a I(days) (R⊕) (M⊕) (au) (◦)b 7.0082 ± 0.000019 1.31 ± 0.17 2.4 0.076 0.28c 8.7193 ± 0.000027 1.19 ± 0.14 1.7 0.088 0.00d 59.7367 ± 0.00038 2.87 ± 0.3 7.9 0.307 0.03e 91.9391 ± 0.00073 2.66 ± 0.29 6.9 0.424 0.11f 124.9144 ± 0.00190 2.88 ± 0.52 8.1 0.520 0.09g 210.6069 ± 0.00043 8.10 ± 0.8 69.1 0.736 0.12h 331.6006 ± 0.00037 11.30 ± 1.0 297.9 0.996 0.08M∗ = 1.2± 0.1M R∗ = 1.2± 0.1R Table 4.1: Nominal orbital elements of the known planets of Kepler-90 (K90) (Lissauer et al.,2014b; Cabrera et al., 2014; Schmitt et al., 2014). The uncertainties provided by Cabreraet al. (2014) for the orbital periods and planetary sizes are the 1σ error of their transit modeloptimization for each planet and their stellar parameter estimation from synthetic spectrafitting. Compared to Table 3.1, this table is missing the columns for e and ω because theseorbital elements have not been estimated from observations.74the planet of the formamig = − 2piτmig(1aua)[3(rˆ · r˙)rˆ + (r× r˙)× rˆ]= − 2piτmig(1aua)[2 (rˆ · r˙) rˆ + r˙] , (4.1)where r is the perturber’s position vector from the star, r˙ its velocity vector, a its semi-majoraxis, and τmig is the timescale for migration at 1 au. Over the region of interest, migrationis nearly constant. For these simulations, we set τmig = 30 Myr30, which is equivalent toa˙ ≈ 0.42 au Myr−1. During the 10 Myr migration, the orbital precession of the perturberincreases by 1 to 2 orders of magnitude. The perturber’s eccentricity also changes, decreasingamong the stable systems from 0.05 to 0.035 at aJ = 3 au and to 0.02 at aJ = 1 au, correspondingto t = 6 and 10 Myrs, respectively.We complement the numerical simulations with secular theory as laid out by Murray andDermott (2000) and described in Sect. 2.2. For systems in which the planets have small eand i, the theory can identify locations at which a test particle would be in an eccentricityor inclination secular resonance, although higher-order methods (e.g., Laskar, 1985, 1986) ordirect N-body calculations are needed to determine the dynamical effects of a given resonance.As mentioned in Sect. 2.2.1, if planets have strong interactions, such as might be expected nearmean motion resonances (MMRs), then additional frequencies can become important and/orthe predicted frequencies can become shifted with respect to the second-order theory. BecauseK11 and K90 exhibit TTVs, we expect deviations from second-order theory, as is the case inthe Solar System with Jupiter and Saturn (Brouwer and van Woerkom, 1950). Nevertheless,such effects are excluded in our calculations. We also exclude the effects of GR on the secularfrequencies. Regardless, as will be shown, secular theory seems to identify the locations ofresonances with reasonable accuracy.Using a given system configuration (e.g., the current STIP and perturber’s orbital elements) andthe secular code, we highlight the resonant structure in the system. First, we calculate the forcedeccentricities and inclinations. Resonances will correspond to distances at which the forcedeccentricity or inclination show very large increases, resulting from a singularity introducedwhen a precession frequency matches a system eigenfrequency. Using this approach, we cancompare the resonant structure from secular theory with the behavior of the N-body simulations.For example, if a given planet becomes orbit crossing or develops large inclination variations, webuild the secular theory for the system, but exclude the highly perturbed planet(s). As shownbelow, in K90(+), the two innermost planets are the most easily excited, so we examine theresonance structure by removing these planets from the secular theory, allowing us to see if theforced eccentricity and inclination at their semi-major axes suggest a resonance, at least for a30It will be shown that this migration timescale is about 800 times longer than the longest secular timescalein either K11 or K90 without the perturber.75test particle. This is only reasonable whenever the removed planets have low masses comparedwith the other planets. K90b and K90c are tightly coupled, so this must be done with somecaution. Nevertheless, as we will show, this approach is reasonable enough to be useful for thepresent situation. As with K90(+), we examine the secular resonant structure in K11(+) byremoving K11b, which is often the first planet to become orbit crossing, or by removing bothK11b and K11c due to their coupling.We also calculate synthetic frequency spectra for select realizations of K11(+) and K90(+)to further explore the secular dynamics of the systems and to examine the effects of strongplanet interactions. This is done by rerunning a given simulation for 0.3 Myr with output every60 years. The synthetic secular spectrum provides the principal apsidal and nodal precessionfrequencies and amplitudes.4.2 Unperturbed vs. perturbed STIPs: dynamical outcomesThe N-body simulations of K11 and K90 (without the perturber) allow us to examine thedynamical stability of these systems generally and to provide a point of comparison for thesimulations with the Jovian perturber. We find that the stability of K11 is very sensitive toperturbations to the nominal argument of pericenter for non-zero eccentricities. In particular,44% of the simulations become unstable using the nominal eccentricities, but random longi-tudes. This is qualitatively consistent with Mahajan and Wu (2014), who found that K11 ispreferentially unstable if any of the eccentricities e > 0.04, using the same planetary massesused here. In contrast, 100% of the realizations in this study are stable if the initial eccentric-ities are zero, which agrees with the zero-eccentricity stability tests of Lissauer et al. (2011a,2013). For K90, 20% of the our realizations present instability despite having initial orbitaleccentricities set to zero for all planets. This result is unexpected at face value, as we mightanticipate the system to be stable (as with K11) under such conditions. This agrees with thestability test presented in K90’s discovery paper (Cabrera et al., 2014). The system’s secularperturbations could be the origin of the noted instability and/or mean motion resonances, aswill be discussed later.The simulations of K11 and K90 are now used as a reference for exploring additional dynamicalperturbations. Using the same realizations for K11 and K90, we now include a Jupiter analoguein each system, the stability of the STIPs is not affected by the presence of the single perturberwhen it is placed at 5.2 au. The same 44% and 20% analogues become unstable in less than 10Myrs for K11+ and K90+, respectively, which suggests that the inner system instability is notaffected, at least for this short time period and for the initial placement of the Jupiter analogue.While the stability of the STIPs is unaffected, the systems do respond to the presence of theouter perturber. This is highlighted in Figures 4.1 and 4.2, which show the inclination evolution76Figure 4.1: Evolution of the orbital parameters of a stable K11 analogue (left) and K11+ (right)system. The top panel shows the inclinations for all of the nominal planets in the system, whilethe bottom panel shows the longitude of ascending node relative to K11g, the outermost planet.Figure 4.2: The same orbital parameters in Figure 4.1 are shown for K90 (left) and K90+(right). Dynamical rigidity is also observed in this system, i.e., all the longitudes of ascendingnode precess at the same rate. K90b and c (magenta and blue lines) are more tightly coupledto each other than to the rest of the planets.for the planets in stable realizations of K11 and K90, respectively, with and without the presenceof the Jupiter analogue. In each case, the planets exhibit variation among their inclinations withtime. However, when the perturber is included there is an additional large-scale variation inall of the orbital inclinations, i.e., each planet’s inclination oscillates with respect to a commonorbital plane, while the orientation of the shared orbital plane changes. The presence of theperturber precesses the longitudes of ascending node of all the inner planets at the same rate,which is also shown in Figures 4.1 and 4.2 (with ∆Ω). Such “dynamical rigidity” has beenseen in simulations of other systems when perturbers are introduced (Kaib et al., 2011) athigh inclination, and has been further explored analytically by Boue´ and Fabrycky (2014a,b).77Figure 4.3: K90+ system with a highly inclined perturber with an I = 50◦ and with a = 5.2 au.The dynamical rigidity is evident in the coherent variation of the inclinations. The subpanelis a close up of the inclination within 0.12 to 0.14 Myrs to show that the mutual inclinationevolution of the planets remains small due to the synchronous nodal precession of the planets,even though the common orbital plane oscillates.Figure 4.4: Orbital evolution of K11+ in the presence of a migrating outer perturber. Toppanel: each planet displays 3 curves: pericenter, q; apocenter, Q; and semi-major axis, a.Only the first 6.2 Myrs of the simulation are shown. The forced migration of the perturber,labeled as “Jup”, can be observed in orange. The dashed line indicates the perturber’s locationaJ = 0.98 au (at t ≈ 10 Myrs) at which the system becomes unstable, possibly due to secularperturbations prompted by the perturber. Bottom panel: the orbital inclination, I, is shown.While the system is stable, the precession rate of the common orbital plane increases as afunction of the proximity of the perturber.78Figure 4.5: Orbital evolution of K90+ in the presence of a migrating outer perturber. Theorbital elements are the same as those in Figure 4.4. The precession behavior of K90+ issimilar to K11+, with K90+ becoming unstable when aJ = 2.5 au (at t ≈ 6.3 Myrs). However,between perturber orbital distances of approximately aJ = 2.5 au and 3.2 au, K90b and cbecome excited together onto a second plane of higher inclination relative to their originalorbits. The dotted line in the top panel indicates the perturber’s position at which K90b andK90c inclinations become about 10◦ larger than the other planets.Test simulations of K11+ and K90+ with the Jupiter analogue at high inclination (50◦), whichshowed that the longitudes of ascending node remained locked, with the STIP undergoing largeand coherent inclination variations (Fig. 4.3), the eccentricities remain small throughout thelength of the simulation. We note that the amplitude of the i oscillation is approximately twicethe perturber’s inclination, although these amplitude is smaller in K90+ which has a largertotal planetary mass than K11+, as shown in Figures 4.1 and 4.2.We have thus far only considered one location for the perturber, arbitrarily introduced at 5.2 au.We now explore a wider range of semi-major axes by forcing the perturber to migrate inwardsfrom 5.2 au to about 1 au (Eq. 4.1). First, the overall behavior of the STIPs is unchangedby simply moving the Jupiter analogue inwards. Figures 4.4 and 4.5, for example, show thatthe systems continue to exhibit coherent changes in the orbital plane, but with an increasingprecession rate of the plane as the Jupiter analogue approaches the inner system. However, forcertain system configurations, the outcomes can be very different. In most cases, this appearsto be due to shifts in the secular frequencies. For example, eccentricity resonances can drivethe STIP towards instability, causing planetary collisions (as occurs for K11). An inclinationresonance can force a planet (or planets) out of the original common orbital plane, effectivelycreating multiple orbital planes in a stable system (as occurs in K90). Specifically, runs with an79initially stable K11+ develop an instability if the perturber reaches a semi-major axis of 1.0 au,which places the perturber in a 3:1 near-MMR with K11g. This commensurability could becausing the instability in K11+, but the secular resonances could also have a strong contributionto the destabilization of K11+, as will be discussed shortly. At this time, 59% of the unstableK11+ realizations result in either K11b crossing the orbit of K11c followed by K11f crossingK11e’s orbit or vice versa.31 In K90+, when aJ . 3.2 au, K90b and K90c evolve together awayfrom the rest of the planets by about 16◦ in inclination, with the system remaining stable.4.3 Secular theory vs. synthetic secular theoryThe behavior of both K11+ and K90+, observed in Figures 4.4 and 4.5, can be understoodby looking at their forced eccentricities and inclinations. Figures 4.6 to 4.8 show the forcedinclination (red solid line) and eccentricity (black dashed line) of K11+ and K90+ at selectedlocations of the perturber. As the Jupiter analogue moves inwards, the locations of the inner-most inclination and eccentricity resonances move outwards in the STIP, eventually crossingthe innermost planets. For K11+, an inclination and eccentricity resonance overlap the semi-major axis of K11b when aJ = 0.98 au (Fig. 4.6), which is when instability occurs in thecorresponding simulations. For this analysis, K11b and K11c were removed from the seculartheory calculation, effectively treating them as test particles.In the case of K90+, when the perturber is at aJ ≈ 3.3 au the two innermost planets appear to betrapped in an inclination secular resonance, as determined by secular theory with planets K90band K90c removed (Fig. 4.8). Together, Figs. 4.5 and 4.8 suggest that the separation of theK90+ system into two distinct orbital planes (for aJ . 3.2 au) is due to an inclination resonanceoverlapping the K90b and K90c positions. This secular resonance excites the inclination ofboth planets, which are strongly coupled, and its strength increases as the perturbing planet’sdistance to K90 decreases. This second orbital plane for K90b and K90c can acquire a maximuminclination ibc = 18o if aJ ≈ 3.0 au. While there is also an inclination resonance overlappingthe location of K11b in K11+ (when aJ ≈ 0.98 au), an eccentricity resonance is also present(Fig. 4.6), possibly leading to instability before large inclination changes can occur.Removing planets b and c from both K11+ and K90+ for the secular calculations assumesthat the planets are massless, which is not true. As such, the method is not guaranteedto highlight resonances, although the agreement with the N-body simulations suggest that theapproximation is valid in this case. Removing the planets nonetheless reduces the complexity ofthe frequency structure in the analytic theory. To highlight this, we show in Fig. 4.7 a secularmap in which we only remove K11b, allowing K11c to contribute to the secular model. Aninclination and eccentricity resonance is present at r ≈ 0.09 au, just inside K11b. The position31K11f crossing the orbit of K11e first and then K11b crossing K11c’s orbit.80Figure 4.6: Secular map of the forced inclination (red solid line) and eccentricity (black dashed line) of K11 in the presence of a Jupiter-like perturber and excluding the planets in parenthesis of the calculations, i.e., K11b and K11c. Two different semi-major axes of theperturber are shown in the top and bottom panels. The locations of the inner planets are shown by vertical dotted lines. Top panel: thegas giant is at aJ = 5.2 au, and none of the resonances coincides with the location of the inner planets. In contrast, when aJ = 0.98 au(bottom panel), inclination and eccentricity resonances are located at the position of K11b. In this case the eccentric resonance appearsto contribute to destabilizing the system.81Figure 4.7: Similar to Fig. 4.6, but excluding only K11b (as indicated in the plot). As before, the panels show the results for two differentperturber locations. There are multiple inclination and eccentricity resonances just interior to K11b’s position. Due to the tight couplingbetween K11b and K11c, Fig. 4.6 might better reflect the onset of instability. Nevertheless, these panels highlight the richness of thesecular structure of the inner system.82Figure 4.8: Secular map of the forced inclination and eccentricity of K90 in the presence of a Jupiter-like perturber. Top panel: theperturber is at aJ = 5.2 au, and none of the resonances coincides with the location of the inner planets. When we consider aJ = 3.0 au(bottom panel), a wide inclination resonance is located at the position of K90b and K90c. This resonance increases the inclination ofboth K90b and K90c, without affecting the stability of the system.83of this resonance is fixed regardless of the perturber’s proximity to the inner system, suggestingthat the location of this resonance is due to the K11+ inner planets. There are additionalresonances that are at even smaller semi-major axes, which are affected by the perturber’slocation. As the perturber moves inwards, this inner resonance structure moves outwards andthe resonance wings overlap. At face value, Fig. 4.7 suggests that K11b might not enter a secularresonance because the resonance just inside K11b’s semi-major axis does not change locationwhen different semi-major axes are explored for the perturber, not even when aJ = 0.98 au.We will later show that K11b and K11c are strongly coupled and Fig. 4.6 might better reflectthe outer system’s influence on K11b and K11c. Nevertheless, Fig. 4.7 shows the proximity ofthese resonances to K11b, which could explain the observed instability of the system and itsapparent fine tuning even when the perturber is absent, because this particular feature seemsto be an intrinsic property of the inner system. A small perturbation in the configuration ofplanets could cause K11b to enter a secular eccentricity resonance and collide with K11c, themost common outcome in our unstable simulations.We investigate the secular structure of the systems further by calculating the synthetic secularfrequencies (precession spectra) with and without the perturber. The precession spectra aredetermined directly from the N-body calculations, using 0.3 Myr of output with a sample timeof 60 years. The total integration time limits the resolution of the lowest frequency, in this casethe lowest frequency is 0.0012◦ yr−1 or 4.3′′ yr−1. In K11+ and K90+, the perturber is set toaJ = 0.98 au and aJ = 3.0 au, respectively. The synthetic precession spectra and phase anglesare compared with the apsidal and nodal secular eigenfrequencies in Figures 4.9 through 4.12.We note that in these four figures, the scaling of the x-axis is divided in two. For frequenciesshorter than 1◦ yr−1, we use a logarithmic scale, and linear otherwise. As expected, the valueand number of eigenfrequencies change from K11 to K11+ and K90 to K90+, respectively,due to the additional planet in the perturbed systems, with a greater displacement for theapsidal eigenfrequencies. The amplitude of the global Fourier spectra, in both eccentricity andinclination, for K11+ and K90+, is higher than in the unperturbed analogues; there is about oneorder of magnitude difference in K11+ and two to three in K90+. The nodal eigenfrequenciesshow very good agreement with the synthetic spectra for both K11(+) and K90(+), althoughthere are some deviations. The apsidal precession frequencies of both systems exhibit strongerdisagreement between the synthetic spectra and the secular theory, but the features in theFourier spectra remain recognizable, particularly for K11(+). The misalignment of the principalapsidal frequencies relative to the estimated secular eigenfrequencies might indicate that thereare strong interactions among the planets, which are neglected in the secular calculation. Thiscould include non-linear secular resonances and near-mean motion resonance (MMR). Malhotraet al. (1989) showed that the e−$ eigenfrequencies are significantly shifted if a near first orderMMR is present in the system because the averaged effect of the term e cos(jλ− (j+ 1)λ′+$)in the secular expansion is non-zero. There is a weaker effect of such a near-resonant term on84Figure 4.9: Comparison of the theoretical and synthetic apsidal precession frequencies and phases of K11 (left) and K11+ (right). Notethat the frequency axis is divided in two different scales; the frequencies < 1◦ yr−1 are scaled logarithmically, and linearly otherwise.Top panels: Apsidal precession spectra, where at least two of the principal frequencies are offset from the nearest eigenfrequency. Thediscrepancy between the theoretical and synthetic principal frequencies indicates that an unaccounted eccentricity secular non-linearresonance is present in K11+. Bottom panels: Phase angle of the periapses of K11 and K11+ planets. The periapses of both systemsare preferentially out of phase and circulating for high precession frequencies, except for K11b, K11c and K11g (on the left). K11b andK11c exchange phases at shorter precession periods.85Figure 4.10: Comparison of the theoretical and synthetic apsidal precession frequencies and phases of K90 (left) and K90+ (right). Thenotation used in Fig. 4.9 is adopted here. For these configurations, there is aliasing at higher $˙ (short periods) and the characteristicprecession frequencies miss most of the estimated eigenfrequencies. The aliasing reflects in the phase of the periapses, for both K90 andK90+ systems. In K90, K90b and K90c exchange pericenters at different timescales. In K90+, for most planets and timescales, thepericenters circulate, with the exception of K90b, K90c and the perturber.86Figure 4.11: Comparison of the theoretical and synthetic nodal precession frequencies and phases of K11 (left) and K11+ (right). Toppanels: Nodal precession spectra, there is a closer alignment of the principal precession frequencies with the eigenfrequencies than in theapsidal case. Bottom panels: Phase angle of the nodes; in both systems it is clear that there is some degree of dynamical rigidity, ornodal locking. This locking is probably induced by K11g in K11. When the perturber is present, it causes the nodal locking. K11b andK11c nodes circulate at high nodal precession frequencies, indicating that their nodes do not align on short timescales.87Figure 4.12: Comparison of the theoretical and synthetic nodal precession frequencies and phases of K90 (left) and K90+ (right). Thenotation used in Fig. 4.11 is adopted here. The characteristic nodal precession frequencies are in close alignment with the eigenfrequenciesand the aliasing is minor to moderate. Similar to that seen in in Fig. 4.11, the nodes seem to be locked, which is evidence of dynamicalrigidity. In the case K90, K90h seems to be the culprit. The antialignment of K90g might be a result of a metastable configuration. InK90+, as expected, the perturber (in orange) originates a “stronger” rigidity within the STIP, as shown by the aligned planetary nodeangles.88i − Ω because the nodal term is associated with i2. Having this in mind, we will search forpossible first-order MMRs in K11 and K90 later in this chapter (Sect. 4.3.1). The MMRs canalso add eigenfrequencies to the system, as occurs in the Solar System due to the 5:2 MMRbetween Jupiter and Saturn (Brouwer and van Woerkom, 1950).Like the precession spectra, the apsidal phases are not as coherent and well defined as thenodal phases. Large phase changes occur at the principal precession frequencies, with thetypical phase changes being near ±pi or an exchange of phases between two or more planets.The latter is often seen between the two innermost planets in K11(+) and K90(+), for both inthe nodes and pericenters.The study of the apsidal phases of a system as a function of precession frequency indicate thealignment of the pericenters and the correspondent timescales. For planets in or near MMRs,this alignment can determine whether the configuration is stable, which also depends on theeccentricity and closest approach of the two planets at conjunction (Murray and Dermott, 2000,Ch. 8). We will describe briefly the apsidal alignments and frequency ranges where these areencountered. The further stability study as a function of pericenter conjunctions is out of thescope of this dissertation. The apsidal phases of K11(+) and K90(+) have similar character-istics. In K11 (left panels in Fig. 4.9), the periapses are aligned at very short frequencies,but change phase drastically when |$˙| > 0.01◦ yr−1. For |$˙| > 1.0◦ yr−1, planets K11b andK11c become antialigned for |$˙| > 1.5◦ yr−1, with the latter being aligned with K11g. Theright panel shows the pericenter phases for K11+ in which the perturber is displayed in orange.There is no such preferred phase alignment in this case. In Fig. 4.10, which shows the resultsfor K90(+), the periapses of K90 change phase even more notoriously than in K11 (Fig. 4.9),with a slight alignment for |$˙| < 0.05◦ yr−1. There is also clear phase swapping between K90band K90c. In K90+, the apsidal phases do not show coherent structure, except for K90b, K90cand K90J, which are aligned for |$˙| > 0.1◦ yr−1.The nodal phases provide information about the orientation of the nodes within a STIP, andcan be used to determine whether a planetary system shows dynamical rigidity. The node ofK11g, in K11, is approximately pi rad out-of-phase with respect to the rest of the planets forshort precession frequencies (see left bottom panel Fig. 4.12) within the region |Ω˙| < 0.2◦ yr−1.At higher nodal frequencies, the node of K11f aligns to K11g’s making them both antialignedto the nodes of the other planets in the system. This major nodal conjunction within a STIPis an indication that one of the planets in the STIP is locking the nodes, in this case K11g. Incomparison, the nodes of the inner planets of K11+ are aligned with each other but pi rad awayfrom the perturber’s node at frequencies |Ω˙| < 1.5◦ yr−1 and |Ω˙| ≥ 0.8◦ yr−1. The nodes ofK11b and K11c circulate at precession frequencies |Ω˙| > 1.0◦ yr−1. K90 nodal phases (Fig. 4.12)are difficult to interpret for precession frequencies below |Ω˙| < 1.0◦ yr−1. Above this thresholdthe planetary nodes are aligned, except for K90g, which is antialigned to the other planets.K90h seems to drive the nodes of most of the planets. Something similar occurs with the nodes89Figure 4.13: Zoom-in of the nodal principal component of planet b in K90 (red line) and K90+(black line). In this case, K90+ was evolved for 0.4 Myrs, which increased the resolution of theFourier transform. The K90+ dominant peak on the top right of Fig. 4.12 splits into two peakswhen the resolution is increased. The vertical black dashed-lines represent the eigenfrequenciesfor K90+ in the given range, while the red dotted line corresponds to the eigenfrequency forK90 in that range.of K90+, but in this case all nodes in the STIP are aligned with K90h (including K90g) andmisaligned to Jupiter’s node. This is another example of internal node locking within a STIP.The main precession period of each planet, in K11(+) and K90(+), and the correspondingamplitude are displayed in Table 4.2. In the four systems, the two innermost planets have thesame main precession period ($˙ and Ω˙) and their precession spectra are almost indistinguish-able, confirming the strong coupling of these planetary pairs. Furthermore, in K11, the nodesof the six planets are dominated by the same period, suggesting that K11g is dominating thenodal precession, which causes a dynamical rigidity in the orbital plane. On the other hand,there are two distinct pericenter precession periods, in which K11b and K11c have an apsidalperiod of PK11b−c ' 1346 yr and K11d through K11g have periods of PK11d−g ' 33340 yr. Theperturber in K11+ changes the frequencies and decouples K11g; K11g’s node precesses at thesame rate as the perturber, which is around 4 times faster than the rest of planets in K11+. Incontrast to K11, the main nodal precession periods of K90 do not suggest an intrinsic dynamicalrigidity, although there is strong coupling between K90b and K90c, as well as K90g and K90h.However in K90+, there are two distinct nodal precession periods that are different by only7%. These nearly overlapping eigenfrequencies (Fig. 4.12) suggest that a secular resonance iscausing the large change in inclinations for K90b and K90c (Fig. 4.5). The amplitude of thebroad frequency peak seen in Fig. 4.12 dominates the spectra by one order of magnitude ormore. When the K90+ simulation is allowed to evolve to 0.4 Myrs,32 it is possible to identify the32The integration time is limited by instability of the system after 0.4 Myrs90Planet K11 K11+ K90 K90+Pap e Pap e Pap e Pap e(yrs) (yrs) (yrs) (yrs)b 1346 0.0225 1266 0.0276 104 0.0008 6667 0.0028c 1346 0.0201 1266 0.0246 104 0.0011 6667 0.0026d 33340 0.0131 3615 0.0158 2703 0.0007 2655 0.0110e 33340 0.0146 3615 0.0157 138 0.0007 1230 0.0078f 33340 0.0194 3615 0.0142 205 0.0012 695 0.0050g 33340 0.0988 2190 0.0908 205 0.0121 216 0.0047h · · · · · · · · · · · · 205 0.0026 6667 0.0039J · · · · · · 60012 0.0436 · · · · · · 60006 0.0398Pnode i Pnode i Pnode i Pnode i(yrs) (◦) (yrs) (◦) (yrs) (◦) (yrs) (◦)b 12503 0.3815 9377 1.0840 21431 0.1738 23079 8.8099c 12503 0.3777 9377 1.0728 21431 0.1732 23079 8.7856d 12503 0.3492 9377 0.9740 2831 0.0831 21431 0.4781e 12503 0.3374 9377 0.9362 1288 0.0988 21431 0.4999f 12503 0.2893 9377 0.8085 785 0.1174 21431 0.4982g 12503 0.2253 2084 1.3211 480 0.0864 21431 0.4880h · · · · · · · · · · · · 480 0.0164 21431 0.4829J · · · · · · 2084 0.0525 · · · · · · 21431 0.3145Table 4.2: Main precession period and amplitude of each planet in K11, K11+, K90, and K90+.Pap and Pnode denote the apsidal and nodal period, respectively. The values of e and i providedhere are the maximum amplitudes from Figures 4.9, 4.11, 4.10 and 4.12 for each planet, andwhenever pertinent, the perturber. The apparent switch of the nodal precession from planetsK90b and K90c to K90+d through K90+h is an artifact of the Fourier frequency spacing, alongwith a modest change in the eigenfrequency structure. Comparing the tabulated periods amongthe planets in a system helps to highlight which planets exhibit tight coupling.91two close frequency components in this dominant broad peak, as expected from secular theoryfor K90+. These are shown in Fig. 4.13, with the peaks corresponding to Pnode ' 23100 yr andPnode ' 21400 yr, which straddle the eigenfrequency associated with K90 in that range. Theamplitude of these two frequencies varies when the spectra is calculated from a shorter or longerintegration, which complicates the identification of the main precession frequency among thesetwo. The whole precession spectra of the longer N-body integration is not shown in Fig. 4.13,due to severe aliasing at large frequencies.The identification of the frequency components in the apsidal precession spectra of K90 andK90+ (top panels in Fig. 4.10) is complicated due to aliasing. Even so, we can discern asubstantial misalignment of the main frequency components compared with the nearest eigen-frequency, as mentioned above. In K90, there is a double-peaked feature corresponding to K90band K90c at $˙ ∼ 0.05◦ yr−1. This is absent in the spectra of the other planets, which couldoriginate from a non-linear secular resonance or from an MMR.4.3.1 Mean motion resonances MMRsThe discrepancies between the precession spectra of K11 and K90 and their secular eigenfre-quencies demonstrate that either additional physics should be included in the secular theoryor that we need to consider a higher order expansion in e, i, and mass. Some of the physicaleffects that we know we are neglecting in the secular code include GR and MMRs. The GRcontribution to the precession rates is small; it only affects the apsidal precession frequenciesof the synthetic spectra by . 2%, with K90 being the most affected. GR does not shift thenodal spectra in either system. A MMR, depending on the order, can affect linear terms ineccentricity and inclination.K11 and K90 exhibit TTVs (Lissauer et al., 2011a; Cabrera et al., 2014), which occur if thereare strong gravitational interactions between the planets. The discovery papers of both K11 andK90 also discussed MMRs for different planet combinations. The period ratios of K11b-K11c,and K90b-K90c suggest a 5:4 near commensurability in both systems. In addition, Cabreraet al. (2014) discussed a possible 2:3:4 near resonance between K90d, K90e and K90f. Theseinitial claims of MMR in the systems were based on the period ratios of neighboring planets, butdynamically an MMR might not occur, due to the resonance’s dependence on the eccentricityand inclination. To confirm the existence of an MMR in a pair of planets, it is necessary toestimate their resonant angle.92Figure 4.14: Librating resonant angles in K90. In blue is the interaction with the externalplanet, ϕ′, and in red with the internal, ϕ. For clarity, the left panels show the time evolutionof the resonant angle only within the first 40 kyr of the simulation; while the right panel (verticalhistograms) summarize the behavior of ϕ′ and ϕ over 0.3 Myrs. The top panel correspondsto the 5:4 near-MMR between K90b and K90c. The bottom panel shows the resonant anglebetween K90g and K90h in a 3:2 near-MMR.The two resonant angles for a first-order MMR between neighboring planets can be calculatedfrom the N-body simulations usingϕ′ = jλ′ + (1− j)λ−$′, (4.2)ϕ = jλ′ + (1− j)λ−$, (4.3)where λ is the mean longitude and $ is the longitude of pericenter. The prime here denotes theexternal planet. If the planets are in an MMR, one or both resonant angles will librate. Usingthe given masses and orbital elements for K11, we find that all the resonant angle combinationscirculate at all times. Thus, we cannot confirm a 5:4 near-MMR between K11b and K11c.In K90’s case, two pairs of neighboring planets display a librating first-order resonant angle(Fig. 4.14). These plots suggest that K90b is in 5:4 MMR with K90c and that K90g is in a 3:2MMR with K90h. In both cases, the resonant angle librates with large-amplitude around 0◦ and180◦. The existence of the 5:4 resonance between K90b-c suggests that some migration couldhave taken place. The possible 3-body resonance between d, e and f noted by Cabrera et al.93(2014) would correspond to a resonant angle ϕ = 2λd−6λe +4λf . After plotting this angle fromour simulations, we find that the resonant angle circulates through the length of the realization.The 3:2 MMR between K90g and K90h could be contributing to the instability present inthe K90 simulations with secular resonances that were not accounted in our secular theorycalculations. This instability would depend on the exact masses of the planets, particularlythose for K90g and K90h (Cabrera et al., 2014), which are the most massive. Although, thenominal masses of K90 planets are unknown, the impact of these massive planets on the systemstability could be strengthened by their apparent MMR.The behavior of the K11 and K90 analogues explored here show a range of dynamical outcomesthat STIPs could have in the presence of an outer (undetected) perturber. To date, it is uncer-tain whether the detected planets among the known exoplanet systems represent a complete setor whether there are undetected planets in those systems. Follow up measurements of orbitalinclinations of STIPs, using transit duration variations (TDVs), could reveal dynamic rigidityin the orbital plane and as a consequence, could be used as an indirect detection method foran outer planet. The observational determination of the longitudes of the nodes of planets ina given STIP could tell us whether there is an external perturber should the ∆Ωs for all theplanets be . 20◦ (this value coming from Figs. 4.1 and 4.2). In a non-perturbed STIP, thenodes would not preferentially align to a specific value, although it might be the case that asubset of the planets present nodal alignment due to coupling. Moreover, an external perturberwill cause a coherent change in the inclination of the planets, causing the duration of theirtransits to change coherently. An outer planet could also cause large breaks in the nominalorbital plane, which could explain some spin-orbit misalignment of low mass planets or reducethe number detected by transits.94Chapter 5DiscussionIn Chapters 3 and 4 we discussed the results of the respective individual studies. Here, weexpand on this discussion, considering the methodologies and results of the work as a whole.We also highlight caveats and future directions throughout the chapter.5.1 Close-in giant planets as an extreme outcome ofplanet-planet scattering.In the simulations that explored the metastability of STIPs and corresponding planet consol-idation, the outcomes were very different between simulations with and without gaseous tidaldamping. Undoubtedly, this damping has an important effect on the stability of the systems.The initial conditions that we used in these studies also likely played a key role in determin-ing the results. In particular, systems that were K11 analogues had a higher mass than thesix-planet synthetic systems, which was an unintended consequence of the experimental design.Specifically, the K11 analogues contained 30.1M⊕ in planets, while the six-planet systems ex-plored with tidal damping had on average a much smaller total planetary mass. Thus threecritical cores would form in principle in the K11 analogues, while typically only one, and inmany cases none, could form in the synthetic systems (Fig. 3.2). This inconsistency betweenstudies was introduced in an attempt to use random planetary masses within a specific rangethat also did not already have essentially a critical core. This approach would seem to bereasonable. However, a multiplicity larger than six should have been used to increase the rangeof possible STIP masses. Moreover, gaseous tidal timescales are inversely proportional to aplanet’s mass, causing small planets to damp quickly. This could justify having a large numberof small planets packed tightly together as an initial condition for proto-STIPs in gas.While not conducted as part of this work, the project’s outcomes could be readdressed by fixinga total system mass, and making many realizations of planet masses and multiplicity to buildsynthetic STIPs. The total mass could then be associated with a surface density profile inside agiven semi-major axis, e.g., 1 au, or orbital period, particularly if a non-solar star is consideredin the calculation.In the synthetic systems, planetary spacings were determined using a fixed number of mutual95Hill radii. Because the planetary masses were randomly assigned, the physical separationbetween the planets varied as a result. Even though actual planetary spacings in multiplesystems follow a distribution with a mean between 12 and 15 mutual Hill radii, the fixed Hillseparations were used here for simplicity and also to investigate how the spacings at the end ofthe simulations would compare with those at the beginning. It is clear that other methodologiescould be used. For example, both the initial masses and Hill separations could be randomlyassigned through, e.g., uniform, Gaussian, or Rayleigh distributions that are, unfortunately,unconstrained at this time.While the current simulations require further tuning to properly address the formation of crit-ical cores through STIP instability, they nonetheless reveal an interesting result for the perioddistributions of the planets. Figures 3.6, 3.7 and 3.9 show that the system evolves to a similarperiod ratio distribution regardless of the initial separations. Moreover, the structure in thedistributions near first-order MMRs has some similarities to the actual period ratio distribu-tion (Fig. 1.2) in that the exact commensurability is depleted, followed by a spike in densityat larger period ratios. It would be worth verifying whether we can reproduce the period ratiodistribution from our simulations using other initial spacings strategies, planet multiplicities,and disk profiles. The results so far lend support to the idea that the features near MMRs in theactual period ratio distribution are due to dissipative processes. In this case, eccentricity exci-tation from mutual gravitational perturbations followed by gaseous tidal damping are scullingthe distribution. The apparent universality in the period ratio distribution at the end of thesimulations suggests that the simulated and known distributions should match, particularly ifthe gaseous tidal damping is the main driving planet-disk interaction. The evident discrepancybetween the observed period ratio distribution (Fig. 1.2) and the one found in Chapter 3 (Fig-ures 3.6, 3.7 and 3.9) points to non-dissipative forces acting on the observed planetary systems,either before or after gas dispersal.Regardless of the initial mass distribution, critical cores did consolidate in the simulation setswith very tight initial conditions, i.e., small mutual Hill spacings of f = 1.5 and f = 3.0. Wetake this apparent efficiency result with caution because the initial separations were within thefeeding zone of planetary embryos, and therefore strong instability was expected. Nonetheless,we note that the instability was not halted instantaneously. The unstable systems had ap-proximately 100 orbits to consolidate into a lower multiplicity configuration with more massivecomponents. The period ratio distribution of these sets preserve the general structure of themore initially spaced experiments (f = 5, 9.5) and the peak delimited by the 3:2 and 2:1 res-onances remains significant. The mutual interaction of among planets and the disk structuresthe separation of the planets, even in unstable configurations.We are aware that the gas surface density used, i.e., a MMEN disk, might not be representativeof primordial protoplanetary disks. It was used in our simulations in the absence of a betterconstraint. Other parametrizations of the surface density could be used in subsequent studies.96It would be interesting to test the values used by Dawson et al. (2015), which found that highersolid surface densities form rocky planets on shorter timescales. This correlates with a highertotal planetary mass within a given period range.At the end of the simulations with tidal damping, the planets are on coplanar and nearly cir-cular orbits (e . 3×10−4). When these systems are evolved to much longer timescales withoutgas, they remain stable. If, on the other hand, the planets are given small eccentricity andinclination perturbations (Sect. 3.3.1), then the systems again exhibit metastability, showingsudden transitions to an unstable state. Should a mechanism exist, whether internal or exter-nal, which prevents the planets from becoming perfectly circular, then instability rates couldbe higher that found in our simulations, even in the presence of gaseous tidal damping. Whilewe have so far only included gaseous eccentricity and inclination damping, there are gas-planetinteractions that could cause eccentricity excitation (Ben´ıtez-Llambay et al., 2015; Eklund andMasset, 2017). In future studies, such excitation mechanisms could be added to IAS15, andshould be explored further, including the parametrizations necessary for capturing the effects.Regardless, if instability does set it, then the formation of critical cores should be treated care-fully, as dusty/metal-rich environments could cause significant delays in runaway gas accretion(Lee et al., 2014).An unexpected outcome of our simulations was the formation of coorbital planets. The occur-rence rate is low, and a few limited studies have also seen the formation of coorbitals. While itmight be easy to dismiss these configurations as flukes of a simulation, should a highly dissipa-tive environment be present during the early evolution of STIPs, coorbital planets are physicallypossible. A big question to answer here is whether we are capable of detecting coorbital planetsshould they exist. We pointed out that several detection techniques for coorbital planets havebeen suggested and some surveys have also been carried out, although no coorbital planets havebeen found so far. However, many of these studies have focused on the tadpole configurationrather than a horseshoe. This has stability reasoning behind it, because tadpole orbits areusually longer-lived in simulations than are horseshoe orbits. The same situation occurs in theSolar System: a larger number of low mass objects are known to orbit in or near the L4 andL5 Lagrange points than exist in horseshoe configurations. A deeper study of horseshoe orbitsis therefore necessary.A short test of the stability mass limit for horseshoe orbits (Laughlin and Chambers, 2002)showed that the existing threshold in the literature may be incorrect, albeit of the right order ofmagnitude. According to previous results, two Saturn-mass planets in a horseshoe configurationaround a solar-mass star should be stable over long timescales. Preliminary test with Mercury6,on the other hand, resulted in planetary collisions when using initial conditions of previousstudies. We did not include this test in the results of Ch. 3, but it could be easily tested andconstrained in future work.97More research for understanding transit lightcurves of coorbital horseshoe planets is needed. Ananalysis similar to that used in studying transit timing variations (TTVs) could help to identifyplanets in this configuration, especially if results have a resemblance to Fig. 3.13. Success usinga TTV analysis would depend on the horseshoe period of the orbits and the observing time.In this sense, it is important to determine an approximate relation for the horseshoe periodand the stellar and planetary parameters so that one can determine what type of horseshoeconfigurations might be detected with a given survey.Coorbital planets are unlikely to appear in the period ratio distribution, although, as mentionedin Ch. 1, we found six pairs of exoplanets with period ratios close to the 1:1 MMR. Half of thesepairs are only candidates and further assessment is challenging. More observations are neededof the K0624, K00521 and K02248 systems. Among the other three pairs, two are thoughtto orbit different stars (i.e., Kepler-132b and c, and Kepler-271 d and b, Lissauer et al. 2012;Rowe et al. 2014; Morton et al. 2016). The remaining pair, Kepler-1625 b and Kepler-1625 bI (Teachey et al., 2018; Teachey and Kipping, 2018; Heller, 2018), has been claimed to be thefirst observed exomoon. However, careful examination of their folded lightcurves shows thatthey are similar to the bottom panel of Fig. 3.16, particularly if the period is mistaken to bedouble its real value. This can occur in folded transit lightcurves when coorbital planets areat their maximum separation: in the periodogram the highest amplitude corresponds to twicethe physical period. The resulting lightcurve folded to the wrong period could look similarto what one might expect for a moon. We could determine whether the published lightcurveof Kepler-1625 b is due to a moon or a coorbital planetary companion by running N-bodysimulations of a moon orbiting a planet and generating the correspondent transit lightcurve,akin to what was done in Sect. 3.4. The comparison between the folded transit lightcurvesof coorbital planets with that of a planet-moon system should show that in the planet-moonsystem the phase of the second transit does not vary significantly, and it is zero on average withrespect to the planetary transit. As was shown, the phase of a coorbital planet pair exploresmost values except for a ∆φ centered at zero.5.2 STIPs and external planetary perturbers: effects onobservability and stability.The dynamical rigidity of planetary systems has also been observed in the N-body simulationsof 55 Cancri (Kaib et al., 2011; Boue´ and Fabrycky, 2014b) and HD 20794 (Boue´ and Fabrycky,2014b). Boue´ and Fabrycky (2014b) further studied in detail the theory behind the dynamicalrigidity of planetary systems in a hierarchical setting, in which the host star is part of a binaryor an outer giant planet is present. They particularly explored the conditions needed to drivespin-orbit misalignments (Boue´ and Fabrycky, 2014a).98Hansen (2017) also found dynamical rigidity while working on the hypothesis that the KeplerSingle Tranet33 Excess (KSTE)34 could be explained by secular resonances driven by long-period giant planets. The KSTE is the fractional surplus of single transiting planets based onthe expected fraction determined from the systems with multiple transiting planets. Hansen(2017) found that a fraction of the excess could be explained by a mixed population of Jovianand Saturn analogues, but only if they had high inclinations. They used in their numericalsimulations a selection of prototype planetary systems from Hansen and Murray (2013), whichincluded systems with a multiplicity of 3 to 10 planets and masses from 1 to 10 M⊕. The massweighted semi-major axis for the systems was 〈a〉M = 0.26-0.5 au, and external perturbers wereplaced between 1 and 5 au.We present a case study of two Kepler systems with high multiplicity, which is complementaryto Hansen (2017). We demonstrate that it is possible to drive high inclinations of low-massplanets through secular resonances without a highly inclined perturber and confirm that theeccentricities of a STIP are not necessarily affected or interchanged with the high inclinations,as occurs in the Lidov-Kozai effect (Kozai, 1962; Lidov, 1962), driven by a massive planetaryouter perturber. In our particular case, the inclination resonance excited the orbit of the twoinnermost planets in K90+; a single planet could be driven to a similar outcome through secularresonances. It is unclear where in a STIP a secular resonance will develop due the presenceof outer planetary companions; more case studies are needed to understand these dynamicaleffects.A subset of planets in a system having high mutual inclinations has observational implications.For example, the transit method is biased toward detecting planets that have low inclinationswith respect to the observer’s line of sight. In the case of K90+, the transit method would detecteither two or five planets instead of all seven. An RV observation of the system could make thedetection of the out-of-sight planets possible, because the independent velocity amplitude ofthe two innermost planets would only decay by about 5% (I ≈ 20◦), depending on the overallalignment of the system. In this particular case, the convolution of RV amplitudes is governedby K90g and K90h due to their mass. A quick analysis of the RV observability of K90+ andsimilar systems could be done by extracting the stellar radial velocity of the N-body simulationsand calculating the amplitude of the signal at different observational inclinations.In these simulations, the Jovian perturber was placed on a low eccentricity orbit (e ≤ 0.05).If the perturber’s orbital eccentricity were to be increased, which affects the width and powerof secular resonances, we would expect additional instability. For example, Clement and Kaib(2017) showed that increasing the nominal eccentricity of Jupiter and Saturn by a factor of twoenhances chaos within the inner Solar System and reduces the system’s stability.33Transiting planet (Tremaine and Dong, 2012).34Also known as the Kepler dichotomy.99As a first step, only a single Jovian perturber was considered in our simulations. Futurework will need to include the effects of multiplicity among outer planetary systems. Havingtwo or more Jovian planets would modify the orbital precession frequencies and increase thecomplexity of the system’ secular frequency spectrum, particularly if the giant planets aremutually interacting. This would have consequences for the dynamics and stability of STIPsand could increase the fraction of inner systems that become unstable. As an example, considerstudies based, at least in part, on the Solar System, which show that the stability of the innerSolar System or inner system analogues can depend sensitively on the orbital configuration ofthe Solar System’s outer planets (e.g., Lithwick and Wu, 2011; Agnor and Lin, 2012; Clementand Kaib, 2017). The lack of planets detected with orbital periods between 1000 and 10000days its mainly due to detection biases. A deeper understanding of the gravitational effects ofunseen outer planetary systems on configurations like the detected STIPs might help to improvethe statistics in this period range. Further research is needed.5.3 STIPs and giant planetsMigration and in-situ formation are not necessarily exclusive mechanisms. There are manyuncertainties in both situations, but there are also physical reasons for thinking that either onecould operate, at least under some situations. Both should be explored.The in-situ formation paradigm could give rise to a rich variety of planetary types on hot,warm or cold orbital distances. The details of their formation may be strongly related tothe disk density and the corresponding gas-to-solid fraction. One important consideration forthe formation models is the fraction of massive, rocky planets with and without a gaseousenvelope, since many planets are close to the typical critical mass (Lee et al., 2014). In our in-situ paradigm, rocky, massive planets would be those that consolidated late, when the gaseousdisk had already dissipated.A present concern is how to determine whether STIPs are just the inner region of a morecomplex system, or if these configurations are complete systems. One option is to determinethe behavior of the inner system with the aid of N-body simulations, and analyze the possibleimplications in the dynamics and subsequent observation, as we did in Chapter 4. Otheroptions would involve simulating the simultaneous formation of STIPs and an outer systemof giant planets. The last option is computationally more difficult to realize because it wouldneed hydrodynamical simulations and depends on many unknown parameters for the initialconditions.100Chapter 6ConclusionsIn Chapter 3, we tested the hypothesis that HJ formation is an extreme outcome of metastabilityin a STIP. For this study, we have run simulations of STIPs and explored the consolidation ofcritical cores in both gas-free and gas-rich environments. The gas-rich simulations are N-bodyrealizations with a prescription for gaseous tidal damping. We find that the consolidation ofcritical cores is only possible in the gas-free environment, because the gaseous tidal dampingefficiently removes any source of instability. However, as discussed in Chapter 5, there are manyother considerations that need to be explored before we can say that STIP instability is notpossible in the presence of gas. The systems that resulted from evolution with gaseous tidaldamping were long-term stable. However, instability would still occur if their eccentricities wereperturbed to values ∼ 0.01.The period-ratio distribution of the damped systems is independent of the initial orbital distri-bution and has a characteristic shape and behavior. The period ratios avoid first-order MMRs,especially the 2:1 resonance. Although the tendency is to avoid i:i+ 1 resonances, the distribu-tion presents overdensities at slightly larger values. The two highest peaks in the distributionare located at 1.57 and 1.38, just beyond the 3:2 and 4:3 ratios.Among the N-body simulations of STIPs with gaseous tidal damping, we found two realizationsin which planets evolved into coorbital configurations. One of these remains in the 1:1 MMRthroughout the length of the simulation and is stable for at least 10 Myr after the gas wasremoved. Following up this discovery, we examined the detection of the pair with the transitmethod by synthesizing their lightcurve at 3 different epochs, mimicking Kepler observations,with the three epochs reflecting the different regions of the coorbital interaction. The transitseparations between the primary and secondary were found to change linearly, non-linearly, orremain roughly constant with time, depending on the epoch. The synthetic transit lightcurveswere then analyzed with the Kepler Transit Planet Search pipeline (TPS). Only the biggerplanet can be detected without modifying or constraining the searching parameters of thepipeline. The detection of the smaller planet (about 10 times smaller) was only possible on thelinear transit lightcurve and if the pipeline was constrain to search for periods within 10% theperiod of the larger planet. Our preliminary tests of the transit method as a viable techniquefor detecting coorbital planets seems promising. More testing is needed, however. If successful,transit observations could provide an initial determination of the total mass of the coorbital101pair and the mass ratio between them, not only of their sizes and periods.In Chapter 4, we studied the stability and observability of two known high-multiplicity STIPs,K11 and K90, in the presence of an outer Jupiter-like planet, through both N-body simulationsand secular theory. The presence of the perturber causes dynamical rigidity about a commonorbital plane for the inner planets, while the stability of the system remains unaltered whencompared with unperturbed realizations for most perturber locations. The observed instabilityseems to be inherent to STIPs, suggesting secular resonances among the planets. The rigidbehavior of the orbital plane occurred for most of the parameter space that we explored as longas no instability developed. The presence of the perturber also caused two possible effects onsystems that are otherwise stable: (1) the orbital plane of the planets could be separated into twodistinct planes, as in K90+; and (2) the system could become unstable for particular perturberlocations. The N-body simulations and secular analysis demonstrate that the instability andmultiple orbital planes are consequences of the eccentricity and inclination secular resonances,respectively. For K11, we suggest that the eccentricity resonance close to K11b is the source ofthe system’s inherent instability.Comparing STIP secular eigenfrequencies to the synthetic counterparts provides a deeper insightinto the coupling and possible presence of MMRs between planets. K11’s nodal precessionfrequencies indicate dynamical rigidity, seemingly due to K11g, although this is dependent onthe actual mass of K11g. In K90, our simulations show a 5:4 MMR between K90b and K90c,as well as a 3:2 MMR between K90g and K90h.Observations of the rigid behavior of a STIP would indicate the existence of an outer planetarysystem. It is possible that some of the detected planetary systems with low multiplicity arepart of higher multiplicity STIPs that are affected by secular resonances.102Chapter 7Future workHere, I briefly describe a number of follow-up studies that could be carried out. Each iscomplementary to the work in this dissertation, and is a natural extension of many of theprojects. The reasoning behind each point is described throughout Chapter 5.1. The total planetary mass of the synthetic planetary systems in Sects. 3.1.2 and 3.3, wason average 15 M⊕, which in the best of circumstances could end up consolidating onecritical core. To further evaluate the feasibility of the in-situ formation of gas short-period giant (gSPG) a higher average planetary per system should be explored. As aconsequence, follow-up of the present work we will need to produce synthetic STIPs witha total mass of at least 30 M⊕ and run simulations with the same methodology as thatused in Ch. 3. Finally, a comparison with the present work would indicate whether thein-situ formation paradigm is feasible.2. Due to the unknown parameters that realistically describe protoplanetary disks, like typ-ical gas surface densities or temperature profiles, different scalings and radial profiles ofτlocal (Eq. 2.5) could be explored in the gaseous tidal damping experiments. I am particu-larly interested in testing how variations in τlocal shape the final period-ratio distribution.3. As mentioned before, the gaseous tidal damping is one of several planet-disk interac-tions that are known to affect the dynamics of planetary systems. I plan to includethe parametrization of additional planet-disk interactions within my IAS15 code. Im-mediately, I am interested in interactions that excite the eccentricity and inclination ofthe planets (e.g., Ben´ıtez-Llambay et al., 2015; Eklund and Masset, 2017). After theparametrizations are included in IAS15, the experiment exploring the formation of criti-cal cores would be repeated and analyzed accordingly.4. Although, the longterm stability of planets in horseshoe configuration has been previouslyevaluated (Laughlin and Chambers, 2002; C´uk et al., 2012). I found an inconsistency withthe published threshold after running preliminary simulations to confirm the validity ofthat threshold. I seek to properly verify the longterm stability for horseshoe configurationswith detailed and systematic N-body simulations.1035. Section 3.4 demonstrated that the detection of horseshoe planets with the transit methodseems promising, but that more research is needed to properly characterize the transitlightcurves. As a logical continuation, I plan to expand the research on coorbital transitlightcurves and detection by synthesizing and analyzing lightcurves produced by differentplanetary sizes, masses and semi-major axes.6. As part of the horseshoe parameter characterization, it is necessary to determine thecorrelation of the horseshoe period with planetary masses and semi-major axes.From Chapter 3, Sect. 3.3, has not yet been published which could be done as a single papertaking into account the lack of critical core consolidation, period distribution and evaluation ofcoorbital transit lightcurves. An alternative is to publish the work in Sect. 3.3 with additionalcomplementary work in two papers. One could focus on the critical core consolidation in thepresence of gas, with its respective period distribution, including an additional spacing setupamong planets to verify the behavior of the period ratio with gaseous tidal damping. Thesecond could be about transit lightcurves of horseshoe planets, and verification of the longtermstability threshold. In the latter paper, I could also include the transit lightcurve of a moonfor comparison and evaluate whether Kepler-1625 bI (Teachey et al., 2018) is likely to be acoorbital planet or an exomoon.104BibliographyAbramowitz, M. and Stegun, I. A.: 1964, Handbook of mathematical functions with formulas,graphs, and mathematical tables, U.S. Govt. Print. Off.Agnor, C. B. and Lin, D. N. C.: 2012, Astrophys. J. 745, 143Agol, E., Steffen, J., Sari, R., and Clarkson, W.: 2005, MNRAS 359, 567Alibert, Y., Mordasini, C., Benz, W., and Winisdoerffer, C.: 2005, Astron. Astrophys. 434,343Armitage, P. J.: 2007, ArXiv e-printsBailes, M., Bates, S. D., Bhalerao, V., Bhat, N. D. R., Burgay, M., Burke-Spolaor, S., D’Amico,N., Johnston, S., Keith, M. J., Kramer, M., Kulkarni, S. R., Levin, L., Lyne, A. G., Milia,S., Possenti, A., Spitler, L., Stappers, B., and van Straten, W.: 2011, Science 333, 1717Ballard, S. and Johnson, J. A.: 2016, Astrophys. J. 816, 66Batalha, N. M., Rowe, J. F., Bryson, S. T., Barclay, T., Burke, C. J., Caldwell, D. A., Chris-tiansen, J. L., Mullally, F., Thompson, S. E., Brown, T. M., Dupree, A. K., Fabrycky, D. C.,Ford, E. B., Fortney, J. J., Gilliland, R. L., Isaacson, H., Latham, D. W., Marcy, G. W.,Quinn, S. N., Ragozzine, D., Shporer, A., Borucki, W. J., Ciardi, D. R., Gautier, Thomas N.,I., Haas, M. R., Jenkins, J. M., Koch, D. G., Lissauer, J. J., Rapin, W., Basri, G. S., Boss,A. P., Buchhave, L. A., Carter, J. A., Charbonneau, D., Christensen-Dalsgaard, J., Clarke,B. D., Cochran, W. D., Demory, B.-O., Desert, J.-M., Devore, E., Doyle, L. R., Esquerdo,G. A., Everett, M., Fressin, F., Geary, J. C., Girouard, F. R., Gould, A., Hall, J. R., Holman,M. J., Howard, A. W., Howell, S. B., Ibrahim, K. A., Kinemuchi, K., Kjeldsen, H., Klaus,T. C., Li, J., Lucas, P. W., Meibom, S., Morris, R. L., Prsˇa, A., Quintana, E., Sanderfer,D. T., Sasselov, D., Seader, S. E., Smith, J. C., Steffen, J. H., Still, M., Stumpe, M. C., Tarter,J. C., Tenenbaum, P., Torres, G., Twicken, J. D., Uddin, K., Van Cleve, J., Walkowicz, L.,and Welsh, W. F.: 2013, Astrophys. J., Suppl. Ser. 204, 24Bate, M. R., Lodato, G., and Pringle, J. E.: 2010, MNRAS 401, 1505Batygin, K.: 2012, Nature 491, 418105Batygin, K., Morbidelli, A., and Tsiganis, K.: 2011, Astron. Astrophys. 533, A7Beauge´, C. and Nesvorny´, D.: 2012, Astrophys. J. 751, 119Beauge´, C., Sa´ndor, Z., E´rdi, B., and Su¨li, A´.: 2007, Astron. Astrophys. 463, 359Becker, J. C., Vanderburg, A., Adams, F. C., Rappaport, S. A., and Schwengeler, H. M.: 2015,Astrophys. J. 812, L18Ben´ıtez-Llambay, P., Masset, F., Koenigsberger, G., and Szula´gyi, J.: 2015, Nature 520, 63Boley, A. C.: 2009, Astrophys. J. 695, L53Boley, A. C., Granados Contreras, A. P., and Gladman, B.: 2016, Astrophys. J., Lett. 817,L17Boley, A. C., Morris, M. A., and Ford, E. B.: 2014, Astrophys. J., Lett. 792, L27Borucki, W. J., Agol, E., Fressin, F., Kaltenegger, L., Rowe, J., Isaacson, H., Fischer, D.,Batalha, N., Lissauer, J. J., Marcy, G. W., Fabrycky, D., De´sert, J.-M., Bryson, S. T.,Barclay, T., Bastien, F., Boss, A., Brugamyer, E., Buchhave, L. A., Burke, C., Caldwell,D. A., Carter, J., Charbonneau, D., Crepp, J. R., Christensen-Dalsgaard, J., Christiansen,J. L., Ciardi, D., Cochran, W. D., DeVore, E., Doyle, L., Dupree, A. K., Endl, M., Everett,M. E., Ford, E. B., Fortney, J., Gautier, T. N., Geary, J. C., Gould, A., Haas, M., Henze,C., Howard, A. W., Howell, S. B., Huber, D., Jenkins, J. M., Kjeldsen, H., Kolbl, R.,Kolodziejczak, J., Latham, D. W., Lee, B. L., Lopez, E., Mullally, F., Orosz, J. A., Prsa,A., Quintana, E. V., Sanchis-Ojeda, R., Sasselov, D., Seader, S., Shporer, A., Steffen, J. H.,Still, M., Tenenbaum, P., Thompson, S. E., Torres, G., Twicken, J. D., Welsh, W. F., andWinn, J. N.: 2013, Science 340, 587Boue´, G. and Fabrycky, D. C.: 2014a, Astrophys. J. 789, 110Boue´, G. and Fabrycky, D. C.: 2014b, Astrophys. J. 789, 111Brouwer, D. and Clemence, G. M.: 1961, Methods of celestial mechanics, Academic Press, NewYorkBrouwer, D. and van Woerkom, A. J. J.: 1950, in Astronomical Papers American Ephemeris,Vol. 13, pp 81–107Bryan, M. L., Knutson, H. A., Howard, A. W., Ngo, H., Batygin, K., Crepp, J. R., Fulton,B. J., Hinkley, S., Isaacson, H., Johnson, J. A., Marcy, G. W., and Wright, J. T.: 2016,Astrophys. J. 821, 89Burke, C. J., Bryson, S. T., Mullally, F., Rowe, J. F., Christiansen, J. L., Thompson, S. E.,Coughlin, J. L., Haas, M. R., Batalha, N. M., Caldwell, D. A., Jenkins, J. M., Still, M.,106Barclay, T., Borucki, W. J., Chaplin, W. J., Ciardi, D. R., Clarke, B. D., Cochran, W. D.,Demory, B.-O., Esquerdo, G. A., Gautier, III, T. N., Gilliland, R. L., Girouard, F. R.,Havel, M., Henze, C. E., Howell, S. B., Huber, D., Latham, D. W., Li, J., Morehead, R. C.,Morton, T. D., Pepper, J., Quintana, E., Ragozzine, D., Seader, S. E., Shah, Y., Shporer,A., Tenenbaum, P., Twicken, J. D., and Wolfgang, A.: 2014, Astrophys. J., Suppl. Ser. 210,19Burrows, A., Hubbard, W. B., Lunine, J. I., and Liebert, J.: 2001, Reviews of Modern Physics73, 719Cabrera, J., Csizmadia, S., Lehmann, H., Dvorak, R., Gandolfi, D., Rauer, H., Erikson, A.,Dreyer, C., Eigmu¨ller, P., and Hatzes, A.: 2014, Astrophys. J. 781, 18Chambers, J. E.: 1999, MNRAS 304, 793Chambers, J. E. and Wetherill, G. W.: 1998, Icarus 136, 304Chatterjee, S. and Ford, E. B.: 2015, Astrophys. J. 803, 33Chatterjee, S., Ford, E. B., Matsumura, S., and Rasio, F. A.: 2008, Astrophys. J. 686, 580Chen, J. and Kipping, D.: 2017, Astrophys. J. 834, 17Chiang, E. and Laughlin, G.: 2013, MNRAS 431, 3444Chiang, E. and Youdin, A. N.: 2010, Annual Review of Earth and Planetary Sciences 38, 493Clement, M. S. and Kaib, N. A.: 2017, Icarus 288, 88Cresswell, P., Dirksen, G., Kley, W., and Nelson, R. P.: 2007, Astron. Astrophys. 473, 329Cresswell, P. and Nelson, R. P.: 2006, Astron. Astrophys. 450, 833Cresswell, P. and Nelson, R. P.: 2008, Astron. Astrophys. 482, 677Crida, A. and Batygin, K.: 2014, Astron. Astrophys. 567, A42C´uk, M., Hamilton, D. P., and Holman, M. J.: 2012, MNRAS 426, 3051Cumming, A.: 2004, MNRAS 354, 1165Cumming, A., Butler, R. P., Marcy, G. W., Vogt, S. S., Wright, J. T., and Fischer, D. A.:2008, Publications of the Astronomical Society of the Pacific 120, 531Dawson, R. I., Chiang, E., and Lee, E. J.: 2015, MNRAS 453, 1471Dawson, R. I. and Murray-Clay, R. A.: 2013, Astrophys. J., Lett. 767, L24107Dones, H., Zahnle, K. J., and Alvarellos, J. L.: 2018, in American Astronomical Society, DDAmeeting #49, id. #102.02Dong, S. and Zhu, Z.: 2013, Astrophys. J. 778, 53Eklund, H. and Masset, F. S.: 2017, MNRAS 469, 206Everhart, E.: 1985, in A. Carusi and G. B. Valsecchi (eds.), Dynamics of Comets: Their Originand Evolution, Vol. 115 of Proceedings of IAU, p. 185, Dordrecht: Reidel, Astrophysics andSpace Science LibraryFabrycky, D. and Tremaine, S.: 2007, Astrophys. J. 669, 1298Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., Rowe, J. F., Steffen, J. H., Agol, E., Barclay, T.,Batalha, N., Borucki, W., Ciardi, D. R., Ford, E. B., Gautier, T. N., Geary, J. C., Holman,M. J., Jenkins, J. M., Li, J., Morehead, R. C., Morris, R. L., Shporer, A., Smith, J. C., Still,M., and Van Cleve, J.: 2014, Astrophys. J. 790, 146Fabrycky, D. C. and Winn, J. N.: 2009, Astrophys. J. 696, 1230Fischer, D. A. and Valenti, J.: 2005, Astrophys. J. 622, 1102Ford, E. B.: 2014, Proceedings of the National Academy of Science 111, 12616Ford, E. B. and Gaudi, B. S.: 2006, Astrophys. J., Lett. 652, L137Ford, E. B. and Holman, M. J.: 2007, Astrophys. J., Lett. 664, L51Fressin, F., Torres, G., Charbonneau, D., Bryson, S. T., Christiansen, J., Dressing, C. D.,Jenkins, J. M., Walkowicz, L. M., and Batalha, N. M.: 2013, Astrophys. J. 766, 81Fulton, B. J., Petigura, E. A., Howard, A. W., Isaacson, H., Marcy, G. W., Cargile, P. A.,Hebb, L., Weiss, L. M., Johnson, J. A., Morton, T. D., Sinukoff, E., Crossfield, I. J. M., andHirsch, L. A.: 2017, Astron. J. 154, 109Gaudi, B. S., Seager, S., and Mallen-Ornelas, G.: 2005, Astrophys. J. 623, 472Gilliland, R. L., Chaplin, W. J., Dunham, E. W., Argabright, V. S., Borucki, W. J., Basri,G., Bryson, S. T., Buzasi, D. L., Caldwell, D. A., Elsworth, Y. P., Jenkins, J. M., Koch,D. G., Kolodziejczak, J., Miglio, A., van Cleve, J., Walkowicz, L. M., and Welsh, W. F.:2011, Astrophys. J., Suppl. Ser. 197, 6Gillon, M., Jehin, E., Lederer, S. M., Delrez, L., de Wit, J., Burdanov, A., Van Grootel, V.,Burgasser, A. J., Triaud, A. H. M. J., Opitom, C., Demory, B.-O., Sahu, D. K., BardalezGagliuffi, D., Magain, P., and Queloz, D.: 2016, Nature 533, 221Ginzburg, S. and Sari, R.: 2016, Astrophys. J. 819, 116108Gladman, B.: 1993, Icarus 106, 247Goldreich, P., Lithwick, Y., and Sari, R.: 2004, Ann. Rev. of Astron. Astrophys. 42, 549Goldreich, P. and Schlichting, H. E.: 2014, Astron. J. 147, 32Goldreich, P. and Tremaine, S.: 1980, Astrophys. J. 241, 425Gonzalez, G.: 1997, MNRAS 285, 403Granados Contreras, A. P. and Boley, A. C.: 2018, Astron. J. 155, 139Gu¨ttler, C., Krause, M., Geretshauser, R. J., Speith, R., and Blum, J.: 2009, Astrophys. J.701, 130Haghighipour, N., Capen, S., and Hinse, T. C.: 2013, Celestial Mechanics and DynamicalAstronomy 117, 75Han, E., Wang, S. X., Wright, J. T., Feng, Y. K., Zhao, M., Fakhouri, O., Brown, J. I., andHancock, C.: 2014, Publ. Astron. Soc. Pac. 126, 827Hands, T. O., Alexander, R. D., and Dehnen, W.: 2014, MNRAS 445, 749Hansen, B. M. S.: 2017, MNRAS 467, 1531Hansen, B. M. S. and Murray, N.: 2012, Astrophys. J. 751, 158Hansen, B. M. S. and Murray, N.: 2013, Astrophys. J. 775, 53Hatzes, A. P. and Rauer, H.: 2015, Astrophys. J. 810, L25Hayashi, C.: 1981, Progress of Theoretical Physics Supplement 70, 35Helled, R., Bodenheimer, P., Podolak, M., Boley, A., Meru, F., Nayakshin, S., Fortney, J. J.,Mayer, L., Alibert, Y., and Boss, A. P.: 2014, in Protostars and Planets VI, p. 643Heller, R.: 2018, Astron. Astrophys. 610, A39Hildebrand, F. B.: 1987a, Introduction to numerical analysis, Dover Publications, New York,2nd editionHildebrand, F. B.: 1987b, Introduction to numerical analysis, Chapt. 8.11, Dover Publications,New York, 2nd editionHippke, M. and Angerhausen, D.: 2015, Astrophys. J. 811109Howard, A. W., Marcy, G. W., Bryson, S. T., Jenkins, J. M., Rowe, J. F., Batalha, N. M.,Borucki, W. J., Koch, D. G., Dunham, E. W., Gautier, III, T. N., Van Cleve, J., Cochran,W. D., Latham, D. W., Lissauer, J. J., Torres, G., Brown, T. M., Gilliland, R. L., Buchhave,L. A., Caldwell, D. A., Christensen-Dalsgaard, J., Ciardi, D., Fressin, F., Haas, M. R.,Howell, S. B., Kjeldsen, H., Seager, S., Rogers, L., Sasselov, D. D., Steffen, J. H., Basri,G. S., Charbonneau, D., Christiansen, J., Clarke, B., Dupree, A., Fabrycky, D. C., Fischer,D. A., Ford, E. B., Fortney, J. J., Tarter, J., Girouard, F. R., Holman, M. J., Johnson, J. A.,Klaus, T. C., Machalek, P., Moorhead, A. V., Morehead, R. C., Ragozzine, D., Tenenbaum,P., Twicken, J. D., Quinn, S. N., Isaacson, H., Shporer, A., Lucas, P. W., Walkowicz, L. M.,Welsh, W. F., Boss, A., Devore, E., Gould, A., Smith, J. C., Morris, R. L., Prsa, A., Morton,T. D., Still, M., Thompson, S. E., Mullally, F., Endl, M., and MacQueen, P. J.: 2012,Astrophys. J., Suppl. Ser. 201, 15Howard, A. W., Marcy, G. W., Johnson, J. A., Fischer, D. A., Wright, J. T., Isaacson, H.,Valenti, J. A., Anderson, J., Lin, D. N. C., and Ida, S.: 2010, Science 330, 653Huang, C., Wu, Y., and Triaud, A. H. M. J.: 2016, Astrophys. J. 825, 98Ida, S. and Lin, D. N. C.: 2004a, Astrophys. J. 604, 388Ida, S. and Lin, D. N. C.: 2004b, Astrophys. J. 616, 567Ikoma, M., Nakazawa, K., and Emori, H.: 2000, Astrophys. J. 537, 1013Janson, M.: 2013, Astrophys. J. 774Johansen, A., Davies, M. B., Church, R. P., and Holmelin, V.: 2012, Astrophys. J. 758, 39Kahan, W.: 1965, Communications of the ACM 8(1), 40Kaib, N. A., Raymond, S. N., and Duncan, M. J.: 2011, Astrophys. J., Lett. 742, L24Kenyon, S. J. and Bromley, B. C.: 2006, Astron. J. 131, 1837Kokubo, E., Kominami, J., and Ida, S.: 2006, Astrophys. J. 642, 1131Kozai, Y.: 1962, Astron. J. 67, 591Lai, D.: 2014, MNRAS 440, 3532Laskar, J.: 1985, Astron. Astrophys. 144, 133Laskar, J.: 1986, Astron. Astrophys. 157, 59Laskar, J.: 1994, Astron. Astrophys. 287, L9110Latham, D. W., Rowe, J. F., Quinn, S. N., Batalha, N. M., Borucki, W. J., Brown, T. M.,Bryson, S. T., Buchhave, L. A., Caldwell, D. A., Carter, J. A., Christiansen, J. L., Ciardi,D. R., Cochran, W. D., Dunham, E. W., Fabrycky, D. C., Ford, E. B., Gautier, Thomas N.,I., Gilliland, R. L., Holman, M. J., Howell, S. B., Ibrahim, K. A., Isaacson, H., Jenkins, J. M.,Koch, D. G., Lissauer, J. J., Marcy, G. W., Quintana, E. V., Ragozzine, D., Sasselov, D.,Shporer, A., Steffen, J. H., Welsh, W. F., and Wohler, B.: 2011, Astrophys. J. 732, L24Laughlin, G. and Chambers, J. E.: 2002, Astron. J. 124, 592Lee, E. J., Chiang, E., and Ormel, C. W.: 2014, Astrophys. J. 797, 95Lee, M. H. and Peale, S. J.: 2002, Astrophys. J. 567, 596Leleu, A., Robutel, P., and Correia, A. C. M.: 2015, Astron. Astrophys. 581Lidov, M. L.: 1962, Planetary and Space Science 9, 719Lillo-Box, J., Barrado, D., Figueira, P., Leleu, A., Santos, N. C., Correia, A. C. M., Robutel,P., and Faria, J. P.: 2018, Astron. Astrophys. 609Lin, D. N. C., Bodenheimer, P., and Richardson, D. C.: 1996, Nature 380, 606Lin, D. N. C. and Papaloizou, J.: 1986, Astrophys. J. 309, 846Lissauer, J.: 2018, in American Astronomical Society, DDA meeting #49, id. #300.04Lissauer, J. J.: 1993, Ann. Rev. of Astron. Astrophys. 31, 129Lissauer, J. J., Dawson, R. I., and Tremaine, S.: 2014a, Nature 513, 336Lissauer, J. J., Fabrycky, D. C., Ford, E. B., Borucki, W. J., Fressin, F., Marcy, G. W., Orosz,J. A., Rowe, J. F., Torres, G., Welsh, W. F., Batalha, N. M., Bryson, S. T., Buchhave,L. A., Caldwell, D. A., Carter, J. A., Charbonneau, D., Christiansen, J. L., Cochran, W. D.,Desert, J.-M., Dunham, E. W., Fanelli, M. N., Fortney, J. J., Gautier, III, T. N., Geary,J. C., Gilliland, R. L., Haas, M. R., Hall, J. R., Holman, M. J., Koch, D. G., Latham,D. W., Lopez, E., McCauliff, S., Miller, N., Morehead, R. C., Quintana, E. V., Ragozzine,D., Sasselov, D., Short, D. R., and Steffen, J. H.: 2011a, Nature 470, 53Lissauer, J. J., Jontof-Hutter, D., Rowe, J. F., Fabrycky, D. C., Lopez, E. D., Agol, E., Marcy,G. W., Deck, K. M., Fischer, D. A., Fortney, J. J., Howell, S. B., Isaacson, H., Jenkins, J. M.,Kolbl, R., Sasselov, D., Short, D. R., and Welsh, W. F.: 2013, Astrophys. J. 770, 131Lissauer, J. J., Marcy, G. W., Bryson, S. T., Rowe, J. F., Jontof-Hutter, D., Agol, E., Borucki,W. J., Carter, J. A., Ford, E. B., Gilliland, R. L., Kolbl, R., Star, K. M., Steffen, J. H., andTorres, G.: 2014b, Astrophys. J. 784, 44111Lissauer, J. J., Marcy, G. W., Rowe, J. F., Bryson, S. T., Adams, E., Buchhave, L. A., Ciardi,D. R., Cochran, W. D., Fabrycky, D. C., Ford, E. B., Fressin, F., Geary, J., Gilliland, R. L.,Holman, M. J., Howell, S. B., Jenkins, J. M., Kinemuchi, K., Koch, D. G., Morehead, R. C.,Ragozzine, D., Seader, S. E., Tanenbaum, P. G., Torres, G., and Twicken, J. D.: 2012,Astrophys. J. 750, 112Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., Steffen, J. H., Ford, E. B., Jenkins, J. M.,Shporer, A., Holman, M. J., Rowe, J. F., Quintana, E. V., Batalha, N. M., Borucki, W. J.,Bryson, S. T., Caldwell, D. A., Carter, J. A., Ciardi, D., Dunham, E. W., Fortney, J. J.,Gautier, III, T. N., Howell, S. B., Koch, D. G., Latham, D. W., Marcy, G. W., Morehead,R. C., and Sasselov, D.: 2011b, Astrophys. J., Suppl. Ser. 197, 8Lithwick, Y. and Wu, Y.: 2011, Astrophys. J. 739, 31Lopez, E. D. and Fortney, J. J.: 2014, Astrophys. J. 792, 1Madhusudhan, N. and Winn, J. N.: 2009, Astrophys. J. 693, 784Mahajan, N. and Wu, Y.: 2014, Astrophys. J. 795, 32Malhotra, R., Fox, K., Murray, C. D., and Nicholson, P. D.: 1989, Astron. Astrophys. 221, 348Mamajek, E. E.: 2009, in T. Usuda, M. Tamura, and M. Ishii (eds.), American Institute ofPhysics Conference Series, Vol. 1158 of American Institute of Physics Conference Series, pp3–10Marzari, F.: 2014, MNRAS 442, 1110Masset, F. S., D’Angelo, G., and Kley, W.: 2006, Astrophys. J. 652, 730Mayor, M., Marmier, M., Lovis, C., Udry, S., Se´gransan, D., Pepe, F., Benz, W., Bertaux,J. ., Bouchy, F., Dumusque, X., Lo Curto, G., Mordasini, C., Queloz, D., and Santos, N. C.:2011, ArXiv e-printsMayor, M. and Queloz, D.: 1995, Nature 378, 355Mizuno, H.: 1980, Progress of Theoretical Physics 64, 544Moorhead, A. V. and Adams, F. C.: 2005, Icarus 178, 517Morton, T. D., Bryson, S. T., Coughlin, J. L., Rowe, J. F., Ravichandran, G., Petigura, E. A.,Haas, M. R., and Batalha, N. M.: 2016, Astrophys. J. 822, 86Murray, C. D. and Dermott, S. F.: 2000, Solar System Dynamics, Cambridge University PressNASA Exoplanet Archive: 2018, Exoplanet and Candidate Statistics, Online112Nobili, A. and Roxburgh, I. W.: 1986, in J. Kovalevsky and V. A. Brumberg (eds.), Relativity inCelestial Mechanics and Astrometry. High Precision Dynamical Theories and ObservationalVerifications, Vol. 114 of IAU Symposium, pp 105–110Obertas, A., Van Laerhoven, C., and Tamayo, D.: 2017, Icarus 293, 52Ogilvie, G. I.: 2014, Ann. Rev. of Astron. Astrophys. 52, 171Paardekooper, S. J., Baruteau, C., Crida, A., and Kley, W.: 2010, MNRAS 401, 1950Paardekooper, S. J., Baruteau, C., and Kley, W.: 2011, MNRAS 410, 293Paardekooper, S. J. and Mellema, G.: 2006, Astron. Astrophys. 459, L17Papaloizou, J. C. B. and Larwood, J. D.: 2000, MNRAS 315, 823Perryman, M. A.: 2011, The exoplanet handbook, Cambridge University Press, Cambridge,New YorkPetigura, E. A., Marcy, G. W., and Howard, A. W.: 2013, Astrophys. J. 770, 69Petigura, E. A., Marcy, G. W., Winn, J. N., Weiss, L. M., Fulton, B. J., Howard, A. W.,Sinukoff, E., Isaacson, H., Morton, T. D., and Johnson, J. A.: 2018, Astron. J. 155, 89Piso, A.-M. A., Youdin, A. N., and Murray-Clay, R. A.: 2015, Astrophys. J. 800Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., and Greenzweig,Y.: 1996, Icarus 124, 62Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: 1992a, NumericalRecipes in C, Chapt. 16.4, Cambridge University PressPress, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: 1992b, NumericalRecipes in C, Cambridge University PressPu, B. and Wu, Y.: 2015, Astrophys. J. 807, 44Rafikov, R. R.: 2006, Astrophys. J. 648, 666Rafikov, R. R.: 2011, Astrophys. J. 727Rasio, F. A. and Ford, E. B.: 1996, Science 274, 954Raymond, S. N.: 2010, Formation and Evolution of Exoplanets, Chapt. Formation of TerrestrialPlanets, pp 123–144, Wiley Online LibraryRaymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., and Walsh, K. J.: 2014, Protostarsand Planets VI pp 595–618113Rein, H. and Spiegel, D. S.: 2015, MNRAS 446, 1424Rogers, L. A.: 2015, Astrophys. J. 801, 41Rogers, L. A. and Seager, S.: 2010a, Astrophys. J. 712, 974Rogers, L. A. and Seager, S.: 2010b, Astrophys. J. 716, 1208Rogers, T. M., Lin, D. N. C., and Lau, H. H. B.: 2012, Astrophys. J., Lett. 758, L6Rogers, T. M., Lin, D. N. C., McElwaine, J. N., and Lau, H. H. B.: 2013, Astrophys. J. 772,21Rowe, J. F., Bryson, S. T., Marcy, G. W., Lissauer, J. J., Jontof-Hutter, D., Mullally, F.,Gilliland, R. L., Issacson, H., Ford, E., Howell, S. B., Borucki, W. J., Haas, M., Huber, D.,Steffen, J. H., Thompson, S. E., Quintana, E., Barclay, T., Still, M., Fortney, J., Gautier,III, T. N., Hunter, R., Caldwell, D. A., Ciardi, D. R., Devore, E., Cochran, W., Jenkins, J.,Agol, E., Carter, J. A., and Geary, J.: 2014, Astrophys. J. 784, 45Santerne, A., Moutou, C., Tsantaki, M., Bouchy, F., He´brard, G., Adibekyan, V., Almenara,J. M., Amard, L., Barros, S. C. C., Boisse, I., Bonomo, A. S., Bruno, G., Courcol, B., Deleuil,M., Demangeon, O., Dı´az, R. F., Guillot, T., Havel, M., Montagnier, G., Rajpurohit, A. S.,Rey, J., and Santos, N. C.: 2016, Astron. Astrophys. 587Santos, N. C., Israelian, G., and Mayor, M.: 2004, Astron. Astrophys. 415, 1153Schmitt, J. R., Wang, J., Fischer, D. A., Jek, K. J., Moriarty, J. C., Boyajian, T. S., Schwamb,M. E., Lintott, C., Lynn, S., Smith, A. M., Parrish, M., Schawinski, K., Simpson, R., La-Course, D., Omohundro, M. R., Winarski, T., Goodman, S. J., Jebson, T., Schwengeler,H. M., Paterson, D. A., Sejpka, J., Terentev, I., Jacobs, T., Alsaadi, N., Bailey, R. C., Gin-man, T., Granado, P., Vonstad Guttormsen, K., Mallia, F., Papillon, A. L., Rossi, F., andSocolovsky, M.: 2014, Astron. J. 148, 28Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., and Zolotukhin, I.: 2011, Astron. Astro-phys. 532, A79Schwarz, R., Bazso´, A´., Funk, B., and Zechner, R.: 2015, MNRAS 453, 2308Shallue, C. J. and Vanderburg, A.: 2018, Astron. J. 155, 94Sinukoff, E., Howard, A. W., Petigura, E. A., Fulton, B. J., Crossfield, I. J. M., Isaacson, H.,Gonzales, E., Crepp, J. R., Brewer, J. M., Hirsch, L., Weiss, L. M., Ciardi, D. R., Schlieder,J. E., Benneke, B., Christiansen, J. L., Dressing, C. D., Hansen, B. M. S., Knutson, H. A.,Kosiarek, M., Livingston, J. H., Greene, T. P., Rogers, L. A., and Le´pine, S.: 2017, Astron.J. 153114Smith, A. W. and Lissauer, J. J.: 2009, Icarus 201, 381Southworth, J.: 2011, MNRAS 417, 2166, Last modified: 2017/05/06Steffen, J. H. and Agol, E.: 2005, MNRAS 364, L96Steffen, J. H., Ragozzine, D., Fabrycky, D. C., Carter, J. A., Ford, E. B., Holman, M. J., Rowe,J. F., Welsh, W. F., Borucki, W. J., Boss, A. P., Ciardi, D. R., and Quinn, S. N.: 2012,Proceedings of the National Academy of Science 109, 7982Stevenson, D. J.: 1982, Planetary and Space Science 30, 755Stoer, J. and Bulirsch, R.: 1980, Introduction to Numerical Analysis, Chapt. 7.2.14, New York:Springer-VerlagTan, J. C., Chatterjee, S., Hu, X., Zhu, Z., and Mohanty, S.: 2015, ArXiv e-printsTanaka, H. and Ward, W. R.: 2004, Astrophys. J. 602, 388Teachey, A. and Kipping, D. M.: 2018, Science Advances 4(10)Teachey, A., Kipping, D. M., and Schmitt, A. R.: 2018, Astron. J. 155, 36Tremaine, S. and Dong, S.: 2012, Astron. J. 143, 94Udry, S., Mayor, M., and Santos, N. C.: 2003, Astron. Astrophys. 407, 369Valencia, D., O’Connell, R. J., and Sasselov, D.: 2006, Icarus 181, 545Vokrouhlicky´, D. and Nesvorny´, D.: 2014, Astrophys. J. 791Volk, K. and Gladman, B.: 2015, Astrophys. J., Lett. 806, L26Wadhwa, M., Amelin, Y., Davis, A. M., Lugmair, G. W., Meyer, B., Gounelle, M., and Desch,S. J.: 2007, in B. Reipurth, D. Jewitt, and K. Keil (eds.), Protostars and Planets V, p. 835Ward, W. R.: 1997, Icarus 126, 261Weidenschilling, S. J.: 1977a, MNRAS 180, 57Weidenschilling, S. J.: 1977b, Astrophys. J., Suppl. Ser. 51, 153Weidenschilling, S. J. and Marzari, F.: 1996, Nature 384, 619Weiss, L. M. and Marcy, G. W.: 2014, Astrophys. J. 783, L6Weiss, L. M., Marcy, G. W., Petigura, E. A., Fulton, B. J., Howard, A. W., Winn, J. N.,Isaacson, H. T., Morton, T. D., Hirsch, L. A., Sinukoff, E. J., Cumming, A., Hebb, L., andCargile, P. A.: 2018, Astron. J. 155, 48115Winn, J. N. and Fabrycky, D. C.: 2015, Ann. Rev. of Astron. Astrophys. 53, 409Winn, J. N., Noyes, R. W., Holman, M. J., Charbonneau, D., Ohta, Y., Taruya, A., Suto, Y.,Narita, N., Turner, E. L., Johnson, J. A., Marcy, G. W., Butler, R. P., and Vogt, S. S.: 2005,Astrophys. J. 631, 1215Wisdom, J. and Holman, M.: 1991, Astron. J. 102, 1528Wisdom, J., Holman, M., and Touma, J.: 1996, Fields Institute Communications 10, 217Wolfgang, A., Rogers, L. A., and Ford, E. B.: 2016, Astrophys. J. 825, 19Wolszczan, A. and Frail, D. A.: 1992, Nature 355, 145Wright, J. T., Fakhouri, O., Marcy, G. W., Han, E., Feng, Y., Johnson, J. A., Howard, A. W.,Fischer, D. A., Valenti, J. A., Anderson, J., and Piskunov, N.: 2011, Publ. Astron. Soc. Pac.123, 412Wright, J. T., Marcy, G. W., Howard, A. W., Johnson, J. A., Morton, T. D., and Fischer,D. A.: 2012, Astrophys. J. 753, 160Wright, J. T., Upadhyay, S., Marcy, G. W., Fischer, D. A., Ford, E. B., and Johnson, J. A.:2009, Astrophys. J. 693, 1084Wu, Y. and Lithwick, Y.: 2013, Astrophys. J. 772, 74Zapolsky, H. S. and Salpeter, E. E.: 1969, Astrophys. J. 158, 809Zsom, A. and Dullemond, C. P.: 2008, Astron. Astrophys. 489, 931116Appendix AComplementary materialA.1 Laplace-Lagrange literal expansion of the disturbingfunction.A.1.1 Direct termThe literal expansion of the direct term in disturbing function to second-order in inclinationand eccentricity, and first in masses, isRD ={12b(k)1/2 +18(e2i + e2j) (−4k2 + 2αD + α2D2) b(k)1/2− 14α(s2i + s2j) (b(k−1)3/2 + b(k+1)3/2)}cos(kλj − kλi)+{14eiej(2 + 6k + 4k2 − 2αD − α2D2) b(k+1)1/2 } cos(kλj − kλi +$j −$i)+{sisj α b(k+1)3/2}cos(kλj − kλi + Ωj − Ωi)+{12ei (−2k − αD) b(k)1/2}cos(kλj − (k − 1)λi −$i)+{12ej (−1 + 2k + αD) b(k−1)1/2}cos(kλj − (k − 1)λi −$j)+{18e2i(−5k + 4k2 − 2αD + 4kαD + α2D2) b(k)1/2} cos(kλj − (k − 2)λi − 2$i)+{14eiej(−2− 6k − 4k2 + 2αD − 4kαD − α2D2) b(k−1)1/2 }× cos(kλj − (k − 2)λi −$j −$i)+{18e2j(2− 7k + 4k2 − 2αD + 4kαD + α2D2) b(k−2)1/2 } cos(kλj − (k − 2)λi − 2$j)+{12s2i α b(k−1)3/2}cos(kλj − (k − 2)λi − 2Ωi)−{sisj α b(k−1)3/2}cos(kλj − (k − 2)λi − Ωj − Ωi)+{12s2j α b(k−1)3/2}cos(kλj − (k − 2)λi − 2Ωj). (A.1)117Here, sk = sinIk2for bodies k = i, j and inclination Ik.A.1.2 Indirect termsSimilarly, the explicit expansion of the external indirect term isRE ≈[−1 + 12(e2i + e2j ) + s2i + s2j]cos(λj − λi)− eiej cos(2λj − 2λi −$j +$i)− 2 sisj cos(λj − λi − Ωj + Ωi)− 12ei cos(λj − 2λi +$i)− 32ei cos(λj −$i)− 2 ej cos(2λj − λi −$j)− 38e2i cos(λj − 3λi + 2$i)−18e2i cos(λj + λi − 2$i)+ 3 eiej cos(2λi −$j −$i)− 18e2j cos(λj + λi − 2$j)− 278e2j cos(3λj − λi − 2$j)− s2i cos(λj + λi − 2Ωi)+ 2 sisj cos(λj + λi − Ωj − λi)− s2j cos(λj + λi − 2Ωj). (A.2)Finally, the internal indirect part isRI ≈[−1 + 12(e2i + e2j ) + s2i + s2j]cos(λj − λi)− eiej cos(2λj − 2λi −$j +$i)− 2 sisj cos(λj − λi − Ωj + Ωi)− 2 ei cos(λj − 2λi +$i) + 32ei cos(λi −$j)− 12ej cos(2λj − λi −$j)− 278e2i cos(λj − 3λi + 2$i)−18e2i cos(λj + λi − 2$i)+ 3 eiej cos(2λi −$j −$i)− 18e2j cos(λj + λi − 2$j)− 38e2j cos(3λj − λi − 2$j)− s2i cos(λj + λi − 2Ωi)+ 2 sisj cos(λj + λi − Ωj − λi)− s2j cos(λj + λi − 2Ωj) (A.3)A.2 Laplace-Lagrange solution for more than two secondarymasses.The disturbing function for more than two secondary masses, using the Laplace-Lagrange so-lution isRi = nia2i[12Aii e2i +12Bii I2i118+N∑j=1,j 6=iAij eiej cos($i −$j) +N∑j=1,j 6=iBij IiIj cos(Ωi − Ωj)], (A.4)withAii =14nimc +miN∑j=1,j 6=imj αij α¯ij b(1)3/2(αij) (A.5)Aij = −14nimc +mimj αij α¯ij b(2)3/2(αij) (A.6)Bii = −14nimc +miN∑j=1,j 6=imj αij α¯ij b(1)3/2(αij) (A.7)Bij =14nimc +mimj αij α¯ij b(1)3/2(αij) (A.8)whereαij =aj/ai internal perturber, if ai > ajai/aj external perturber, if ai < aj , (A.9)andα¯ij =1 internal perturber, if ai > ajai/aj external perturber, if ai < aj , (A.10)Using the eccentricity and inclination vector definitions, Eq. A.4 transforms toRi = nia2i[12Aii (h2i + k2i ) +12Bii (p2i + q2i )+N∑j=1,j 6=iAij (hihj + kikj) +N∑j=1,j 6=iBij (pipj + qiqj)], . (A.11)According to Brouwer and Clemence (1961) the derivatives of the eccentricity and inclinationvector in Eq. A.11 areh˙i = Aii ki +N∑j=1,j 6=iAij kj (A.12)k˙i = −Aii hi −N∑j=1,j 6=iAij hj (A.13)p˙i = Bii qi +N∑j=1,j 6=iBij qj (A.14)q˙i = −Bii pi −N∑j=1,j 6=iBij pj (A.15)119and the solution changes tohi =N∑j=1eij sin(gj t+ βj), ki =N∑j=1eij cos(gj t+ βj) (A.16)pi =N∑j=1Iij sin(fj t+ γj), qi =N∑j=1Iij cos(fj t+ γj) (A.17)The matrices A and B seen in Sect. 2.2 are still square matrices, but with N ×N dimension.The process to find the eigenfrequencies, eigenvectors, forced inclination and eccentricity is thesame that in Sect. 2.2 and in Murray and Dermott (2000). The subindex i indicates a body,and j the mode.120A.3 Orbital evolution of coorbital realizations in the corotating frame(a) r-054 (b) r-628Figure A.1: Phase-on corotational evolution of the realizations in Fig. 3.12 in the frame of the larger planet. The angular location ofthe L4 and L5 Lagrange points are indicated as green dotted-lines. The richness of the orbital evolution in function of time is evident.In panel (a), the smaller planet, initially at a higher semi-major axis, due the eccentricity and inclination damping is forced to migrateinwards and catches up with the planet in red, resulting in a 1:1 MMR capture and shifts between tadpole and horseshoe orbits. (b)Starting at a lower semi-major axis, the smaller planet moves outwards and is trapped in a tadpole orbit than then evolves into ahorseshoe.121A.4 IAS15 integration procedure.Start integrationInitial conditionsStart time loopCalculate forceson bodies fromy0 and y′0Estimate g-valuesfrom b-valuesSet all bn = 0Convergence reached?Start h-loopPredict y and y′with b-values and hnfor all bodiesCalculate forces usingnew r and v forall bodiesCalculate b-valuesfrom new value ofgnCalculate new gnat each hnEstimateb˜6Estimate convergencecriterionδ˜b6bn calculationand convergenceloopAll hn used?Adaptive timestepON?Calculate newtimestepdtneed < dtused?dt = dtneedt = t+ dtCalculate newy and y′from b-valuesB-values estimationdt = dtneedt = tfinal?Stop integrationSave t, y and y′y0 = yandy′0 = y′122Start integrationInitial conditionsStart time loopCalculate forceson bodies fromy0 and y′0Estimate g-valuesfrom b-valuesSet all bn = 0Convergence reached?Start h-loopPredict y and y′with b-values and hnfor all bodiesCalculate forces usingnew r and v forall bodiesCalculate b-valuesfrom new value ofgnCalculate new gnat each hnEstimateb˜6Estimate convergencecriterionδ˜b6bn calculationand convergenceloopAll hn used?Adaptive timestepON?Calculate newtimestepdtneed < dtused?dt = dtneedt = t+ dtCalculate newy and y′from b-valuesB-values estimationdt = dtneedt = tfinal?Stop integrationSave t, y and y′y0 = yandy′0 = y′Figure A.2: IAS15 integration procedure. The red lines represent a negative answer to the questions inside each of the decision takingrhombus.123"@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2019-02"@en ; edm:isShownAt "10.14288/1.0373613"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Astronomy"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@* ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@* ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Orbital outcomes of STIPs and consequences for hot-Jupiter formation and planet diversity"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/67781"@en .