A Simplified Photodynamical Model for Planetary Mass Determination in Low-eccentricity Multitransiting Systems

Transit timing variations (TTVs) caused by the gravitational interactions in multiplanet systems (Agol et al. 2005; Holman 2005) may be modeled to determine the governing system parameters, particularly masses. This is especially valuable in cases where no additional observations beyond photometry, such as radial velocity data, are available. 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 techniques from the literature.

Global model: The standard technique for the detection of TTVs follows two steps: first, one optimizes an empirical light-curve model using common shape parameters for all of the transit events (in the usual Mandel & Agol 2002 framework these are $a/{r}_{\star }$, $b/{r}_{\star }$, and ${r}_{p}/{r}_{\star }$, the planetary semimajor axis, impact parameter, and radius, each normalized by the stellar radius) but assigns an individual time for each transit event. Second, the list of best-fitting transit times is modeled with a dynamical description. There are shortcomings to this approach: the transits must be deep enough to be individually significant, they must be long enough to have more sufficient points that determine their timing, and the number of fitting parameters grows with the number of transit, greatly diminishing the sensitivity of this approach. To rectify these shortcomings Ofir et al. (2018) recently introduced a more sensitive technique for the detection of TTVs: one assumes sinusoidal shape for the TTVs, and a global optimization allows determination of a single set of parameters describing the full light-curve variations. This approach is effective for small transit depths or durations (even in cases where individual transits cannot be seen), and has significantly fewer degrees of freedom—enhancing its sensitivity, as discussed next.

Parsimonious model: The problem of inverting the observed TTVs for the component masses is nontrivial. The degeneracy between the effects of masses and orbital eccentricities has been appreciated in early works on the subject (Agol et al. 2005; Holman 2005) and a solution based on a full n-body integration is highly multidimensional, even more degenerate than just explained. It is thus difficult to initialize and optimize, and is computationally expensive or even prohibitive. More recently, it was shown that the pattern of the TTVs does not depend strongly on the individual eccentricities but rather only on the eccentricity-vector difference between adjacent planetary pairs (Hadden & Lithwick 2016; Jontof-Hutter et al. 2016; Linial et al. 2018). Fitting for the eccentricity differences allows us to capture most of the TTV variability using two fewer degrees of freedom.

We use TTVFaster (Agol & Deck 2016) as our dynamical tool of choice. TTVFaster is a semianalytic model accurate to first order in eccentricity, which approximates TTVs using a series expansion. An important simplification in using this tool arises from its reliance only on average Keplerian elements, which are well determined from previous processing 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 \omega ,e\sin \omega$), where ω is the argument of periapse and e is the eccentricity magnitude. Consequently, there are fewer correlations among parameters, producing better determined values, and initial guesses are more easily selected. Additionally, TTVFaster offers a considerably faster 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 multitransiting Kepler systems.

An important disadvantage of using TTVFaster is that its underlying approximations are inaccurate for systems with high eccentricities. Therefore, we require that our model pass a series of tests to verify the validity of our model in the relevant region of the parameter space. We also note that performing fits in flux space (the light curve) does not 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 that is both global and minimal in degrees of freedom to achieve better precision and sensitivity than previous studies, and in a uniform fashion to multiple relevant Kepler targets. We also use the recently updated stellar parameters based on Gaia DR2 (Fulton & Petigura 2018), which allows us to better derive the planetary physical properties. The paper is divided as follows: Section 2 describes light-curve processing, in Section 3 we describe our photodynamical model, in Section 4 we describe our fitting procedure, in Section 5 we present and discuss the results of the modeling, and in Section 6 we conclude.

The sample of Kepler systems we processed is all the multitransiting systems of which at least one component was identified 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. The gravitational 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 on Ofir & Dreizler (2013). All the planets in the system were fitted simultaneously using a common value $a/{r}_{\star }$, where a is the semimajor axis and ${r}_{\star }$ is the stellar radius. This value was scaled between planets using Kepler's Third Law. TTVs, if they exist, were accounted for by a phase shift of each event. This is a circular-orbit approximation of the TTVs. The detrending curve itself was the better of either a Savitzky–Golay filter or a cosine filter—applied to each section of continuous data independently. Each of the filters was applied several times—as many as there are known transiting planets in the system—since each transit has a unique duration. In every continuous section the transit with the longest duration appearing on that section determined the critical timescale for filtering, and each of the filters (Savitzky–Golay or cosine) used a timescale no shorter than three times the critical timescale on that section.

In order to account for all TTVs we started by using the Holczer et al. (2016) catalog of TTVs. In cases where measuring individual timings is difficult (shallow or short transits) individual planets were designated as not TTV-bearing. In these cases a periodic model is also likely a good approximation of the true signal (both would be difficult to detect in the first place had they exhibited large amplitude TTVs). The individual transit times were manually checked to verify that no transits were missed (e.g., Holczer et al. 2016 do not report partial transits or transit close to data gaps), that no transits were unnecessarily added (some reported transit timings are at times when there is a gap in the data, for example), and in certain cases whole planets were added or removed from systems based on literature later than Holczer et al. (2016). We note that even when starting from a comprehensive catalog like Holczer et al. (2016), minimizing errors at this stage is labor intensive. For example, by the definition of TTVs their times are not exactly known and sometimes transits are "unexpectedly" visible in the data while the linear ephemeris suggest that the transit should fall in a data gap. Or, when transits of two similarly sized planets that both exhibit TTVs happen to be close in time, or even overlap, correctly assigning a specific transit to a particular planet is sometimes not immediately obvious. We note this since the need for the most accurate detrending of the light curve stems from the goal of detecting low-amplitude TTVs of shallow transits—these are inherently faint signals, and any contamination may wash away the sought-after signal altogether.

The end result of this stage is a normalized, detrended light curve as well as a set of parameters describing the linear ephemeris for each planet: the mean period P, a mean reference epoch t_lin (i.e., the constant term in the linear fit to the times of transit), $b/{r}_{\star }$, ${r}_{p}/{r}_{\star }$, and $a/{r}_{\star }$ for the innermost planet. Note that the empirically fitted individual times of midtransit were not used in the dynamical modeling below.

For each multiplanet system, we compute a photodynamical model in three steps: first, we simulate a TTV signal with TTVFaster. The required input for each planet i are the average Keplerian elements, in addition to the planetary masses and the first time of midtransit for a linearly approximated orbit (_Pi , i_i , ${e}_{x,i}$, ${e}_{y,i}$, ${t}_{\mathrm{lin},i}$, and m_i ), where i refers to the index of the planet, _Pi is the orbital period, i is the orbital inclination that can be converted from b with Equation (7) in Winn et al. (2011), ${e}_{x,i}$ and ${e}_{y,i}$ are the perpendicular eccentricity components, respectively, ${t}_{\mathrm{lin},i}$ is the mean reference epoch and m_i is the mass. Note that throughout this study the Keplerian parameters refer to their time averages over Kepler long-cadence observations (as opposed to their osculating values).

Excluding masses and eccentricities, these parameters are well constrained during the data processing stage (Section 2). Thus remain $3{n}_{\mathrm{pl}}$ degrees of freedom for the dynamical fit. The output obtained from TTVFaster is a vector of times of midtransit, from which TTVs may be deduced.

Each individual transit event is then shifted from a strict periodic regime and generates the individual-event transit light-curve Mandel–Agol model (Mandel & Agol 2002), using PyAstronomy's implementation of EXOFAST ⁴ (Eastman et al. 2013). Assuming small orbital eccentricity allows us to further simplify our model and approximate the shape of the transit curve to be that of an arc of a circular orbit (referred to as a quasi-circular approximation). An additional assumption throughout this study is that all the TTVs can be attributed to mutual perturbations among adjacent pairs of known transiting planets. While TTVs can be sensitive to mutual inclinations (Payne et al. 2010), Kepler multiplanet systems are usually characterized by small values of these (Xie et al. 2016), such that this effect does not make a notable contribution to TTVs of planets near low-order MMR (Payne et al. 2010).

Finally, we perform nonlinear optimization of the dynamical parameters and estimate their uncertainties using MultiNest (Feroz et al. 2008), a multimodal 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 the efficient calculation of the Bayesian evidence. Using MultiNest allows for the exploration of regions of interest in the parameter space with a reduced risk of being trapped in local likelihood maxima, which is especially important due to the known mass-eccentricity degeneracy (Lithwick et al. 2012).

For any multiplanet model, we fit $3{n}_{\mathrm{pl}}$ free parameters and also use $4{n}_{\mathrm{pl}}+5$ input constants. The latter are the linear ephemeris (see Section 2) and stellar parameters—${m}_{\star }$ and ${r}_{\star }$, two stellar limb darkening coefficients, and a of the innermost planet. Stellar parameters are adopted from Fulton & Petigura (2018).

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 adjacent planetary pair, rather than on the individual eccentricities. Attempting to fit for both individual eccentricities yields highly degenerate results, e.g., Jontof-Hutter et al. (2016, their Figures 3, 5, 7, 9, 11, 15, 17, and 19). Fitting for the eccentricity difference, however, results in considerably weaker mass-eccentricity degeneracy (see correlation plots in the Appendix). Therefore, we fit the absolute orbital eccentricity (${e}_{x,0}$, ${e}_{y,0}$) for only the innermost planet (expecting it to remain poorly constrained), and for the remaining planets we fit for the eccentricity differences (${\rm{\Delta }}{e}_{x,i+1}={e}_{x,i+1}-{e}_{x,i}$, ${\rm{\Delta }}{e}_{y,i+1}={e}_{y,i+1}-{e}_{y,i}$).

4.1. Statistical Model

We evaluated the goodness of fit to the multiplanet light curve for each simulated model with parameter set Θ using the log of Student's t-distribution function with 4 degrees of freedom, whose log is

$$\begin{eqnarray}&&\mathrm{log}({\mathscr{L}})=-\displaystyle \sum _{i=1}^{n}\mathrm{log}\left(2\sqrt{2}{\sigma }_{i}\right)-\displaystyle \frac{5}{2}\displaystyle \sum _{i=1}^{n}\mathrm{log}\left(1+\displaystyle \frac{{\left({y}_{i,{\rm{\Theta }}}-{x}_{i}\right)}^{2}}{4{\sigma }_{i}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

where there are n observed photometric measurements x_i with respective uncertainties ${\sigma }_{i}$, and y_i is the simulated multiplanet 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 two normalized and symmetric statistical distributions, a Gaussian, and Student's t-distribution with an optimized number of 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. In our analysis we find that similarly to Jontof-Hutter et al. (2016) and MacDonald et al. (2016), the wings of the distribution contain more outliers than would be expected by Gaussian uncertainties, such that the Student's t-distribution scored better in the Kolmogorov–Smirnov test.

4.2. Priors

4.3. Posteriors

4.4. Validity Map

To validate our fits for systems with significant nonzero eccentricities, we compare the output of TTVFaster to that of TTVFast—a symplectic n-body integrator with a Keplerian interpolator for the calculation of transit timings, itself verified 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 ${\rm{\Delta }}{e}_{x}$ and ${\rm{\Delta }}{e}_{y}$ for every planet pair, keeping the remaining parameters constant at their best-fit values (producing ${n}_{\mathrm{pl}}-1$ validity grids). At each point in the validity map the squared difference between TTVFast and TTVFaster normalized by the data's error estimation are used to construct a ${\chi }^{2}$-like statistic:

$$\begin{eqnarray}&&{\chi }_{\mathrm{validity}}^{2}=\displaystyle \sum _{i=1}^{n}{\left(\displaystyle \frac{{y}_{i,\mathrm{TTVFast}}-{y}_{i,\mathrm{TTVFaster}}}{{\sigma }_{i}}\right)}^{2}.\end{eqnarray} \tag{ 2 }$$

This allows evaluation of the region in parameter space in which TTVFast and TTVFaster are consistent to within the measurement error. A difficulty in performing the analysis lies in the fact that only the average Keplerian parameters are known (with uncertainty), whereas TTVFast requires the instantaneous values at the beginning of integration. We therefore use the average Keplerian elements as initial conditions, from which revised average values are computed from the osculating Keplerian elements. These computed averages are then used for the TTVFaster comparison simulation. Agreement thus also implies that the initial conditions used by the first are a good approximation of the average values used by the latter.

4.5. Libration Width and Proximity to MMR

4.6.  PyDynamicaLC

5.1. Kepler-36

We analyze Kepler-36 (KOI-277), a two-planet system close to a 6:7 MMR for which exceptionally precise mass constraints were previously inferred through TTV inversion by Carter et al. (2012): both masses were determined to high precision (to better than $\sim 14\sigma$). We reproduce these results in Table 2. We chose to analyze this system as a particularly difficult test for our analysis to recover.

Table 1. Kepler-36 Geometrical Model Best-fit Average Keplerian Parameters, in Addition to Stellar Parameters, Planetary Radii, and Their Uncertainties Reported in Fulton & Petigura (2018)

Parameter	Kepler-36 b	Kepler-36 c
P [days]	13.849039 ± 7.8 × 10−5	16.2318713 ± 9.4 × 10−6
$a/{r}_{\star }$	14.76 ± 0.19	16.41 ± 0.22
${t}_{\mathrm{lin}}$ [days]	141.7185 ± 0.0047	171.67304 ± 4.5 × 10−4
$b/{r}_{\star }$	0.359 ± 0.041	0.244 ± 0.046
rp $[{r}_{\oplus }]$	1.582 ± 0.015	3.692 ± 0.010
${m}_{\star }$ $[{m}_{\odot }]$	1.034${}_{-0.012}^{+0.022}$
${r}_{\star }$ $[{r}_{\odot }]$	1.629+0.020 −0.019

Download table as:  ASCIITypeset image

Table 2. Kepler-36 Derived Absolute Planetary Masses, Absolute Eccentricity Components of the Inner Planet and Eccentricity Component Differences for the Outer Planets

${\sigma }_{e}$	0.01	0.02	0.05	Carter+12
mb $[{m}_{\oplus }]$	5.65+0.35 −0.48	5.61+0.36 −0.46	5.56+0.41 −0.45	4.45+0.33 −0.27
mc $[{m}_{\oplus }]$	9.76+0.66 −0.91	9.78+0.70 −0.88	9.61+0.79 −0.89	8.08+0.6 −0.46
${e}_{x,b}$	0.0035+0.0076 −0.0090	0.011+0.014 −0.017	0.045+0.023 −0.043	<0.039
${e}_{y,b}$	−0.0031+0.0092 −0.0082	−0.007+0.018 −0.014	−0.034+0.045 −0.033	<0.033
${\rm{\Delta }}{e}_{x,c}$	0.0117+0.0069 −0.0082	0.019+0.013 −0.016	0.045+0.023 −0.043	⋯
${\rm{\Delta }}{e}_{y,c}$	−0.0140+0.0084 −0.0076	−0.018+0.016 −0.013	−0.043+0.041 −0.030	⋯
$-\mathrm{log}({\mathscr{L}})$	496641	495704	491722	⋯

Note. Adopted values are in bold and results of Carter et al. (2012) are quoted in the rightmost column for comparison (see also Figure 1).

Download table as:  ASCIITypeset image

The mean Keplerian parameters of the system are listed in Table 1. We performed three optimization rounds with different eccentricity priors, as described in Section 4. The best-fit values, their uncertainties and associated likelihoods, in addition to the results of Carter et al. (2012) are listed in Table 2. We find the system to be nearly circular. In all three priors, eccentricities remain consistent with zero to within $\sim 2\sigma$.

We plot our adopted absolute derived planetary masses versus their radii, and those of Carter et al. (2012) in Figure 1. We find that our results are consistent with those of Carter et al. (2012) both in value and error estimates despite some differences in the underlying data sets (due to reprocessing of Kepler data here, use of partial data by Carter et al. 2012 versus full data set here, and inclusion of short-cadence data by Carter et al. 2012 versus only long-cadence here)—and obviously in the data analysis techniques. We also plot the resultant posterior distributions (see the Appendix) and find that our parameterization of the eccentricity (as a single value e₀ and a series of differences ${\rm{\Delta }}{e}_{i}$) is indeed nearly uncorrelated with the masses, though some correlations remain between other pairs of parameters.

Zoom In Zoom Out Reset image size

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 average Keplerian parameters for the system are listed in Table 3. We performed three optimization rounds as described in Section 4.2. The adopted values, their uncertainties, and associated likelihoods, in addition to the adopted values of Jontof-Hutter et al. (2014) are listed in in Table 4.

Table 3. Kepler-79 Geometrical Model Best-fit Average Keplerian Parameters, in Addition to Stellar Parameters, Planetary Radii and their Uncertainties Reported in Fulton & Petigura (2018)

Parameter	Kepler-79 b	Kepler-79 c	Kepler-79 d	Kepler-79 e
P [days]	13.484542 ± 2.1 × 10−5	27.402300 ± 4.2 × 10−5	52.090767 ± 3.2 × 10−5	81.063633 ± 4.4 × 10−4
$a/{r}_{\star }$	19.37 ± 0.25	31.08 ± 0.40	47.69 ± 0.61	64.05 ± 0.82
${t}_{\mathrm{lin}}$ [days]	136.6209 ± 0.0012	133.6276 ± 0.0013	158.74645 ± 5.4 × 10−4	139.5377 ± 0.0045
$b/{r}_{\star }$	0.360 ± 0.035	0.057 ± 0.077	0.043 ± 0.076	0.934 ± 0.016
rp $[{r}_{\oplus }]$	3.338 ± 0.028	3.553 ± 0.027	6.912 ± 0.014	3.414 ± 0.129
${m}_{\star }$ $[{m}_{\odot }]$	1.244+0.027 −0.042
${r}_{\star }$ $[{r}_{\odot }]$	1.283+0.029 −0.015

Download table as:  ASCIITypeset image

Table 4. Kepler-79 Derived Absolute Planetary Masses, Absolute Eccentricity Components of the Inner Planet, and Eccentricity Component Differences for the Outer Planets

${\sigma }_{e}$	0.01	0.02	0.05	JH+14
mb $[{m}_{\oplus }]$	17.8+5.4 −3.8	15.0+4.5 −3.6	12.5+4.5 −3.6	10.9+7.4 −6.0
mc $[{m}_{\oplus }]$	12.4+2.1 −1.8	11.0+2.1 −1.9	9.5+2.3 −2.1	5.9+1.9 −2.3
md $[{m}_{\oplus }]$	13.2+2.1 −1.9	12.3+1.9 −1.9	11.3+2.2 −2.2	6.0+2.1 −1.6
me $[{m}_{\oplus }]$	6.85+0.91 −0.88	6.53+0.86 −0.87	6.3+1.0 −1.0	4.1+1.4 −1.1
${e}_{x,b}$	−0.0125+0.0026 −0.0030	−0.0148+0.0033 −0.0039	−0.0173+0.0046 −0.0061	${0.032}_{-0.032}^{+0.036}$
${e}_{y,b}$	0.0014+0.0020 −0.0019	0.0007+0.0026 −0.0025	−0.0007+0.0038 −0.0038	${0.032}_{-0.029}^{+0.059}$
${\rm{\Delta }}{e}_{x,c}$	−0.0145+0.0040 −0.0044	−0.0171+0.0051 −0.0056	−0.0195+0.0072 −0.0081	⋯
${\rm{\Delta }}{e}_{y,c}$	−0.0065+0.0051 −0.0053	−0.0123+0.0072 −0.0079	−0.021+0.011 −0.015	⋯
${\rm{\Delta }}{e}_{x,d}$	−0.0108+0.0075 −0.0077	−0.014+0.012 −0.012	−0.015+0.021 −0.021	⋯
${\rm{\Delta }}{e}_{y,d}$	0.0028+0.0076 −0.0074	−0.004+0.013 −0.012	−0.018+0.021 −0.023	⋯
${\rm{\Delta }}{e}_{x,e}$	−0.0129+0.0062 −0.0066	−0.015+0.010 −0.010	−0.017+0.017 −0.017	⋯
${\rm{\Delta }}{e}_{y,e}$	−0.0067+0.0063 −0.0063	−0.014+0.011 −0.011	−0.026+0.018 −0.020	⋯
$-\mathrm{log}({\mathscr{L}})$	443838	443819	443806	⋯

Note. Adopted values are boldfaced and results of Jontof-Hutter et al. (2014) are quoted in the rightmost column.

Download table as:  ASCIITypeset image

For this system, mass constraints were previously inferred through TTV inversion by Jontof-Hutter et al. (2014). In their analysis, Jontof-Hutter et al. (2014) find that the four planets have low bulk densities, in particular Kepler-79 d. In this case, our approach allows us to decrease the mass uncertainty of Kepler-79 b by a factor of ∼2 for the innermost and smallest Kepler-79 b, and increase the mass of Kepler-79 d, moving it to a more physically plausible location on the mass–radius diagram (Figure 2). We also plot the resultant posterior distributions (see the Appendix) and find some, but not strong, mass-eccentricity correlations.

Zoom In Zoom Out Reset image size

Similarly to Kepler-36 (Section 5.1), we find all planets in this system have eccentricities that are consistent with zero. For 12 degrees of freedom, the variations in $-\mathrm{log}({\mathscr{L}})$ of the different optimization runs are statistically insignificant. We choose the ${\sigma }_{e}=0.05$ optimization run as our adopted values based on its lowest $-\mathrm{log}({\mathscr{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 our knowledge, 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. We plot our adopted planetary masses versus their radii in Figure 3.

Table 5. Kepler-450 Geometrical Model Best-fit Average Keplerian Parameters, in Addition to Stellar Parameters, Planetary Radii, and Their Uncertainties Reported in Fulton & Petigura (2018)

Parameter	Kepler-450 d	Kepler-450 c	Kepler-450 b
P [days]	7.5144279 ± 5.5 × 10−6	15.4131395 ± 2.9 × 10−6	28.4548844 ± 1.0 × 10−6
$a/{r}_{\star }$	10.859 ± 0.012	17.531 ± 0.020	26.383 ± 0.031
${t}_{\mathrm{lin}}$ [days]	136.17796 ± 5.0 × 10−4	136.94162 ± 1.1 × 10−4	176.705831 ± 2.1 × 10−5
$b/{r}_{\star }$	0.5214 ± 0.0033	0.3007 ± 0.0027	0.3038 ± 0.0032
rp $[{r}_{\oplus }]$	0.9408 ± 0.0018	2.5958 ± 0.0017	6.0834 ± 0.0022
${m}_{\star }$ $[{m}_{\odot }]$	1.334 ${}_{-0.022}^{+0.023}$
${r}_{\star }$ $[{r}_{\odot }]$	1.6 ${}_{-0.008}^{+0.028}$

Download table as:  ASCIITypeset image

Table 6. Kepler-450 Derived Absolute Planetary Masses, Absolute Eccentricity Components of the Inner Planet and Eccentricity Component Differences for the Outer Planets

${\sigma }_{e}$	0.01	0.02	0.05
md $[{m}_{\oplus }]$	17.6 ${}_{-6.7}^{+9.5}$	${11.4}_{-5.1}^{+6.8}$	${7.5}_{-3.5}^{+6.8}$
mc $[{m}_{\oplus }]$	12.5 ${}_{-2.6}^{+3.2}$	${55}_{-31}^{+30}$	${33}_{-24}^{+35}$
mb $[{m}_{\oplus }]$	19.4+11.1 −6.8	30+16 −13	25+14 −11
ex, d	0.0079+0.0049 −0.0040	0.0037+0.0054 −0.0031	0.0058+0.0158 −0.0050
${e}_{y,d}$	0.0128+0.0051 −0.0061	0.0303+0.0084 −0.0069	0.039+0.022 −0.014
${\rm{\Delta }}{e}_{x,c}$	0.0194+0.0077 −0.0073	0.0059+0.0083 −0.0031	0.0094+0.0254 −0.0059
${\rm{\Delta }}{e}_{y,c}$	0.0108+0.0085 −0.0094	0.0451+0.0034 −0.0104	0.0422+0.0055 −0.0283
${\rm{\Delta }}{e}_{x,b}$	0.0071+0.0097 −0.0096	0.0046+0.0082 −0.0059	0.0076+0.0239 −0.0093
${\rm{\Delta }}{e}_{y,b}$	0.0245+0.0095 −0.0134	0.0541+0.0057 −0.0075	0.056+0.011 −0.011
$-\mathrm{log}({\mathscr{L}})$	500095	500099	500098

Note. Adopted values are boldfaced.

Download table as:  ASCIITypeset image

In the case of Kepler-450, the results from the optimization runs assuming the two wider eccentricity prior distributions indicate small but significantly nonzero eccentricities. Therefore, we perform our validity map analysis, varying ${\rm{\Delta }}{e}_{x},{\rm{\Delta }}{e}_{y}$ near their best-fit values. This analysis shows that the fit lies, in fact, in a region showing 2σ–5σ disagreement between the TTVFast and TTVFaster models (see Figure 4). This suggests that in these regions of parameter space the TTVFaster model may not be reliable. The ${\sigma }_{e}=0.01$ optimization run does, however, lie well within a favorable region, with $\lt 1\sigma$ disagreement between the two calculations (Figure 4). Consequently, we choose this solution as our preferred values for this system. Note that while the masses of Kepler-450 d and Kepler-450 b appear to be physically implausible—too dense in the case of the first and too light in the case of the latter, both have large error bars, and hence we do not attach significance to their mass constraints. We also plot the resultant posterior distributions (see the Appendix) and here too find little to no mass-eccentricity correlation using our parameterization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.4. Kepler-595

We analyze Kepler-595 (KOI-547), which exhibits three sets of transit signals with pairs close to the 5:3 and 2:1 MMRs. 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 as candidates. The average Keplerian parameters for the system are listed in Table 7. The best-fit values of the three optimization runs, their uncertainties, and associated likelihoods are listed in Table 8.

Table 7. Kepler-595 Geometrical Model Best-fit Average Keplerian Parameters, in Addition to Stellar Parameters, Planetary Radii, and Their Uncertainties Are from NexSci

Parameter	KOI-547.02	KOI-547.03	Kepler-595 b
P [days]	7.3467925 ± 5.9 × 10−6	12.38602 ± 5.0 × 10−6	25.3029092 ± 6.0 × 10−6
$a/{r}_{\star }$	19.613 ± 0.035	27.783 ± 0.049	44.730 ± 0.080
${t}_{\mathrm{lin}}$ [days]	134.5069 ± 0.0052	132.8192 ± 0.0037	188.06096 ± 1.0 × 10−4
$b/{r}_{\star }$	0.593 ± 0.052	0.710 ± 0.063	0.206 ± 0.051
rp $[{r}_{\oplus }]$	0.952 ± 0.015	1.009 ± 0.024	3.708 ± 0.010
${m}_{\star }$ $[{m}_{\odot }]$	0.93 ± 0.059
${r}_{\star }$ $[{r}_{\odot }]$	0.821 ± 0.242

Download table as:  ASCIITypeset image

Table 8. Kepler-595 Derived Absolute Planetary Masses, Absolute Eccentricity Components of the Inner Planet, and Eccentricity Component Differences for the Outer Planets

${\sigma }_{e}$	0.01	0.02	0.05
${m}_{547.02}$ $[{m}_{\oplus }]$	29 ${}_{-19}^{+53}$	${35}_{-25}^{+75}$	126+42 −41
${m}_{547.03}$ $[{m}_{\oplus }]$	3.3 ${}_{-1.0}^{+1.7}$	2.34+1.57 −0.85	2.59+0.94 −0.81
mb $[{m}_{\oplus }]$	17.4+7.1 −3.8	14.9+9.4 −5.0	23.1+8.3 −7.8
${e}_{x,547.02}$	−0.0108+0.0074 −0.0063	−0.021+0.014 −0.012	−0.060+0.037 −0.034
${e}_{y,547.02}$	0.0055+0.0070 −0.0080	0.017+0.014 −0.015	0.094+0.031 −0.037
${\rm{\Delta }}{e}_{x,547.03}$	−0.0219+0.0088 −0.0082	−0.032+0.014 −0.016	−0.042+0.027 −0.025
${\rm{\Delta }}{e}_{y,547.03}$	0.0027+0.0054 −0.0052	0.0078+0.0098 −0.0090	0.028+0.022 −0.022
${\rm{\Delta }}{e}_{x,b}$	−0.007+0.011 −0.011	−0.015+0.023 −0.022	−0.043+0.070 −0.067
${\rm{\Delta }}{e}_{y,b}$	0.007+0.012 −0.012	0.018+0.021 −0.023	0.079+0.051 −0.052
$-\mathrm{log}({\mathscr{L}})$	407075	407072	406992

Note. Adopted values are boldfaced.

Download table as:  ASCIITypeset image

Kepler-595 also requires the validity map test for results with significant nonzero eccentricities. Figure 5 shows that the two lowest $-\mathrm{log}({\mathscr{L}})$ solutions, for ${\sigma }_{e}=0.02,0.05$, lie in regions of the parameter space, where TTVFast and TTVFaster models are not in agreement. The ${\sigma }_{e}=0.01$ optimization run, however, lies in a more favorable region. We therefore 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 with zero at 2σ).

Zoom In Zoom Out Reset image size

We plot our preferred absolute derived planetary masses versus their radii in Figure 6, and the resultant posterior distributions (see the Appendix). Here the mass-eccentricity correlation is evident even in our eccentricity parameterization—possibly reflecting the dearth of true mass information.

Zoom In Zoom Out Reset image size

We implemented a simplified photodynamical model to infer planetary masses and eccentricities in multiplanet systems. The approach couples dynamical output simulated with TTVFaster (Agol & Deck 2016), a semianalytic TTV model, with a Mandel & Agol (2002) analytical light-curve model. Optimization of planetary parameters was done by MultiNest (Feroz et al. 2008), a multimodal Bayesian inference algorithm, applied using a unique eccentricity parameterization: the eccentricity differences between adjacent planets as our primary degrees of freedom, which enables the optimization to be significantly less prone to correlations. The median CPU-time that was required for the systems analyzed is ∼6.5 hr per system (single-threaded). The dynamical optimization is thus not CPU-intensive for these few-planet systems. Overall, the reduced number of degrees of freedom in TTVFaster, the combination of the choice of jump parameters and the fitting of the light curve and not the transit timings improved the mass sensitivity here, which allowed the analysis of a sample of Kepler multiplanet systems, and in particular planets that were thus far too small to have well-constrained transit timings. Finally, we quality-check our best-fit values by ensuring that they are at a region of parameter space where TTVFaster is consistent with TTVFast, and that the suggested TTVFaster solution is close to, but not in, resonance.

We were able to provide significant planetary mass (and eccentricity) constraints in several planetary systems: we decreased the uncertainty of Kepler-79 b by a factor of ∼2, and our derived mass values of Kepler-79 d suggests more physically plausible bulk density. We also provide the first mass constraints to Kepler-595 b, Kepler-450 c, and KOI-547.03—Neptune-, sub-Neptune, and Earth-sized planets respectively. The mass of Kepler-450 b is marginally detected to $2.9\sigma$. 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 multitransiting systems from the Kepler and TESS data sets with known TTVs.

We wish to thank Dr. Trifon Trifonov and Yair Judkovsky for helpful discussions, and Prof. Eric Agol for advice on TTVFaster and TTVFast. This study was supported by the Helen Kimmel Center for Planetary Sciences and the Minerva Center for Life Under Extreme Planetary Conditions #13599 at the Weizmann Institute of Science.

In Figures 7–10 we provide a plot for each system considered, in which pairwise correlations are displayed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size