Mass and Eccentricity Constraints in the WASP-47 Planetary System from a Simultaneous Analysis of Radial Velocities & Transit Timing Variations
Lauren M. Weiss, Katherine Deck, Evan Sinukoff, Erik A. Petigura, Eric Agol, Eve J. Lee, Juliette C. Becker, Andrew W. Howard, Howard Isaacson, Ian J. M. Crossfield, Benjamin J. Fulton, Lea Hirsch
DD RAFT VERSION M ARCH
28, 2017
Preprint typeset using L A TEX style AASTeX6 v. 1.0
NEW INSIGHTS ON PLANET FORMATION IN WASP-47 FROM A SIMULTANEOUS ANALYSIS OF RADIALVELOCITIES AND TRANSIT TIMING VARIATIONS L AUREN
M. W
EISS , K
ATHERINE
M. D
ECK , E VAN S INUKOFF , E
RIK
A. P
ETIGURA , E
RIC A GOL , E VE J. L EE ,J ULIETTE
C. B
ECKER , A NDREW
W. H
OWARD , H OWARD I SAACSON , I AN J. M. C
ROSSFIELD , B ENJAMIN
J. F
ULTON , L EA H IRSCH B JÖRN B ENNEKE Institut de Recherche sur les Exoplanètes, Université de Montréal, Montréal, QC, Canada California Institute of Technology, Pasadena, CA, USA Institute for Astronomy, University of Hawai‘i at M¯anoa, Honolulu, HI, USA Astronomy Department, University of Washington, Seattle, WA, USA NASA Astrobiology Institute Virtual Planet Laboratory, Seattle, WA, USA Astronomy Department, University of California, Berkeley, CA, USA Astronomy Department, University of Michigan, Ann Arbor, MI, USA Astronomy Department, University of California, Santa Cruz, CA, USA Trottier Fellow NSERC Fellow NASA Hubble Fellow
ABSTRACTMeasuring precise planet masses, densities, and orbital dynamics in individual planetary systems is an impor-tant pathway toward understanding planet formation. The WASP-47 system has an unusual architecture thatmotivates a complex formation theory. The system includes a hot Jupiter (“b”) neighbored by interior (“e”)and exterior (“d”) sub-Neptunes, and a long-period eccentric giant planet (“c”). We simultaneously modeledtransit times from the Kepler K2 Mission and 118 radial velocities to determine precise masses, densities, andKeplerian orbital elements of the WASP-47 planets. Combining RVs and TTVs provides a better estimate ofthe mass of planet d (13 . ± . M ⊕ ) than obtained with only RVs (12 . ± . M ⊕ ) or TTVs (16 . ± . M ⊕ ).Planets e and d have high densities for their size, consistent with a history of photo-evaporation and/or forma-tion in a volatile-poor environment. Through our RV and TTV analysis, we find that the planetary orbits haveeccentricities similar to the solar system planets. The WASP-47 system has three similarities to our own solarsystem: (1) the planetary orbits are nearly circular and coplanar, (2) the planets are not trapped in mean motionresonances, and (3) the planets have diverse compositions. None of the current single-process exoplanet for-mation theories adequately reproduce these three characteristics of the WASP-47 system (or our solar system).We propose that WASP-47, like the solar system, formed in two stages: first, the giant planets formed in agas-rich disk and migrated to their present locations, and second, the high-density sub-Neptunes formed in situin a gas-poor environment. INTRODUCTIONOne of the key questions driving exoplanet science is theformation of planetary systems in general and the solar sys-tem in particular. The
Kepler
Mission (Borucki et al. 2010;Koch et al. 2010) led to a wealth of statistical measurementsthat provide valuable insight into planet formation. Smallplanets close to their stars are a common outcome of planetformation (Howard et al. 2012; Fang & Margot 2012; Batalhaet al. 2013; Fressin et al. 2013; Petigura et al. 2013b; Dress-ing & Charbonneau 2015), such that half of sun-like starshave at least one planet smaller than Neptune within the or-bital distance of Mercury (Petigura et al. 2013a; Burke et al.2015). Many small planets close to their stars are in multi- planet systems (Lissauer et al. 2012; Fabrycky et al. 2014;Rowe et al. 2014; Lissauer et al. 2014). Given the abundanceof small, compact planetary systems around other stars, is ourown solar system unusual in that it is barren from Mercury’sorbit inward? Because the
Kepler
Mission only obtained 4years of continuous photometry and only observed 150,000stars, it had poor sensitivity to long-period planets, which areunlikely to transit. If
Kepler were pointed at our solar systemand were lucky enough to discover the inner planets, it stillmost likely would have missed the planets from Mars out.How can we reconcile planet formation theory for the close-in exoplanets with planet formation theory for our spaciouslyspread solar system? a r X i v : . [ a s t r o - ph . E P ] M a r WASP-47 is a system of unusual architecture that mightbe a Rosetta Stone for linking the exoplanet population tothe solar system. WASP-47 contains a transiting, Jupiter-sizeplanet with an orbital period of 4.2 days (a “hot Jupiter”)that was detected from the ground-based WASP-South tran-sit survey (Hellier et al. 2012, WASP-47 b). What makesWASP-47 b unusual is that, contrary to the vast majority ofhot Jupiters, which do not have nearby planetary companions(Steffen et al. 2012; Bryan et al. 2016), WASP-47 b has twonearby neighbors: an interior, transiting planet with an or-bital period of less than a day (WASP-47 e) and an exterior,transiting planet with an orbital period of 9.0 days (WASP-47 d). The system also has a distant, moderately eccentricplanet (WASP-47 c). While the compactness of the WASP-47 inner planetary system is comparable to other Kepler sys-tems, especially those that contain ultra-short period planets(Sanchis-Ojeda et al. 2013), the combination of the compactplanetary system with a hot Jupiter is unprecedented amongthe 2217 planetary systems studied to date. The architectureof WASP-47 was not predicted by planet formation theory,and so uncovering a physically plausible formation mecha-nism for WASP-47 will deepen our understanding of planetformation in general.To better understand which of the various physical mod-els of planet formation and evolution were important in theWASP-47 system, we would like to measure the masses, den-sities, bulk compositions, and orbital dynamics of all theplanets as precisely as the current data permit. The compo-sitions of the planets might provide clues about where theyformed within the proto-planetary disk–for instance, if theplanets are rich in water or other high mean-molecular weightvolatiles, they might have formed beyond a molecular snow-line. Furthermore, the present-day orbital elements for theplanets can be related to their dynamical history: the semi-major axes and eccentricities of the planets today relate tohow they have exchanged energy and angular momentum inthe past.Several analyses of this system have already character-ized various dynamical properties of the WASP-47 planets.The discovery paper (Hellier et al. 2012) used two yearsof ground-based photometry to find WASP-47 b in transit,and also obtained 19 low-precision radial velocities (RVs)to measure the planet’s mass. Becker et al. (2015, hereafterB15) discovered two additional transiting planets (e and d)in transits from K2, characterized the planet masses withtransit timing variations (TTVs), and used the Mercury N-body integrator (Chambers 1999) to explore the dynamicalstability of the planets. Sanchis-Ojeda et al. (2015) measuredthe projected spin-orbit obliquity of the hot Jupiter via theRossiter-McLaughlin effect, finding that the planetary orbital Based on a 2016-11-24 query of Exoplanets.org. axis and the stellar spin axis are not strongly misaligned. Daiet al. (2015) obtained high-cadence precision RVs of the sys-tem, precisely characterizing the mass of the giant planet andplacing new mass constraints on the other transiting planets.Neveu-VanMalle et al. (2016) discovered a long-period giantplanet with a multi-year baseline of radial velocities. Alme-nara et al. (2016) simultaneously modeled the K2 light curveand the literature RVs, arriving at planet masses that were de-termined to a precision of ∼ M p sin i measurements( < MEASUREMENTSThe measurements we use in this analysis are all availablein the literature. We combine the 108 transit times (TTs orTTVs) of WASP-47 e, b, and d (B15) with 118 measurementsof the radial velocity (RV) of the WASP-47 host star fromHellier et al. (2012), Dai et al. (2015), Neveu-VanMalle et al.(2016), and S17.To combine the RV measurements, we use the valuesfor the RV zero-point offset ( γ ) and jitter ( σ jit ) determinedin S17. The zero-point offset is added to each RV mea-surement, and the jitter is added to each RV uncertainty inquadrature. For the Hellier et al. (2012) CORALIE RVs,these values are γ = 27070 . − , σ jit = 5 . − . For theNeveu-VanMalle et al. (2016) CORALIE RVs, these valuesare γ = 27085 . − , σ jit = 6 . − . For the Dai et al.(2015) Magellan-PFS RVs, these values are γ = 20 . − , σ jit = 6 . − . For the S17 Keck-HIRES RVs, these valuesare γ = 6 . − , σ jit = 3 . − . For simplicity, we keep thevalues of the zero-point offset and jitters fixed at the valuesdetermined in S17. The zero-point offsets and jitter are statis-tical properties of the RVs, and so we do not expect the TTVsto provide any new information about these parameters. TTVS-ONLY ANALYSIS WITH N-BODY ANDANALYTIC APPROACHESIn this section, we fit the WASP-47 TTVs as measuredin B15 using two different approaches. First, we do a fullN-body simulation of the three transiting planets to repro-duce the observed TTVs using the publicly available code
TTVFast (Deck et al. 2014). Then, we use
TTVFaster (Agol & Deck 2016a,b), a publicly available code that an-alytically models orbits to first order in eccentricity. The
TTVFaster code was designed to reproduce both low-frequency sinusoidal features and high-frequency “chop-ping” patterns in the TTVs. In the tests below, we deter-mine that
TTVFaster , which is orders of magnitude fasterthan
TTVFast , is appropriate for modeling the TTVs in theWASP-47 system.3.1.
Modeling transit times with an N-body Integrator
Table 1 . Priors on Dynamical Pa-rameters
Parameter Priors
M M > 0, Hill criterion
P P > 0, Hill criterion TT None √ e cos ω e < .
06, Hill criterion √ e sin ω e < .
06, Hill criterion √ e c cos ω c e <
1, Hill criterion √ e c sin ω c e <
1, Hill criterion
We used the publicly available N-body integrator
TTVFast (Deck et al. 2014) with the python wrapper ttvfast-python to forward-model the transit times ofthe inner three planets. Unlike in the B15 analysis, we al-lowed all of the initial osculating elements, particularly theorbital periods and initial times of transits, to vary. Thus,the variables for each planet k are: the mass of the planet M k , the orbital period P k , the first time of transit tt k , and theeccentricity parametrization √ e k cos ω k , √ e k sin ω k . We lim-ited e < .
06 for the three inner planets, in accordance withthe 10 Myr stability analysis in B15. For each adjacent pairof planets, we satisfied the Hill criteria for stability for low-eccentricity orbits, as described in Equations 24 and 28 ofGladman (1993): p − q > . µ + µ ) / (1) p − q > (cid:113) / e + e ) + × max( µ , µ ) / (2)where the subscript 1 refers to the inner planet and 2 refersto the outer planet, q is the apoapse distance of the innerplanet q = 1 + e , p is the periapse distance of the outer https://github.com/mindriot101/ttvfast-python planet p = (1 − e ) a a , e is the eccentricity, and µ is theplanet-to-star mass ratio. We also required the planet massesand orbital periods to have positive values. Because theorbits of the planets are very nearly coplanar (Becker et al.2015; Almenara et al. 2016), we fixed the orbital inclinationand longitude of ascending node in a manner consistent withcoplanar, edge-on orbits. See Table 3.1 for a summary of thepriors and constraints. Table 2 . Dynamical Parameters from Best N-Body Fit to TTVs Only (TTVFast)
Parameter Median ± Std. Dev. Units M e ± M ⊕ M b ± M ⊕ M d ± M ⊕ P e ± P b ± P d ± TT e ± TT b ± TT d ± √ e e cos ω e -0.006 ± √ e e sin ω e -0.008 ± √ e b cos ω b ± √ e b sin ω b ± √ e d cos ω d -0.09 ± √ e d sin ω d ± OTE —These are the initial astrocentric Keplerian orbitalelements, reported at epoch BJD 2456979.4961. Theyare not the time-averaged orbital properties of the plan-ets.
We used the Markov Chain Monte Carlo (MCMC) Pythonpackage emcee (Foreman-Mackey et al. 2013) to explore theposteriors of various combinations of the dynamical param-eters. We ran 60 walkers 5 × steps each, throwing awaythe first 10 steps as burn-in, and checked that the multivari-ate extension of the potential scale reduction factor (PSRF)statistic ( ˆ R < to the transit times us-ing TTVFast is shown in Figure 1. Table 3.1 summarizesour results from this N-body fit to the TTVs.We find that, for planets e and b, the masses and eccen-tricities of the planets are highly covariant with the initialosculating orbital period and the initial transit time of neigh-boring planets (see Figure 2). This is because the initial oscu-lating orbital period is translated to an instantaneous velocity from the MCMC maximum likelihood Figure 1 . The observed minus linear-ephemeris calculated tran-sit times of WASP-47 e (blue points, top panel), b (green points,middle panel), and d (brown points, bottom panel). The best-fitN-body model to the TTVs alone (using
TTVFast , black con-nected dots) and the best-fit analytic model to the TTVs alone (using
TTVFaster , gray connected dots) are shown. and acceleration, and the planet’s acceleration depends onthe mass and position of its neighbor. Therefore, it is criticalto allow the initial orbital periods and times of transit of allthe planets to vary in order to explore the full range of pos-sible planet masses. Furthermore, only one super-period ofthe TTVs is observed, and so an average orbital period and arepresentative time of transit are not as well determined forplanets b and d as they might be with the observation of mul-tiple super-periods. Thus, we find that the TTVs do not con-strain the masses of WASP-47 e or WASP-47 b as tightly aswhat is reported in B15. Whereas B15 find M b = 341 + − M ⊕ ,we find M b = 549 ± M ⊕ , using the exact same TTV mea- surements. The mass for planet e reported in B15 stems fromtheir choice of prior: they find M e < M ⊕ , whereas wefind M e = 176 ± M ⊕ . However, the TTVs of planet bplace strong constraints on the mass of planet d. B15 find M d = 15 . ± M ⊕ , and we find M d = 16 . ± . M ⊕ . Whilean analytic covariance between planet masses and the freeeccentricity exists (Lithwick et al. 2012), our stability con-straint that e < .
06 minimizes the effects of this degeneracy.3.2.
Modeling transit times analytically
We used the publicly available analytic TTV package
TTVFaster to model the transit times of the inner threeplanets observed in B15. The
TTVFaster code analyticallytransforms the planet mass M k and average orbital elements P k , tt k , e k cos ω k , e k sin ω k into a TTV pattern, to first orderin eccentricity, to a user-specified precision in the disturbingfunction. We found that modeling to sixth order in the expan-sion of the Laplace coefficient b j / ( j = { , , , , , , } ),(Murray & Dermott 2000) was sufficient to reproduce the ob-served TTV signature with the same fidelity as produced inthe N-body analysis (see Figure 1). In general, TTVFaster is designed to work for planets that (1) are not extremelyclose to a mean motion resonance, (2) have low eccentrici-ties, (3) have low masses. Because WASP-47 b is a Jupiter-mass planet, we wanted to see if
TTVFaster modeled theorbital dynamics correctly.Incorporating the priors from Table 3.1, we used Pythonpackages lmfit (Newville et al. 2014) and emcee to ex-plore the posteriors of the dynamical parameters. We ran 100walkers 20,000 steps, throwing away the first 4000 steps asburn-in. We note that the chains converged much faster (ac-cording to the PSRF statistic) when we used
TTVFaster than when we used
TTVFast , because the TTVs providebetter constraints on the average orbital parameters than theydo on the initial orbital parameters.The mass and eccentricity distributions we determinedfrom the analytic solution to the observed TTVs are consis-tent with the N-body model (Figure 3). Since fitting an ana-lytic model to the TTVs is orders of magnitude faster than afull N-body analysis (especially when the long time baselinefor RVs is required), we use the analytic modeling techniquein the rest of this paper. COMBINING TTVS AND RVS WITH
TTVFASTER
We combined the python packages
TTVFaster and radvel (Fulton & Petigura in prep. ) to simultaneouslyfit the RVs and TTVs, resulting in refined masses and or-bital properties of all four known planets. To fit the RV andTTV data simultaneously, we maximized the following log- We used jump parameters √ e k cos ω k , √ e k sin ω k to avoid an eccentricitybias and speed convergence. https://github.com/California-Planet-Search/radvel Figure 2 . Posteriors of the planet masses, initial osculating orbital periods, and initial times of transit for WASP-47 e, b, and d in the
TTVFast (N-body) analysis. The planet masses are highly covariant with the initial orbital elements. Note that the orbital periods and times of transithere are initial osculating elements, not the time-averaged orbital elements. The MCMC chains shown are available as the Data behind theFigure. likelihood function while satisfying our priors:ln L = − N RV (cid:88) i ( RV obs , i − RV Kep , i ) σ , i − N pl (cid:88) k N TT , k (cid:88) j ( T T obs , k , j − T T model , k , j ) σ , k , j (3)where N RV is the number of RVs; N pl is the number of plan-ets; N TT , k is the number of transit times for planet k ; RV obs , i isthe i th observed RV, including the instrument-specific fixed zero-point offset γ determined in S17; RV Kep , i is the i thKeplerian-modeled RV; σ RV , i is the uncertainty in RV obs , i , in-cluding a constant jitter for each spectrometer determined inS17; T T obs , k , j is the j th observed transit time for planet k ; T T model , k , j is the j th modeled transit time for planet k ; and σ TT , k , j is the uncertainty in T T obs , k , j .Our model included all four known planets. The variableparameters for each planet k are: the mass of the planet M k ,the orbital period of the planet P k , a representative time oftransit tt k , and the eccentricity parametrization √ e k cos ω k , Figure 3 . The posteriors of the analyses of the TTVs alone. Blue: from the
TTVFast
N-body integrator MCMC; red: from the analytic TTVmodeler TTVFaster MCMC. In both analyses, we limited e < .
06 for the three inner planets, in accordance with the stability analysis fromB15. We also required Hill stability. We used jump parameters of the form √ e cos ω to avoid high eccentricity biases. The mass and eccentricityposteriors obtained from the analytic TTV analysis are in good agreement with those obtained with N-body modeling. The values above thehistograms correspond to the N-body posterior median and 1 σ bounds. The MCMC chains shown are available as the Data behind the Figure √ e k sin ω k . Note that the argument of periapse passage, ω k ,is for the planet, not the star. This formulation is consistentwith both the TTVfaster definition and the definition inMurray & Correia (2010) and Lovis & Fischer (2010).These parameters are transformed into the appropriate ba-sis to drive a Keplerian RV model (for comparison to theRVs) and the basis used for
TTVFaster computations.Note that this scheme is not possible for an N-body integra-tor, since the initial orbital elements are not the same as the time-averaged orbital elements used in a Keplerian prescrip-tion. We also required Hill stability for all the planets, asdescribed in Equations 1 and 2. In addition, we allowed thestellar mass to vary, using the prior M (cid:63) = 0 . ± . M (cid:12) from S17, in case the combined RV and TTV data addednew information about the stellar mass. The best simulta- The TTVs constrain M k / M (cid:63) , whereas the RVs constrain M k / M / (cid:63) . Figure 4 . The best simultaneous fit to the TTVs and RVs of theWASP-47 system, and residuals. From top to bottom, the panelsshow the TTVs of the transiting planets e ( P = 0 .
79 days, blue),b ( P = 4 .
16 days, green), and d ( P = 9 . neous fit to the TTVs and RVs is shown in Figures 4 (TTVs)and 5 (RVs).Incorporating the priors in Table 3.1, we used emcee toexplore the posteriors and covariances of the dynamical pa-rameters. We ran 100 walkers 10,000 steps each, throwingaway the first 4000 steps as burn-in, and found that our chainshad converged based on the PSRF statistic. (The inclusion ofRV data helped the chains converge faster.) The result of ourMCMC analysis is shown in Figure 6. Our MCMC resultsare summarized in Table 3.To compute planet densities, we utilized the precise stellardensity determined by the photodynamical analysis in Alme- nara et al. (2016). Because the transit of the giant planet hasvery high signal-to-noise, the transit ingress and egress arewell-resolved. The clear ingress and egress and the nearlycircular orbit of the giant planet enable a precise characteri-zation of the stellar limb darkening parameters and the planetimpact parameter. Knowledge of these physical quantities al-lows a precise determination of the stellar density.We translated the precise stellar density into precise planetdensities and radii in the following manner. For each MCMCtrial, we drew a random stellar density from a normal distri-bution N (0 . , . R p = R p R (cid:63) R (cid:63) (4) ρ p = ρ (cid:63) (cid:18) M p M (cid:63) (cid:19) (cid:18) R p R (cid:63) (cid:19) − . (5)The small uncertainty in the stellar density constrains theplanet densities, since the stellar mass and radius (and henceplanet mass and radius) are correlated. Including the stellardensity information reduces the uncertainties in the planetdensities by ∼ INDEPENDENT MULTIPLIED POSTERIORS (IMPS):A QUICK AND DARING WAY TO COMBINEDATASETSIn this section we offer a sanity check of our combinedRV+TTV dynamical solution. Since the RVs and TTVs areindependent observations, the marginalized posteriors fromtheir separate analyses can be multiplied together to estimatethe joint probability distribution of a parameter of interest.This is codified in probability theory as P ( A ∩ B ) = P ( A ) × P ( B ) (6)if A and B are independent events. Assuming that the RVtime series is independent of the TTV times series , the pos-teriors of the RV-only analysis and the TTV-only analysis canbe multiplied together to estimate their joint posterior. Wecall the product of multiplying the posteriors from indepen-dent data sets an Independent, Multiplied Posterior (IMP).The IMP loses some of the information of a simultaneous This is easier to satisfy than the claim that each RV and each TTV is anindependent measurement
Table 3 . Dynamical Parameters from Simultaneous Fit toTTVs (ttvfaster) + RVs (Keplerian)
Parameter Median ± Std. Dev. 95% U.L. Units Ref.
MCMC jump parametersM (cid:63) . ± . M (cid:12) A M e . ± . M ⊕ A M b ± M ⊕ A M d . ± . M ⊕ A M c sin i c ± M ⊕ A P e . ± . P b . ± . P d . ± . P c ± TT e . ± . b A TT b . ± . TT d . ± .
001 KJD A TT c ± √ e e cos ω e . ± .
13 A √ e b cos ω b . ± .
03 A √ e d cos ω d − . ± .
06 A √ e c cos ω c − . ± .
04 A √ e e sin ω e . ± .
13 A √ e b sin ω b . ± .
04 A √ e d sin ω d . ± .
08 A √ e c sin ω c − . ± .
06 A
Parameters from photodynamical analysis ρ (cid:63) . ± .
015 gcm − B R e / R (cid:63) . ± . R b / R (cid:63) . ± . R d / R (cid:63) . ± . Derived ParametersR (cid:63) . ± . R (cid:12) A,B R e . ± . R ⊕ A,B R b . ± . R ⊕ A,B R d . ± . R ⊕ A,B ρ e . ± .
06 gcm − A,B ρ b . ± .
02 gcm − A,B ρ d . ± .
23 gcm − A,B e e . ± . < b A e b . ± . < e d . ± . < e c . ± .
02 A e d cos ω d - e b cos ω b − . ± .
005 A e d sin ω d - e b sin ω b . ± .
007 A ω e . ± . ω b . ± . ω d . ± . ω c . ± . OTE —Results from the MCMC analysis of the TTVs + RVs. The columnsare: parameter, median value plus-or-minus standard deviation, 95% upperlimit (if interesting), and units. A–Derived in this analysis. B–incorporating ρ (cid:63) = 0 . ± .
015 from Almenara et al. (2016). a KJD = BJD - 2454833.0. b –Note that the upper limit on the eccentricity of planet e is determined fromorbital stability requirements, not the measurements. Figure 5 . Left, top: Radial velocities of WASP-47 from four observational campaigns: CORALIE before 2012 (blue squares), CORALIE before2016 (green pentagons), PSF (purple diamonds), and HIRES (red circles). The best-fit model to the TTVs and RVs (fine gray line) is shown.Left, bottom: RV residuals (observations minus the
TTVFaster+radvel model values). The RMS of the residuals is 8.5 m s − , which iscomparable to the mean jitter-enhanced RV uncertainty over all the telescopes (7.7 m s − ). Right: the RVs phase-folded to the orbital periodsof planet e (top), b (second from top), d (second from bottom), and c (bottom). The HIRES RVs are the only single dataset that constrain thesemi-amplitudes of all the planets, because they have the precision (3 m s − ) to capture the small amplitudes of planets e and d, and also thebaseline to capture the amplitude of the long-period planet c. Figure 6 . The mass and eccentricity posteriors of the WASP-47 planets based on a simultaneous fit to the TTVs and RVs (see Equation 3)using the analytic TTV package
TTVFaster and the Keplerian RV package radvel . The masses, orbital periods, times of transit, √ e cos ω , √ e sin ω , and stellar mass were allowed to vary. The MCMC chains shown are available as the Data behind the Figure TTV+RV analysis because the RVs and TTVs are interleavedin time, and because any subtle, slight covariances in the pos-teriors are not captured in the IMP.In Figure 7, we show the posteriors of the planet massesfrom the RV-only (S17) and TTV-only (using the N-body in-tegrator) analyses, and the result of multiplying these pos-teriors together. For planets e and b, the TTVs provide nonew information, and so the IMPs reflect the mass posteri-ors from the RV-only analysis. However, the TTVs alone domeasure the mass of planet d. For planet d, the IMP performsas we would expect: it peaks at a value between the RV-only and TTV-only analyses, and its width is narrower than ei-ther analysis is alone. The IMP for the mass of planet d gives M d = 13 . ± . M ⊕ . This is in good agreement with what wedetermined in the simultaneous modeling of the RVs+TTVs( M d = 13 . ± . M ⊕ ). As expected, the constraint we getfrom the simultaneous modeling of the RVs+TTVs is slightlytighter than the constraint from the IMP. Also, the result ofsimultaneous modeling is slightly closer to the RVs-only so-lution ( M d = 12 . ± . M ⊕ ) than the TTVs-only solution( M d = 16 . ± . M ⊕ ).In cases where one is computation-limited or short on time,1 Figure 7 . Top: posteriors of the mass of WASP-47 e from analysesof the RVs only (blue), the TTVs only (green, using the N-bodyintegrator TTVFast), and their product (i.e. IMP, red). The grayline shows a Gaussian fit to the IMP, with the mean (solid blackline) and 1 σ interval (dashed black lines) shown. Middle: sameas the top, but for WASP-47 b. Bottom: same as the top, but forWASP-47 d. Note that the result of computing the IMP for planetd is in good agreement with the simultaneous RV+TTV analysis(13 . ± . M ⊕ ). the IMP provides an approximate answer. However, if theposteriors are highly covariant in both data sets, the IMPmight grossly overestimate the uncertainties and might alsolose accuracy. This method for combining data sets is con-venient and potentially useful for combining TTV and/or RVdata sets with data from GAIA or WFIRST in the future, butshould be used with great caution. DISCUSSION Here we examine how our joint analysis of the WASP-47 TTVs and RVs provides information about the composi-tions, orbital dynamics, and formation history of the WASP-47 planets.6.1.
Relative Information in Dynamical Analyses
Table 4 summarizes the relative information in various dy-namical analyses of the planet masses and eccentricities. B15modeled the K2 TTVs in a 3-planet N-body analysis in whichthe orbital periods and inital times of transit were fixed, re-sulting in narrow posteriors for the mass of planet b. Themass constraint for planet e comes from the choice of prior,rather than the TTVs. B15 also forward-modeled the systemfor 10 Myr using Mercury (Chambers 1999) to ensure stabil-ity, which resulted in the tight eccentricity constraints for theinner planets: e k < . Table 4 . Relative Information in Literature Dynamical Analyses
Parameter B15 A16.1 A16.2 S16 TTVs-Nbody RVs+TTVs M (cid:63) [ M (cid:12) ] 1 . ± .
08 1 . + . − . . ± .
031 0 . ± .
05 0 . ± .
05 1 . ± . M e [ M ⊕ ] < P . + . − . . + . − . . ± .
17 176 ±
118 9 . ± . M b [ M ⊕ ] 341 + − + − . ± . ±
12 549 ±
252 358 ± M d [ M ⊕ ] 15 . ± . + − . ± . . ± .
70 16 . ± . . ± . M c [ M ⊕ ] – 500 + − + − ±
18 – 416 ± e e < . < .
11 – = 0 P < . P < . P e b < . < .
01 – < . < . < . e d < . < .
024 – = 0 P < . < . e c . ± .
12 – 0 . ± .
04 0 . ± . OTE —B15 - Becker et al. (2015), K2 TTVs, fixed P , TT for each planet, 10 Myr stability enforced. A16.1 -Almenara et al. (2016), photodynamical analysis of K2 TTVs and 71 RVs, A16.2 - including stellar models.S16 - Sinukoff et al. 2016, 118 RVs. TTVs-Nbody - TTVFast N-body analysis (presented herein), K2 TTVs.RVs+TTVs - simultaneous analysis of K2 TTVs and 118 RVs. All upper limits are 95% confidence. P - resultscome from a prior. Simultaneously modeling the TTVs and RVs of planet dyields a more precise determination of the mass of planet dthan can be obtained from either analysis alone: the uncer-tainties shrink from 3 . M ⊕ (TTVs) and 2 . M ⊕ (RVs) to2 . M ⊕ (TTVs + RVs). The TTVs provide no informationabout the mass of planet c, which has a very long period com-pared to the inner planetary system and thus has no effect onthe TTVs. Thus, the RVs provide the majority of the infor-mation about planet masses, although the TTVs contributesubstantially to the mass measurement of planet d.The information about planet eccentricities comes fromstability constraints, the TTVs, and the RVs. The eccentricityof planet e is not constrained by either the TTVs or the RVs,and so its eccentricity varies from 0 to 0.06 (the upper limitfrom stability requirements). The RVs constrain e b < . e b < .
02 (95%confidence), and the combined RVs+TTVs further constrain e b < .
011 (95% confidence). The RVs alone do not providea strong constraint for the eccentricity of planet d (S17 fixed e d = 0). The TTVs alone constrain e d < .
05 (95% confi-dence), and the combined TTVs+RVs constrain e d < . Masses and densities of the planets
WASP-47 is unusual in that the masses of its planets spanalmost two orders of magnitude. The low-mass planets eand d are shown in a density-radius and mass-radius plots for small planets (see Figure 8). The lowest-mass planet,WASP-47 e, is 9 . ± . M ⊕ . At 1 . R ⊕ , this planet has adensity of 9 . ± .
06 g cm − . Planets larger than approx-imately 1 . R ⊕ are unlikely to be rocky (Weiss & Marcy2014; Rogers 2015). Yet, planet e is small and dense enoughthat a rocky composition is likely based on an extrapolationof the empirical relationship for rocky planets smaller than1.5 Earth radii (Weiss & Marcy 2014) and theoretical predic-tions of planet mass and radius for an Earth-like composition(Seager et al. 2007). However, the density of planet e is alsoconsistent with slightly lower densities that might correspondto a rocky interior overlaid with a thin, low-mass envelope ofhigh mean molecular weight materials. Such a compositionhas also been hypothesized for 55 Cancri e, which has a verysimilar mass, radius, and bulk density to WASP-47 e (Lopez2016, S17). Like 55 Cnc e, WASP-47 e is also an ultra-shortperiod planet cohabiting an orbital system with giant plan-ets, which underscores the question of how planetary systemarchitecture and planet compositions are related.By contrast, WASP-47 d, which is 3 . R ⊕ , has a mass of13 . ± . M ⊕ . This is a slightly higher mass than was re-ported in S17 because the TTVs add mass information. Theadditional information from the TTVs also narrows the massposterior, shrinking the uncertainty from 2 . . M ⊕ . Thedensity of WASP-47 d is 1 . ± .
23 g cm − , making it ahigh-density member of the population of sub-Neptune sizedplanets with volatile envelopes. The mass-radius diagram inFigure 8 shows that the mass, radius, and density of WASP-47 d make it one of the most Neptune-like planets discov-ered to date. While its composition could be explained bya two-layer model of a H/He envelope atop a silicate-ironcore, a Neptune-like composition that includes a thick layerof super-ionic water might also explain the bulk properties ofWASP-47 d.WASP-47 b is a Jupiter-mass planet (358 ± M ⊕ ) that re-3ceives 440 ±
70 times as much incident stellar irradiation asthe Earth does. At 13 . ± . R ⊕ , the planet has a typicaldensity (1 . ± .
02 g cm − ) for its mass and incident stellarflux (see Figure 9), consistent with various theories (Baty-gin et al. 2011; Fortney & Nettelmann 2010, and referencestherein) that stellar irradiation inflates the planet and/or pre-vents the planet from cooling.6.3. Eccentricities of the planets
The orbits of the three inner planets are profoundly cir-cular. B15 found that eccentricities of < .
06 were re-quired for stability. Here, we tighten the eccentricities to e b < .
011 and e d < .
025 (95% confidence). In the highest-eccentricity cases for planets b and d, they tend to be ap-sidally aligned. We compute e d cos ω d − e b cos ω b = − . ± .
005 and e d sin ω d − e b sin ω b = 0 . ± . ∼ − years, depending on the tidal Q valuefor the planet, our N-body analysis revealed that the neigh-boring giant planet (b) perturbs the eccentricity of planet eon a timescale of 1.26 days. This timescale happens to berelated to the orbital periods of both e and b by 1 / P kick =(1 / P e + / P b ) − . Since there is no commensurability betweenthe orbital periods of b and e, we expect that over longtimescales, the kicks from planet b average out, and so the ar-gument of periastron of planet e should not have a preferredvalue. However, planet e could still have an average eccen-tricity that is higher than zero. Likewise, planet d shouldnot be assumed to have zero eccentricity because it exhibitsTTVs. Because planet d is near the 2:1 resonance with planetb, its argument of periastron circulates at the frequency of theTTV super-period.The very low eccentricities of the WASP-47 planets areremarkably like those of the solar system planets (see Fig-ure 10). The average eccentricity of the detected planets inWASP-47 is < .
09; in the solar system, the average eccen-tricity of the planets and Pluto is 0.08.6.4.
Observational Constraints on Formation Theories
Compact systems like the WASP-47 inner planets mightform in situ from protoplanetary disks that are more mas-sive than the minimum mass solar nebula (Chiang & Laugh-lin 2013). Scenarios in which either gas-poor sub-Neptunesor gas-rich Jupiters form in situ have been proposed (Leeet al. 2014; Lee & Chiang 2015; Lee & Chiang 2016; Batyginet al. 2016). If all the WASP-47 planets formed in situ fromthe same nebular material, the challenge is to explain howWASP-47 b managed to achieve runaway growth, whereasits immediate neighbors remained gas-poor. Whether aJupiter-mass or sub-Saturn mass planet forms depends onthe core mass, atmospheric opacity, and the disk lifetime(Hori & Ikoma 2011; Venturini et al. 2015). Cores of mass ∼ − M ⊕ and larger undergo runaway gas accretion (withthe exact value of critical core mass depending on the at- mospheric opacities, metallicities and the disk gas dissipa-tion timescale; Lee & Chiang 2016; Batygin et al. 2016).However, gas-poor sub-Neptunes are formed instead of giantplanets if the cores form and accrete their envelopes in a gas-poor (transitional) disk. The formation of sub-Neptune coresby giant impact requires at least four orders of magnitude ofgas-depletion with respect to the minimum mass extrasolarnebula (Lee & Chiang 2016), leaving ∼ . M ⊕ of gas in thedisk. This gas budget is insufficient to form the ∼ M ⊕ ofgas in WASP-47 b.A modification of in situ formation is inside-out growthvia pebble accretion, in which pebbles are transported in-ward through the disk until they are stopped by a pressuremaximum at the boundary between the magneto-rotationalinstability (MRI) zone and the magnetic dead zone (Chatter-jee & Tan 2014). At this boundary, the pebbles may becomeToomre unstable or may coalesce via core accretion. As thegas disk clears and the boundary between the MRI and deadzone moves outward, the site of planet formation graduallymoves out through the disk. Because the hot Jupiter is situ-ated between two sub-Neptunes, the accretion rate of the diskwould likely need to increase by roughly an order of magni-tude, and then decrease again, to explain the high mass of theJupiter compared to its neighbors. Furthermore, the efficacyof forming hot Jupiters via inside-out planet formation hasnot been well studied.Alternatively, the planets could have formed elsewhere inthe disk and then migrated via interactions with the disk totheir present locations. In both Type I and Type II migra-tion, the gas disk damps planet eccentricity, allowing planetsto maintain circular orbits in a manner consistent with thenearly circular orbits of the three inner planets. However,slow migration can trap the planets in mean motion reso-nances. While WASP-47 b and d are near the 2:1 mean mo-tion resonance, they are not trapped. Thus, their migrationhistory must either include a mechanism to prevent planets band d from entering the 2:1 resonance, or remove them fromthe resonance (e.g., Adams et al. 2008; Goldreich & Schlicht-ing 2014). In particular, Deck & Batygin (2015) find that ifthe inner planet of a pair near a first-order mean motion reso-nance is the more massive (as is the case for WASP-47 b andd), escape from resonance is unlikely. Furthermore, migra-tion in the disk does not explain the eccentricity of planet c(0 . ± . ◦ . The planets swapangular momentum, causing dramatic variations in the in-clinations and eccentricities of the planets over time (Kozai4 Figure 8 . Left: planet density versus planet physical radius for 94 transiting planets smaller than 4 . R ⊕ . The gray circles have massesdetermined from RVs; the gold circles have masses determined from TTVs. The size of the circle corresponds to 1/ σ ρ . Blue squares show theweighted mean density in bins of 0 . R ⊕ to guide the eye. The blue diamonds are the solar system planets. The red dashed line is an empiricallinear fit to planet density versus radius for the exoplanets and solar system planets smaller than 1 . R ⊕ , extended to predict the densities ofpotentially rocky planets larger than 1 . R ⊕ . For comparison, we show the predicted density-radius curve for a polytropic equation of state ofan Earth-composition planet (Seager et al. 2007, green dotted line). The black line is an empirical power-law fit to planet mass versus radius forplanets larger than 1 . R ⊕ . WASP-47 e (1 . R ⊕ ) sits on the red line and therefore is consistent with a rocky composition, but could also havea thin volatile envelope. WASP-47 d (3 . R ⊕ ) has a density that requires significant volatiles, including the possibility of high-density volatilesbased on the similarity of its bulk properties to Uranus and Neptune. Right: Same as the left panel, but showing planet mass versus radius, andthe circle sizes correspond to 1 / σ m . v sin i and the planet-star obliquity allowsmoderate misalignments ( | λ | < ◦ with 1 σ confidence). Ifplanet c is coplanar and the spin-orbit alignments of the plan-ets are small, a coplanar, high-eccentricity migration (HEM)scenario such as the one proposed in Petrovich (2015) mightbest explain the present locations of the planets. However,the current nearly-circular, coplanar orbits of the inner plan-etary system are not consistent with a past modest eccentric-ity and/or inclination for planet b. Recall from the stabilityanalysis in B15 that the nearly all of the scenarios in whichthe eccentricity of planet b exceeds 0.06 become unstablewithin 10 Myr. Therefore, although high eccentricity migra-tion mechanisms can explain the present positions of planetsb and c, such mechanisms would have destroyed the small,close-in planets. Indeed, simulations of giant planets migrat-ing inward find that the giant planets collide with low-massinner planets in both high-eccentricity (Mustill et al. 2015) and low-eccentricity (Batygin & Laughlin 2015) scenarios.6.5. The Two Stage Planet Formation Hypothesis
One way to explain the present orbits and compositionsof the WASP-47 planets is that the planets did not all format the same time. As a point of reference, consider the so-lar system. In our own solar system, the giant planets musthave formed early (within 1 to 10 Myr), when the proto-planetary disk was still gas-rich (de Pater & Lissauer 2001).In contrast, the formation of terrestrial-mass planets by plan-etesimal accretion can take as long as ∼ years in N-bodysimulations and is more efficient in gas-poor disks where theplanetesimal eccentricities can grow, leading to more colli-sions (Lissauer 1987; Pollack et al. 1996; Ida & Lin 2004;Lee & Chiang 2016). Thus, it is possible that the solar sys-tem giant planets formed and migrated to their present loca-tions before the terrestrial planets formed. The Nice Modeldemonstrates that early formation and migration of Jupiterand Saturn can reproduce the current orbital architectures ofthe gas giants and ice giants while also accounting for theperiod of Late Heavy Bombardment on the terrestrial plan-ets, and Trojan asteroids (Gomes et al. 2005; Tsiganis et al.2005; Morbidelli et al. 2005). The Grand Tack Model shows5 Figure 9 . Left: planet radius vs. mass for exoplanets with measured masses and radii, as determined by querying exoplanets.org on 2017/02/19and including radii and masses from (Hadden & Lithwick 2016) and (Gillon et al. 2017). The sample is divided into those that receive morethan the median incident flux ( F > F ⊕ , red points) and those that receive less than the median incident flux ( F < F ⊕ , blue points). Thesolar system planets are labeled. The WASP-47 planets are shown (yellow stars, the mass and radius error bars are smaller than the symbols).The sub-Neptunes WASP-47 e and d are high-density for their size. WASP-47 b is a typical-sized hot Jupiter for its mass and incident flux. Figure 10 . The mean eccentricities of the WASP-47 and solar systemplanets versus the orbital distances, scaled such that a b = a J to facil-itate comparison. The point size corresponds to M / , to illustratethe mass range while keeping all the planets visible. The spread ofthe each WASP-47 planet eccentricity is illustrated with a cloud of1000 small red points drawn from the eccentricity posterior. The av-erage eccentricity of the W47 planets is comparable to the averageeccentricity of the solar system planets, if we include Pluto. that a reversal in the direction of Jupiter’s migration due totorques from Saturn can explain the small size of Mars anddetailed compositional features of the asteroids (Walsh et al.2011; Morbidelli et al. 2012). Both of these models featuretwo stages of planet formation:1. Early giant planet growth combined with disk and/orplanet-induced migration produces the current compo- sitions of giant planets and sculpts the radial distribu-tion of planetesimals.2. When or after the gas disk clears, planetesimal accre-tion proceeds in the sculpted planetesimal disk, result-ing in the formation of low-mass, predominantly rockyplanets.Likewise, the formation and evolution of the planets inWASP-47 might be best explained by a two-stage process.The giant planets in the WASP-47 system might have formedearly in the gas-rich disk and exchanged energy with eachother, additional bodies, or the disk to migrate to their presentlocations. If the gas disk cleared while WASP-47 b and cwere finishing their migration, the final years of the giantplanets’ migration could have influenced the solid material inthe disk, exciting protoplanetary solids to higher eccentrici-ties and inducing a second stage of core accretion. Such anepoch could have formed WASP-47 e and d, and perhaps ad-ditional yet-undetected low-mass planets dominated by rockyinteriors.6.6. Predictions of Two-Stage Planet Formation
In WASP-47, the orbital architecture and compositions ofthe low-mass planets provide insight into the formation his-tory of the giant planets. Because the low-mass planets needa small amount of gas to form (about 0 . M ⊕ ), their forma-tion must occur before the gas disk completely dissipates.This puts a deadline on when planet b can finish its migra-tion; it must park at its current location before planets e and6d form (otherwise it destroys them). In other words, planetb must arrive at its present location before the gas disk dissi-pates.Because a gas disk damps planet eccentricities on thetimescale of 10 years (Papaloizou & Larwood 2000), high-eccentricity migration cannot proceed in a gas disk. Further-more, the timescale for type I migration in a gas disk is muchfaster than the timescale for HEM, and so gas disk migra-tion dominates while the gas disk is present. If we considera mostly gas-depleted disk, we can try to migrate WASP-47b inward via HEM in a short window of time before the diskcompletely dissipates. However, some force must circular-ize planet b’s orbit quickly enough to allow the neighboringsmall planets to form. In the absence of gas, the only circu-larizing force is tides raised on the planet by the star. Thetimescale for eccentricity damping due to tides (Murray &Dermott 2000, Equation 4.198) is: τ e = − e ˙ e = 463 M p M (cid:63) (cid:0) aR p (cid:1) µ pl Q pl π/ P . (7)where µ pl is the ratio of the elastic to gravitational forces inthe planet ≈ (10 km / R p ) , and Q pl is the tidal dissipationfactor of the planet. For a = 0 .
05 AU and Q pl = 10 . (typicalfor hot Jupiters, Jackson et al. 2008), the circularization timeis τ e ≈ yr, which is an order of magnitude longer than thedisk lifetime. There is not enough time to move WASP-47 bto its present position via HEM and then circularize its orbitbefore the little planets form. Therefore, a history of HEM ishighly unlikely for planet b.On the other hand, planet b could also have formed in situ or nearly in situ in the gas-rich disk. Antonini et al. (2016)find through N-body simulations that, for the majority ofobserved Jupiter pairs (i.e., systems with one Jupiter-massplanet inside 1 A.U. and one Jupiter-mass planet outside 1A.U.), the inner planet is unlikely to have formed beyond 1A.U. Thus, it seems probable that WASP-47 b formed inside1 A.U., and thus inside the snow line.Since planet c has modest eccentricity (0 . ± . e < . Kepler will enable studies of the occurrence of planetary systemswith diverse compositions like our solar system, illuminatingthe dominant physical processes in their formation.Continued observations of this system and other systemswith diverse planet compositions from the James WebbSpace Telescope, the WFIRST mission, and other facilitieswill provide evidence for or against a two-stage planet for-mation scenario. CONCLUSIONWe combine 118 RVs and 108 K2 TTVs to present someof the most precise masses, densities, and orbital dynam-ics of the WASP-47 planetary system to date. For thetransiting inner planetary system, we obtain M e = 9 . ± . M ⊕ ( ρ e = 9 . ± .
06 g cm − ), M b = 358 ± M ⊕ ( ρ b =1 . ± .
02 g cm − ), and M d = 13 . ± . M ⊕ ( ρ d = 1 . ± .
23 g cm − ).We place tight upper limits on the eccentricities of the threeinner planets: e e < . e b < . e d < .
025 (95%confidence). The mean eccentricity of the WASP-47 planets(0.09) is very similar to the mean eccentricity of the solarsystem planets plus Pluto (0.08).Of the four planets, only WASP-47 e has bulk propertiesconsistent with a rocky composition. WASP-47 d is a sub-Neptune sized planet that is just a little smaller than Neptunebut has the density of Neptune, meaning that its bulk prop-erties consistent with a Neptune-like composition. However,WASP-47 d could also be a rocky interior overlaid with ahydrogen-rich, water-poor envelope. WASP-47 b is a hotJupiter with a typical size for its mass and incident stellarflux.The non-resonant architecture, diverse masses, and pro-foundly circular orbits in the WASP-47 resemble our own so-lar system. We briefly review the highlights of in situ planetformation, inside-out planet formation, disk migration, andhigh-eccentricity migration, noting that none of these mech-anisms alone can reproduce all the physical attributes of theWASP-47 system. We propose that WASP-47, like the solarsystem, formed in two stages. In stage one, the giant planetsformed in a gas-rich disk and migrated via disk migration totheir present locations. In stage two, the high-density sub-Neptunes formed in situ in a gas-poor environment.7L. M. W. acknowledges the Trottier Family Foundationfor their support. K. D. acknowledges the support of theJCPA fellowship at Caltech. E. S. is supported by a post-graduate scholarship from the Natural Sciences and Engi-neering Research Council of Canada. E. A. acknowledgessupport from NASA grants NNX13AF20G, NNX13AF62G,and NASA Astrobiology Institutes Virtual Planetary Lab-oratory, supported by NASA under cooperative agreementNNH05ZDA001C. A. W. H. acknowledges support from aNASA Astrophysics Data Analysis Program grant, supportfrom the K2 Guest Observer Program, and a NASA KeyStrategic Mission Support Project. We thank Simon Walkerfor his help and generosity in using ttvfast-python. The au-thors thank Geoff Marcy, Eva Culakova, Daniel Fabrycky,Jack Lissauer, Jason Rowe, the Kepler-TTV working group, and the KITP Planet Formation and Dynamics workshop foruseful discussions. We thank NASA and the Kepler teamfor the outstanding photometry that contributed to this paper.The authors wish to extend special thanks to those of Hawai-ian ancestry on whose sacred mountain of Maunakea we areprivileged to be guests. Without their generous hospitality,the Keck observations analyzed herein would not have beenpossible.
Facilities:
Keck:I (HIRES), Kepler
Software:
TTVFaster (Agol & Deck 2016a,b),TTVFast (Deck et al. 2014), ttvfast-python(https://github.com/mindriot101/ttvfast-python), emcee(Foreman-Mackey et al. 2013), lmfit (Newville et al. 2014),radvel (https://github.com/California-Planet-Search/radvel)REFERENCES
Adams, F. C., Laughlin, G., & Bloch, A. M. 2008, ApJ, 683, 1117Agol, E., & Deck, K. 2016a, ApJ, 818, 177—. 2016b, TTVFaster: First order eccentricity transit timing variations(TTVs), Astrophysics Source Code Library, ascl:1604.012Almenara, J. M., Díaz, R. F., Bonfils, X., & Udry, S. 2016, A&A, 595, L5Antonini, F., Hamers, A. S., & Lithwick, Y. 2016, AJ, 152, 174Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, The AstrophysicalJournal Supplement Series, 204, 24Batygin, K., Bodenheimer, P. H., & Laughlin, G. P. 2016, ApJ, 829, 114Batygin, K., & Laughlin, G. 2015, Proceedings of the National Academy ofSciences, 112, 4214Batygin, K., Stevenson, D. J., & Bodenheimer, P. H. 2011, ApJ, 738, 1Becker, J. C., Vanderburg, A., Adams, F. C., Rappaport, S. A., &Schwengeler, H. M. 2015, ApJL, 812, L18Borucki, W. J., Koch, D. G., Brown, T. M., et al. 2010, ApJL, 713, L126Brooks, S. P., & Gelman, A. 1998, Journal of computational and graphicalstatistics, 7, 434Bryan, M. L., Knutson, H. A., Howard, A. W., et al. 2016, ApJ, 821, 89Burke, C. J., Christiansen, J. L., Mullally, F., et al. 2015, ApJ, 809, 8Chambers, J. E. 1999, MNRAS, 304, 793Chatterjee, S., & Tan, J. C. 2014, ApJ, 780, 53Chiang, E., & Laughlin, G. 2013, Monthly Notices of the RoyalAstronomical Society, 431, 3444Dai, F., Winn, J. N., Arriagada, P., et al. 2015, ApJL, 813, L9de Pater, I., & Lissauer, J. J. 2001, Planetary Sciences, 544Deck, K. M., Agol, E., Holman, M. J., & Nesvorný, D. 2014, TheAstrophysical Journal, 787, 132Deck, K. M., & Batygin, K. 2015, ApJ, 810, 119Dotter, A., Chaboyer, B., Jevremovi´c, D., et al. 2008, ApJS, 178, 89Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP,125, 306Fortney, J. J., & Nettelmann, N. 2010, SSRv, 152, 423Fressin, F., Torres, G., Charbonneau, D., et al. 2013, The AstrophysicalJournal, 766, 81Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542,456Gladman, B. 1993, Icarus, 106, 247Goldreich, P., & Schlichting, H. E. 2014, AJ, 147, 32Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature,435, 466 Hadden, S., & Lithwick, Y. 2016, ArXiv e-prints, arXiv:1611.03516Hellier, C., Anderson, D. R., Collier Cameron, A., et al. 2012, MNRAS,426, 739Hori, Y., & Ikoma, M. 2011, MNRAS, 416, 1419Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388Jackson, B., Greenberg, R., & Barnes, R. 2008, ApJ, 678, 1396Kipping, D. M. 2014, MNRAS, 440, 2164Koch, D. G., Borucki, W. J., Rowe, J. F., et al. 2010, ApJL, 713, L131Kozai, Y. 1962, AJ, 67, 591Lee, E. J., & Chiang, E. 2015, ApJ, 811, 41Lee, E. J., & Chiang, E. 2016, The Astrophysical Journal, 817, 90Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95Lidov, M. L. 1962, Planet. Space Sci., 9, 719Lissauer, J. J. 1987, Icarus, 69, 249Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, ApJ, 750, 112Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122Lopez, E. D. 2016, ArXiv e-prints, arXiv:1610.01170Lovis, C., & Fischer, D. 2010, Radial Velocity Techniques for Exoplanets,ed. S. Seager, 27–53Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005, Nature,435, 462Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N., & Walsh,K. J. 2012, Annual Review of Earth and Planetary Sciences, 40, 251Murray, C. D., & Correia, A. C. M. 2010, Keplerian Orbits and Dynamicsof Exoplanets, ed. S. Seager, 15–23Murray, C. D., & Dermott, S. F. 2000, in Solar System Dynamics(Cambridge University Press), 321–408, cambridge Books OnlineMustill, A. J., Davies, M. B., & Johansen, A. 2015, ApJ, 808, 14Neveu-VanMalle, M., Queloz, D., Anderson, D. R., et al. 2016, A&A, 586,A93Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT:Non-Linear Least-Square Minimization and Curve-Fitting for Python½u,doi:10.5281/zenodo.11813Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013a, Proceedings of theNational Academy of Sciences, 110, 19273Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013b, ApJ, 770, 69Petrovich, C. 2015, ApJ, 805, 75Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62Rogers, L. A. 2015, ApJ, 801, 41Rowe, J. F., Bryson, S. T., Marcy, G. W., et al. 2014, ApJ, 784, 45Sanchis-Ojeda, R., Rappaport, S., Winn, J. N., et al. 2013, TheAstrophysical Journal, 774, 548