A Simplified Photodynamical Model for Planetary Mass Determination in Low-Eccentricity Multi-Transiting Systems
DDraft version November 10, 2020
Typeset using L A TEX default style in AASTeX63
A Simplified Photodynamical Model for Planetary Mass Determination in Low-EccentricityMulti-Transiting Systems
Gideon Yoffe,
1, 2
Aviv Ofir, and Oded Aharonson
1, 3 Department of Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot, 76100, Israel Now at: Max Planck Institute for Astronomy, K¨onigstuhl 17, 69117 Heidelberg, Germany Planetary Science Institute, Tucson, AZ, 85719, USA (Accepted November 6, 2020)
Submitted to ApJABSTRACTInferring planetary parameters from transit timing variations is challenging for small exoplanetsbecause their transits may be so weak that determination of individual transit timing is difficult orimpossible. We implement a useful combination of tools which together provide a numerically fastglobal photodynamical model. This is used to fit the TTV-bearing light-curve, in order to constrainthe masses of transiting exoplanets in low eccentricity, multi-planet systems - and small planets inparticular. We present inferred dynamical masses and orbital eccentricities in four multi-planet systemsfrom
Kepler ’s complete long-cadence data set. We test our model against Kepler-36 / KOI-277, asystem with some of the most precisely determined planetary masses through TTV inversion methods,and find masses of 5.56 +0 . − . and 9.76 +0 . − . m ⊕ for Kepler-36 b and c, respectively – consistent withliterature in both value and error. We then improve the mass determination of the four planets inKepler-79 / KOI-152, where literature values were physically problematic to 12.5 +4 . − . , 9.5 +2 . − . , 11.3 +2 . − . and 6.3 +1 . − . m ⊕ for Kepler-79 b, c, d and e, respectively. We provide new mass constraints where noneexisted before for two systems. These are 12.5 +3 . − . m ⊕ for Kepler-450 c, and 3.3 +1 . − . and 17.4 +7 . − . m ⊕ for Kepler-595 c (previously KOI-547.03) and b, respectively. The photodynamical code used here,called PyDynamicaLC , is made publicly available.
Keywords:
Exoplanets,
Kepler , Photodynamics INTRODUCTIONTransit timing variations (TTVs) caused by the gravitational interactions in multi-planet systems (Holman 2005;Agol et al. 2005) may be modeled to determine the governing system parameters, in particular masses. This isespecially valuable in cases where no additional observations beyond photometry, such as radial velocity data, areavailable. In practice, performing this inversion is challenging in the case of small planets and low-amplitude TTVs,although these commonly occur in nature and
Kepler data. To address this situation we combine several techniquesfrom the literature.
Global model : The standard technique for the detection of TTVs follows two steps: firstly, one optimizes anempirical light-curve model using common shape parameters for all of the transit events, (in the usual Mandel & Agol(2002) framework these are a/r (cid:63) , b/r (cid:63) , and r p /r (cid:63) , the planetary semi-major axis, impact parameter and radius, eachnormalized by the stellar radius) but assigns an individual time for each transit event. Secondly, the list of best-fittingtransit times is modeled with a dynamical description. There are shortcomings to this approach: the transits mustbe deep enough to be individually significant, they must be long enough to have more sufficient points that determinetheir timing, and the number of fitting parameters grows with the number of transit, greatly diminishing the sensitivity Corresponding author: Gideon Yoffeyoff[email protected] a r X i v : . [ a s t r o - ph . E P ] N ov Yoffe et al. of this approach. To rectify these shortcomings Ofir et al. (2018) recently introduced a more sensitive technique forthe detection of TTVs: one assumes sinusoidal shape for the TTVs, and a global optimization allows determination ofa single set of parameters describing the full light-curve variations. This approach is effective for small transit depthsor durations (even in cases where individual transits cannot be seen), and has significantly fewer degrees of freedom -enhancing it sensitivity, discussed next.
Parsimonious model : The problem of inverting the observed TTVs for the component masses is non-trivial. Thedegeneracy between the effects of masses and orbital eccentricities has been appreciated in early works on the subject(Holman 2005; Agol et al. 2005) and a solution based on a full n-body integration is highly multi-dimensional, evenmore degenerate than just explained. It is thus difficult to initialize and optimize, and is computationally expensive oreven prohibitive. More recently, it was shown that the pattern of the TTVs does not depend strongly on the individualeccentricities but rather only on the eccentricity-vector difference between adjacent planetary pairs (Hadden & Lithwick2016; Jontof-Hutter et al. 2016; Linial et al. 2018). Fitting for the eccentricity differences allows capturing most of theTTV variability using two fewer degrees of freedom.We use
TTVFaster (Agol & Deck 2016) as our dynamical tool of choice.
TTVFaster is a semi-analytic modelaccurate to first order in eccentricity which approximates TTVs using a series expansion. An important simplificationin using this tool arises from its reliance only on average
Keplerian elements, which are well determined from previousprocessing steps, thus significantly decreasing the number of degrees of freedom in the model to just three per planet– its mass and its two perpendicular eccentricity components ( e x , e y , or e cos ω, e sin ω ), where ω is the argument ofperiapse and e is the eccentricity magnitude. Consequently, there are less correlations among parameters, producingbetter determined values, and initial guesses are more easily selected. Additionally, TTVFaster offers a considerablyfaster computation time than n-body integrators (e.g. about an order of magnitude faster than
TTVFast (Deck et al.2014)), resulting in sampler convergence timescales on the order of hours, rather than weeks (Tuchow et al. 2019).This will allow future expansion of the scope of our study to the entire sample of TTV-bearing multi-transiting
Kepler systems.An important disadvantage of using
TTVFaster is that its underlying approximations are inaccurate for systemswith high eccentricities. Therefore, we require that our model pass a series of tests to verify the validity of our modelin the relevant region of the parameter space. We also note that performing fits in flux space (the light-curve) doesnot permit using the linear decomposition that is available in the timing space (Linial et al. 2018).In this work we therefore build on both lessons above, and fit a model which is both global and minimal in degreesof freedom to achieve better precision and sensitivity than previous studies, and in a uniform fashion to multiplerelevant
Kepler targets. We also use the recently updated stellar parameters based on Gaia DR2 (Fulton & Petigura2018), which allows to better derive the planetary physical properties. The paper is divided as follows: § § § § § DATA PROCESSINGThe sample of
Kepler systems we processed is all the multi-transiting systems of which at least one component wasidentified as exhibiting TTVs by Ofir et al. (2018) - a total of 159 systems and 466 planets. In the first processing stage
Kepler ’s photometery was iteratively detrended and the light-curve fitted using empirical TTVs (as detailed below)- to ensure the planet signals themselves do not interfere with detrending. Importantly, the goal of this stage is not to produce a list of transit times, but only to produce the best possible detrended and normalized light-curve. Thegravitational interactions among the objects are modeled in the subsequent step.For the light-curve detrending, each of our target systems was modeled empirically in an evolved method based onOfir & Dreizler (2013). All the planets in the system where fitted simultaneously using a common value a/r (cid:63) , where a is the semi-major axis and r (cid:63) is the stellar radius. This value was scaled between planets using Kepler’s ThirdLaw. TTVs, if exist, were accounted for by a phase shift of each event. This is a circular-orbit approximation of theTTVs. The detrending curve itself was the better of either a Savitzky-Golay filter or a cosine filter - applied to eachsection of continuous data independently. Each of the filters was applied several times - as many as there are knowntransiting planets in the system - since each transit has a unique duration. In every continuous section the transitwith the longest duration appearing on that section determined the critical time-scale for filtering, and each of thefilters (Savitzky-Golay or cosine) used a time-scale no shorter than three times the critical time-scale on that section. Simplified Photodynamical Model for Planetary Mass Determination P , a mean reference epoch t lin (i.e. the constant term in the linear fit tothe times of transit), b/r (cid:63) , r p /r (cid:63) and a/r (cid:63) for the innermost planet. Note, the empirically-fitted individual times ofmid-transit were not used in the dynamical modeling below. PHOTODYNAMICAL MODELFor each multi-planet system, we compute a photodynamical model in three steps: first, we simulate a TTV signalwith
TTVFaster . The required input for each planet i are the average Keplerian elements, in addition to the planetarymasses and the first time of mid-transit for a linearly approximated orbit ( P i , i i , e x,i , e y,i , t lin ,i , m i ), where i refers to the index of the planet, P i is the orbital period, i is the orbital inclination that can be converted from b witheq. 7 in Winn et al. (2011), e x,i and e y,i are the perpendicular eccentricity components, respectively, t lin ,i is the meanreference epoch and m i is the mass. Note, throughout this study the Keplerian parameters refer to their time averagesover Kepler long-cadence observations (as opposed to their osculating values).Excluding masses and eccentricities, these parameters are well constrained during the data processing stage ( § n pl degrees of freedom for the dynamical fit. The output obtained from TTVFaster is a vector of timesof mid-transit, from which TTVs may be deduced.Each individual transit event is then shifted from a strict periodic regime and generate the individual-event transitlight-curve Mandel-Agol model (Mandel & Agol 2002), using
PyAstronomy ’s implementation of
EXOFAST (Eastmanet al. 2013). Assuming small orbital eccentricity allows us to further simplify our model and approximate the shape ofthe transit curve to be that of an arc of a circular orbit (referred to as a quasi-circular approximation). An additionalassumption throughout this study is that all the TTVs can be attributed to mutual perturbations among adjacentpairs of known transiting planets. While TTVs can be sensitive to mutual inclinations (Payne et al. 2010), Kepler multi-planet systems are usually characterized by small values of these (Xie et al. 2016), such that this effect does notmake a notable contribution to TTVs of planets near low-order MMR (Payne et al. 2010).Finally, we perform non-linear optimization of the dynamical parameters and estimate their uncertainties using
MultiNest (Feroz et al. 2008), a multi-modal Markov Chain Monte Carlo algorithm with implementation in Python(Buchner et al. 2014).
MultiNest implements nested sampling (Skilling 2004), a Monte Carlo method targeted at theefficient calculation of the Bayesian evidence. Using
MultiNest allows for the exploration of regions of interest in theparameter space with a reduced risk of being trapped in local likelihood maxima, which is especially important dueto the known mass-eccentricity degeneracy (Lithwick et al. 2012). MODEL OPTIMIZATIONFor any multi-planet model, we fit 3 n pl free parameters and also use 4 n pl + 5 input constants. The latter are thelinear ephemeris (see §
2) and stellar parameters - m (cid:63) and r (cid:63) , two stellar limb darkening coefficients and a of theinnermost planet. Stellar parameters are adopted from Fulton & Petigura (2018). Yoffe et al.
The free parameters for each planet are its mass and its two perpendicular eccentricity vector components, e y and e x (along the line of sight, and perpendicular to it along the planetary motion direction, respectively). As mentioned,in the case of small eccentricities, TTVs depend chiefly on the difference in the eccentricity vectors of each adjacentplanetary pair, rather than on the individual eccentricities. Attempting to fit for both individual eccentricities yieldshighly degenerate results (e.g. Jontof-Hutter et al. (2016, their figures 3,5,7,9,11,15,17,19)). Fitting for the eccentricitydifference, however, results in considerably weaker mass-eccentricity degeneracy (see correlation plots in the Appendix).Therefore, we fit the absolute orbital eccentricity ( e x, , e y, ) for only the innermost planet (expecting it to remainpoorly constrained), and for the remaining planets we fit for the eccentricity differences (∆ e x,i +1 = e x,i +1 − e x,i ,∆ e y,i +1 = e y,i +1 − e y,i ). 4.1. Statistical Model
We evaluated the goodness of fit to the multi-planet light-curve for each simulated model with parameter set Θ usingthe log of Student’s t -distribution function with 4 degrees of freedom, whose log islog( L ) = − n (cid:88) i =1 log (cid:16) √ σ i (cid:17) − n (cid:88) i =1 log (cid:18) y i, Θ − x i ) σ i (cid:19) , (1)where there are n observed photometric measurements x i with respective uncertainties σ i , and y i is the simulatedmulti-planet light-curve associated with the model parameters Θ.We justify the choice of our likelihood function – Student’s t-distribution function with 4 degrees of freedom –by comparing the distribution of residuals for all measured and best-fit model transit times in this study with twonormalized and symmetric statistical distributions, a Gaussian, and Student’s t -distribution with an optimized numberof degrees of freedom. We then use the Kolmogorov-Smirnov test to determine the preferred likelihood distribution,parameterized by the smallest value of D - the Kolmogorov-Smirnov statistic for a given cumulative distribution. Inour analysis we find that similarly to Jontof-Hutter et al. (2016); MacDonald et al. (2016), the wings of the distributioncontain more outliers than would be expected by Gaussian uncertainties, such that the Student’s t-distribution scoredbetter in the Kolmogorov-Smirnov test. 4.2. Priors
For the planetary masses, we use linear priors spanning the range 10 − m ⊕ ≤ m p ≤ m ⊕ .For the eccentricities, we test three Gaussian priors in the individual eccentricity components of the inner planet anddifferences in eccentricity vector components ( e x, , e y, , ∆ e x,i , ∆ e y,i ) , with zero mean and three different variances σ e = 0.01, 0.02 and 0.05. These correspond to Rayleigh distributions of the eccentricity magnitude with scale-widths of0.01, 0.02 and 0.05, respectively. Note that while the optimized parameters are the eccentricity differences, and indeedthe TTV signal depends mostly on these alone, TTVFaster requires the absolute eccentricities of all components.The choice of eccentricity priors is consistent with the distribution of
Kepler multi-planet systems eccentricitiesreported by Moorhead et al. (2011) and Xie et al. (2016), both conducting independent analyses through measuredtransit durations. 4.3.
Posteriors
Our quoted values are the median values of the posterior distribution, and the uncertainties are the 50 ± . < − times that of the best-fit). We refer to the range of uncertainties as σ , defined irrespectiveof the distribution. For a normal distribution, σ corresponds to the standard deviation.Performing three optimizations with three eccentricity priors enables us to test the sensitivity of our optimizationresults and their respective uncertainties to the choice of eccentricity priors (following e.g. Hadden & Lithwick (2014,2016, 2017)). A result for a system in which the best-fit eccentricity values are all small across all prior widths indicatesthat the true eccentricity is indeed small and our result is not limited by the prior.4.4. Validity map https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.kstest.html Simplified Photodynamical Model for Planetary Mass Determination
TTVFaster to thatof
TTVFast – a symplectic n-body integrator with a Keplerian interpolator for the calculation of transit timings, itselfverified to be in good agreement with
Mercury (Chambers 1999). The analysis is done using a validity map we devised,in the form of a two-dimensional grid in ∆ e x and ∆ e y for every planet pair, keeping the remaining parameters constantat their best-fit values (producing n pl − TTVFast and
TTVFaster normalized by the data’s error estimation are used to construct a χ -like statistic: χ = n (cid:88) i =1 (cid:18) y i, TTVFast − y i, TTVFaster σ i (cid:19) (2)This allows evaluating the region in parameter space in which TTVFast and
TTVFaster are consistent to within themeasurement error. A difficulty in performing the analysis lies in the fact that only the average Keplerian parametersare known (with uncertainty), whereas
TTVFast requires the instantaneous values at the beginning of integration. Wetherefore use the average Keplerian elements as initial conditions, from which revised average values are computed fromthe osculating Keplerian elements. These computed averages are then used for the
TTVFaster comparison simulation.Agreement thus also implies that the initial conditions used by first are a good approximation of the average valuesused by the latter. 4.5.
Libration width and proximity to MMR
The approximations utilized in
TTVFaster are based on the underlying assumptions that a system satisfies threeconditions: . small planet-to-star mass ratios (cid:0) m/m (cid:63) (cid:46) − (cid:1) . . small eccentricities. . adjacent planets are closeto, but are not in MMR (Agol & Deck 2016; Weiss et al. 2017; Tuchow et al. 2019).The first requirement, of small mass ratio, is verified by examination of the resulting best-fit values, which are allsmaller than 10 − for the systems presented.The second requirement, of small e , is met in systems with inferred non-zero eccentricities by performing the validitytest described above ( § δn max (their Eq. 8.74), the resonance width,is the maximum libration of the mean motion of a planet in resonance, which can be computed using the fitted valuesand choice of the appropriate resonance. For a given inner planet and resonance, n res is the exactly-resonant meanmotion of the outer planet. The observed mean motion n obs of that outer planet is then compared to this value andwe require that | n res − n obs | /δn max > i.e. that the planets are out of resonance. Solutions that do not meet thiscriterion are discarded as they violate the underlying assumptions of TTVFaster . We also verify the planets are nearMMR in the sense that | n res − n obs | /δn max is not large (typically of order 10 for the systems presented).4.6. PyDynamicaLC
We make available our dynamical light-curve generator code for public use. The code integrates several existingroutines –
TTVFaster , TTVFast (a modified C version which is also available) and
PyAstronomy ’s implementation of
EXOFAST , sewn together to generate a light-curve with planetary dynamics in one of three configurations:
Quasi-Circular : dynamics are simulated with
TTVFaster , a single transit shape is approximated as originatingfrom a circular orbit. Each transit event is only shifted in time by the corresponding TTV determined by the dynamicalsimulation.
Eccentric : dynamics are again simulated with
TTVFaster , the transit shapes are constant in time but with aprescribed eccentricity ( e = e ). Each transit event is shifted in time by the corresponding TTV determined by thedynamical simulation. Osculating : dynamics are simulated with
TTVFast , from which the instantaneous osculating Keplerian elements areextracted at every mid-transit instant. The shape and time of each transit is then determined from its correspondinginstantaneous elements, with the underlying assumption that Keplerian elements during each transit remain constant.While all configurations were validated against synthetic light-curves generated with
Mercury (Chambers 1999) and
EXOFAST , the work presented in this paper is based only on the quasi-circular configuration.We further provide a coupling of PyDynamicaLC to
MultiNest , allowing optimization of planetary masses andeccentricities following the methodology of this study, using the first two-modes of
PyDynamicaLC ( i.e. quasi-circularand eccentric) and analyze the resulting posterior distributions. Yoffe et al.
Parameter Kepler-36 b Kepler-36 c P [days] 13.849039 ± · − ± · − a/r (cid:63) ± ± t lin [days] 141.7185 ± ± · − b/r (cid:63) ± ± r p [ r ⊕ ] 1.582 ± ± m (cid:63) [ m (cid:12) ] 1.034 +0 . − . r (cid:63) [ r (cid:12) ] 1.629 +0 . − . Table 1.
Kepler-36 geometrical model best-fit average Keplerian parameters, in addition to stellar parameters, planetary radiiand their uncertainties reported in Fulton & Petigura (2018). σ e b [ m ⊕ ] 5.65 +0 . − . +0 . − . +0 . − . +0 . − . m c [ m ⊕ ] 9.76 +0 . − . +0 . − . +0 . − . +0 . − . e x,b +0 . − . +0 . − . +0 . − . < y,b -0.0031 +0 . − . -0.007 +0 . − . -0.034 +0 . − . < e x,c +0 . − . +0 . − . +0 . − . -∆ e y,c -0.0140 +0 . − . -0.018 +0 . − . -0.043 +0 . − . - − log( L ) 496641 495704 - Table 2.
Kepler-36 derived absolute planetary masses, absolute eccentricity components of the inner planet and eccentricitycomponent differences for the outer planets. Adopted values are in bold and results of Carter et al. (2012) are quoted in therightmost column for comparison (see also Fig. 1.) 5.
RESULTS5.1.
Kepler-36
We analyze Kepler-36 (KOI-277), a two-planet system close to a 6:7 MMR for which exceptionally precise massconstraints were previously inferred through TTV inversion by Carter et al. (2012): both masses were determined tohigh precision (to better than ∼ σ ). We reproduce these results in Table 2. We chose to analyze this system as aparticularly difficult test for our analysis to recover.The mean Keplerian parameters of the system are listed in Table 1. We performed three optimization rounds withdifferent eccentricity priors, as described in §
4. The best-fit values, their uncertainties and associated likelihoods, inaddition to the results of Carter et al. (2012) are listed in Table 2. We find the system to be nearly circular. In allthree priors, eccentricities remain consistent with zero to within ∼ σ .We plot our adopted absolute derived planetary masses versus their radii, and those of Carter et al. (2012) in Fig. 1.We find that our results are consistent with those of Carter et al. (2012) both in value and error estimates despitesome differences in the underlying data sets (due to reprocessing of Kepler data here, use of partial data by Carteret al. (2012) vs. full data set here, and inclusion of short-cadence data by Carter et al. (2012) vs. only long-cadencehere) - and obviously in the data analysis techniques. We also plot the resultant posterior distributions (Fig. 7) andfind that our parameterization of the eccentricity (as a single value e and a series of differences ∆ e i ) is indeed nearlyuncorrelated with the masses, though some correlations remain between other pairs of parameters.5.2. Kepler-79
We analyze Kepler-79 (KOI-152), a four-planet system close to the 2:1, 2:1 and 3:2 MMRs, respectively. The averageKeplerian parameters for the system are listed in Table 3. We performed three optimization rounds as described in § ∼ Simplified Photodynamical Model for Planetary Mass Determination g c m ! " g c m ! " g c m ! " cbb c Figure 1.
Kepler-36 mass-radius diagram displaying our adopted mass values (in green) and those of Carter et al. (2012) (inred), compared to three constant bulk density curves.Parameter Kepler-79 b Kepler-79 c Kepler-79 d Kepler-79 e P [days] 13.484542 ± · − ± · − ± · − ± · − a/r (cid:63) ± ± ± ± t lin [days] 136.6209 ± ± ± · − ± b/r (cid:63) ± ± ± ± r p [ r ⊕ ] 3.338 ± ± ± ± m (cid:63) [ m (cid:12) ] 1.244 +0 . − . r (cid:63) [ r (cid:12) ] 1.283 +0 . − . Table 3.
Kepler-79 geometrical model best-fit average Keplerian parameters, in addition to stellar parameters, planetary radiiand their uncertainties reported in Fulton & Petigura (2018). and smallest Kepler-79 b, and increase the mass of Kepler-79 d, moving it to a more physically plausible location onthe mass-radius diagram (Fig. 2). We also plot the resultant posterior distributions (Fig. 8) and find some, but notstrong, mass-eccentricity correlations.Similarly to Kepler-36 ( § − log( L ) of the different optimization runs are statistically insignificant. Wechoose the σ e = 0 .
05 optimization run as our adopted values based on its lowest − log( L ) value, though as mentioned,the improvement is small. 5.3. Kepler-450
We analyze Kepler-450 (KOI-279), a three-planet system close to the 2:1 MMR for both adjacent pairs. To ourknowledge, there exist no previous mass or eccentricity constraints in the literature for any of the planets in the system.The average Keplerian parameters for the system are listed in Table 5. We again performed three optimization runs.The fitting results, their uncertainties and associated likelihoods are listed in Table 6.In the case of Kepler-450, the results from the optimization runs assuming the two wider eccentricity prior distribu-tions indicate small but significantly non-zero eccentricities. Therefore, we perform our validity map analysis, varying∆ e x , ∆ e y near their best-fit values. This analysis shows that the fit lies, in fact, in a region showing 2-5 σ disagreementbetween the TTVFast and
TTVFaster models (see Fig. 4). This suggests that in these regions of parameter space
TTVFaster model may not be reliable . The σ e = 0 .
01 optimization run does, however, lie well within a favorableregion, with < σ disagreement between the two calculations (Fig. 4). Consequently, we choose this solution as ourpreferred values for this system. We plot our adopted planetary masses versus their radii in Fig. 3. Note, that whilethe masses Kepler-450 d and Kepler 450 b appear to be physically implausible – too dense in the case of the first andtoo light in the case of the latter, both have large error bars, and hence we do not attach significance to their mass Yoffe et al. σ e b [ m ⊕ ] 17.8 +5 . − . +4 . − . +4 . − . +7 . − . m c [ m ⊕ ] 12.4 +2 . − . +2 . − . +2 . − . +1 . − . m d [ m ⊕ ] 13.2 +2 . − . +1 . − . +2 . − . +2 . − . m e [ m ⊕ ] 6.85 +0 . − . +0 . − . +1 . − . +1 . − . e x,b -0.0125 +0 . − . -0.0148 +0 . − . -0.0173 +0 . − . +0 . − . e y,b +0 . − . +0 . − . -0.0007 +0 . − . +0 . − . ∆ e x,c -0.0145 +0 . − . -0.0171 +0 . − . -0.0195 +0 . − . -∆ e y,c -0.0065 +0 . − . -0.0123 +0 . − . -0.021 +0 . − . -∆ e x,d -0.0108 +0 . − . -0.014 +0 . − . -0.015 +0 . − . -∆ e y,d +0 . − . -0.004 +0 . − . -0.018 +0 . − . -∆ e x,e -0.0129 +0 . − . -0.015 +0 . − . -0.017 +0 . − . -∆ e y,e -0.0067 +0 . − . -0.014 +0 . − . -0.026 +0 . − . - − log( L ) 443838 443819 - Table 4.
Kepler-79 derived absolute planetary masses, absolute eccentricity components of the inner planet and eccentricitycomponent differences for the outer planets. Adopted values are boldfaced and results of Jontof-Hutter et al. (2013) are quotedin the rightmost column. g c m ! " g c m ! " c be d bce d . g c m ! " Figure 2.
Kepler-79 mass-radius diagram displaying our adopted mass values (in green) and those of Jontof-Hutter et al. (2013)(in red), compared to three constant bulk density curves.Parameter Kepler-450 d Kepler-450 c Kepler-450 b P [days] 7.5144279 ± · − ± · − ± · − a/r (cid:63) ± ± ± t lin [days] 136.17796 ± · − ± · − ± · − b/r (cid:63) ± ± ± r p [ r ⊕ ] 0.9408 ± ± ± m (cid:63) [ m (cid:12) ] 1.334 +0 . − . r (cid:63) [ r (cid:12) ] 1.6 +0 . − . Table 5.
Kepler-450 geometrical model best-fit average Keplerian parameters, in addition to stellar parameters, planetary radiiand their uncertainties reported in Fulton & Petigura (2018). constraints. We also plot the resultant posterior distributions (Fig. 9) and here too find little to no mass-eccentricitycorrelation using our parameterization.
Simplified Photodynamical Model for Planetary Mass Determination σ e d [ m ⊕ ] +9 . − . +6 . − . +6 . − . m c [ m ⊕ ] +3 . − . +30 − +35 − m b [ m ⊕ ] +11 . − . +16 − +14 − e x,d +0 . − . +0 . − . +0 . − . e y,d +0 . − . +0 . − . +0 . − . ∆ e x,c +0 . − . +0 . − . +0 . − . ∆ e y,c +0 . − . +0 . − . +0 . − . ∆ e x,b +0 . − . +0 . − . +0 . − . ∆ e y,b +0 . − . +0 . − . +0 . − . − log( L ) Table 6.
Kepler-450 derived absolute planetary masses, absolute eccentricity components of the inner planet and eccentricitycomponent differences for the outer planets. Adopted values are boldfaced. g c m ! " g c m ! " g c m ! " bc d Figure 3.
Kepler-450 mass-radius diagram displaying our adopted mass values compared to three constant bulk density curves.The best-fit masses of Kepler-450 b and Kepler-450 d are anomalous for their sizes, but are poorly constrained.
Kepler-595
We analyze Kepler-595 (KOI-547), which exhibits three sets of transit signals with pairs close to the 5:3 and 2:1MMRs. To our knowledge, there exist no previous mass or eccentricity constraints for any of the planets in the system,and while the outer planet is statistically validated (Morton et al. 2016) the two inner planets are still labeled ascandidates. The average Keplerian parameters for the system are listed in Table 7. The best-fit values of the threeoptimization runs, their uncertainties and associated likelihoods are listed in Table 8.Kepler-595 too requires the validity map test for results with significant non-zero eccentricities. Figure 5 showsthat the two lowest − log( L ) solutions, for σ e = 0 . , .
05, lie in regions of the parameter-space where
TTVFast and
TTVFaster models are not in agreement. The σ e = 0 .
01 optimization run, however, lies in a more favorable region. Wetherefore choose the values from this run as our adopted values for this system. Similarly to the case of Kepler-450 d,the best-fit mass of KOI-547.02 appears to be physically implausible, but it too has large error bars (consistent withzero at 2 σ ).We plot our preferred absolute derived planetary masses versus their radii in Fig. 6, and the resultant posteriordistributions (Fig. 10). Here the mass-eccentricity correlation is evident even in our eccentricity parameterization -possibly reflecting the dearth of true mass information. CONCLUSIONSWe implemented a simplified photodynamical model to infer planetary masses and eccentricities in multi-planetsystems. The approach couples dynamical output simulated with
TTVFaster (Agol & Deck 2016), a semi-analytic0
Yoffe et al. (outermost) (outermost)(outermost) (outermost) ! ! = 0.01! ! = 0.02! ! = 0.05 Figure 4.
Validity maps for Kepler-450. Best-fit solutions and their uncertainties are marked in green. Each row refers to anoptimization run with a different eccentricity prior scale-width ( σ e ). Here, the best-fit parameters are kept constant (Table 6)except the ∆ e x , ∆ e y values of each of the outer planets scanned in each column. Contours are of the χ ( § TTVFast and
TTVFaster generated light-curves. The values are converted to the equivalent number ofstandard deviations assuming a normal distribution with 3 n pl degrees of freedom. Best-fit values from each optimization runand their uncertainties are marked with a green cross.Parameter KOI-547.02 KOI-547.03 Kepler-595 b P [days] 7.3467925 ± · − ± · − ± · − a/r (cid:63) ± ± ± t lin [days] 134.5069 ± ± ± · − b/r (cid:63) ± ± ± r p [ r ⊕ ] 0.952 ± ± ± m (cid:63) [ m (cid:12) ] 0.93 ± r (cid:63) [ r (cid:12) ] 0.821 ± Table 7.
Kepler-595 geometrical model best-fit average Keplerian parameters, in addition to stellar parameters, planetary radiiand their uncertainties are from
NexSci . TTV model, with a Mandel & Agol (2002) analytical light-curve model. Optimization of planetary parameters wasdone by
MultiNest (Feroz et al. 2008), a multi-modal Bayesian inference algorithm, applied using a unique eccentricityparametrization: the eccentricity differences between adjacent planets as our primary degrees of freedom, which enablesthe optimization to be significantly less prone to correlations. The median CPU-time that was required for the systemsanalyzed is ∼ TTVFaster , the combination of the choice ofjump parameters and the fitting the light-curve and not the transit timings improved the mass sensitivity here, whichallowed the analysis of a sample of
Kepler multi-planet systems, and in particular planets that were thus far too smallto have well-constrained transit timings. Finally, we quality-check our best-fit values by ensuring that they are at a
Simplified Photodynamical Model for Planetary Mass Determination σ e . [ m ⊕ ] +53 − +75 − +42 − m . [ m ⊕ ] +1 . − . +1 . − . +0 . − . m b [ m ⊕ ] +7 . − . +9 . − . +8 . − . e x, . -0.0108 +0 . − . -0.021 +0 . − . -0.060 +0 . − . e y, . +0 . − . +0 . − . +0 . − . ∆ e x, . -0.0219 +0 . − . -0.032 +0 . − . -0.042 +0 . − . ∆ e y, . +0 . − . +0 . − . +0 . − . ∆ e x,b -0.007 +0 . − . -0.015 +0 . − . -0.043 +0 . − . ∆ e y,b +0 . − . +0 . − . +0 . − . − log( L ) Table 8.
Kepler-595 derived absolute planetary masses, absolute eccentricity components of the inner planet and eccentricitycomponent differences for the outer planets. Adopted values are boldfaced. ! ! = 0.01! ! = 0.02! ! = 0.05 Figure 5.
Validity maps for Kepler-595, similar to Fig. 4 region of parameter space where
TTVFaster is consistent with
TTVFast , and that the suggested
TTVFaster solution isclose to, but not in, resonance.We were able to provide significant planetary mass (and eccentricity) constraints in several planetary systems: wedecreased the uncertainty of Kepler-79 b by a factor of ∼
2, and our derived mass values of Kepler-79 d suggestsmore physically plausible bulk density. We also provide the first mass constraints to Kepler-595 b, Kepler-450 c andKOI-547.03 - a Neptune-sized , sub-Neptune and an Earth-sized planets respectively. The mass of Kepler-450 b ismarginally detected to 2 . σ . The simplified photodynamical code used in this study is made publicly available .Looking forward, we plan to use this novel and now validated technique to analyze all multi-transiting systems fromthe Kepler and
TESS data sets with known TTVs. https://github.com/AstroGidi/PyDynamicaLC Yoffe et al. c g c m ! " g c m ! " g c m ! " KOI-547.02b
Figure 6.
Kepler-595 mass-radius diagram displaying our adopted mass values compared to three constant bulk density curves.Error-bars denote 1 σ uncertainty in m and r p . Note the mass of KOI-547.02 is poorly constrained. ACKNOWLEDGEMENTSWe wish to thank Dr. Trifon Trifonov and Yair Judkovsky for helpful discussions, and Prof. Eric Agol for adviceon
TTVFaster and
TTVFast . This study was supported by the Helen Kimmel Center for Planetary Sciences and theMinerva Center for Life Under Extreme Planetary Conditions
Agol E., Deck K., 2016, Astrophys. J., 818, 177Agol E., Steffen J., Sari R., Clarkson W., 2005, Month.Not. Royal Acad. Soc., 359, 567Buchner J., et al., 2014, Astro. Astrophys., 564, A125Carter J. A., et al., 2012, Science, 337, 556Chambers J. E., 1999, Month. Not. Royal Acad. Soc., 304,793Deck K. M., Agol E., Holman M. J., Nesvorn´y D., 2014,Astrophys. J., 787, 132Eastman J., Gaudi B. S., Agol E., 2013, Publications of theAstronomical Society of the Pacific, 125, 83Feroz F., Hobson M. P., Bridges M., 2008, Mon. Not. Roy.Astron. Soc. 398: 1601-1614,2009Fulton B. J., Petigura E. A., 2018, Astronomical Journal,156, 264Hadden S., Lithwick Y., 2014, Astrophys. J., 787, 80Hadden S., Lithwick Y., 2016, Astrophys. J., 828, 44Hadden S., Lithwick Y., 2017, Astronomical Journal, 154, 5Holczer T., et al., 2016, Astrophys. J. Supplement Series,225, 9Holman M. J., 2005, Science, 307, 1288Jontof-Hutter D., Lissauer J. J., Rowe J. F., FabryckyD. C., 2013, Astrophys. J.Jontof-Hutter D., et al., 2016, Astrophys. J., 820, 39Linial I., Gilbaum S., Sari R., 2018, Astrophys. J., 860, 16 Lithwick Y., Xie J., Wu Y., 2012, Astrophys. J., 761, 122MacDonald M. G., et al., 2016, Astrophys. J., 152, 105Mandel K., Agol E., 2002, Astrophys. J., 580, L171Moorhead A. V., et al., 2011, Astrophys. J. SupplementSeries, 197, 1Morton T. D., Bryson S. T., Coughlin J. L., Rowe J. F.,Ravichandran G., Petigura E. A., Haas M. R., BatalhaN. M., 2016, Astrophys. J., 822, 86Murray C. D., Dermott S. F., 2000, Solar SystemDynamics. Cambridge University Press,doi:10.1017/cbo9781139174817,https://doi.org/10.1017/cbo9781139174817Ofir A., Dreizler S., 2013, Astro. Astrophys., 555, A58Ofir A., Xie J.-W., Jiang C.-F., Sari R., Aharonson O.,2018, Astrophys. J., Supp., 234, 9Payne M. J., Ford E. B., Veras D., 2010, Astrophys. J., 712,L86Skilling J., 2004, in AIP Conference Proceedings. AIP,doi:10.1063/1.1835238,https://doi.org/10.1063/1.1835238Tuchow N. W., Ford E. B., Papamarkou T., Lindo A.,2019, Month. Not. Royal Acad. Soc., 484, 3772Weiss L. M., et al., 2017, Astronomical Journal, 153, 265Winn J. N., et al., 2011, Astrophys. J., 737, L18Xie J.-W., et al., 2016, Proceedings of the NationalAcademy of Sciences, 113, 11431
Simplified Photodynamical Model for Planetary Mass Determination Figure 7.
Kepler-36 correlation plot. The color bar is the relative likelihood of each point, relative to the best-fit point in theposterior sample (Relative Likelihood ∝ exp( − ∆ log L )). APPENDIX A. CORRELATION PLOTSIn the following section we provide a plot for each system considered, in which pairwise correlations are displayed.4
Yoffe et al.
Figure 8.
Kepler-79 correlation plot, similarly to Fig. 7.
Figure 9.
Kepler-450 correlation plot, similarly to Fig. 7.