Diffuser-assisted Infrared Transit Photometry for Four Dynamically Interacting Kepler Systems

In this study, we focused on the set of multiplanet systems from the original Kepler survey. We began by estimating the expected TTV signal strength for all planet pairs in order to identify the systems most likely to exhibit strong TTVs. We estimated the minimum mass of a planet from its radius, and then estimated the chopping signal and near-first-order resonant TTV signal for planet pairs given their orbital periods. We then use the number of transits and the transit timing uncertainty to estimate a minimum TTV signal-to-noise ratio (S/N) in the limit of circular orbits. For systems exhibiting TTVs with high S/Ns, we performed dynamical fits to the long-cadence transit times in Rowe & Thompson (2015). We fit five parameters per planet, including the orbital period and phase at a chosen epoch, the two eccentricity vector components, and the dynamical mass. We then mapped the resulting posterior using Differential Evolution Markov Chain Monte Carlo sampling (Jontof-Hutter et al. 2016). Because mutual inclinations are a second-order effect for the TTV amplitude, we assumed coplanarity in our models (Lithwick et al. 2012; Nesvorný & Vokrouhlický 2014; Jontof-Hutter et al. 2016). We then forward-modeled sample solutions for each system in order to identify those with the most strongly diverging TTV predictions. A detailed report of our forward modeling is in preparation.

We selected targets for our WIRC program from the subset of systems with strongly detected TTVs and dynamical solutions that diverged measurably in the years following the end of the primary Kepler mission. We excluded systems where the 1σ range of predicted transit times at the epoch of our proposed WIRC observation was greater than one hour, as this meant that there was a significant possibility that the transit might occur outside our window of observability. In order to ensure that the measured transit time was likely to provide a useful constraint on the dynamical fit, we also calculated the expected timing precision of a new WIRC observation and excluded systems where this uncertainty was greater than the 1σ range in predicted transit times.

Within this sample of systems, we searched for targets with an ingress and/or egress visible from Palomar between 2017 August and 2018 May. We then ranked the targets in our sample based on the predicted S/N scaled from early WIRC commissioning data (Stefansson et al. 2017) and prioritized observations of the highest S/N targets. We ultimately obtained high-quality light curves for four confirmed and candidate planets from this ranked list, including Kepler-29b, Kepler-36c, KOI-1783.01, and Kepler-177c. The predicted midtransit times for these planets are shown in Table 1.

We observed our four selected systems in the J band with WIRC, which is located at the prime focus of the Hale 200 inch telescope at Palomar Observatory (Wilson et al. 2003). The current 2048 × 2048 pixel Hawaii-II HgCdTe detector was installed in 2017 January, along with 32-channel readout electronics that allow for a read time of 0.92 s (Tinyanont et al. 2019). The instrument has an 87 × 87 field of view with a pixel scale of 02487, ensuring that (at least for the magnitude range in our sample) there are always on the order of 10 stars with comparable brightness contained within the same field of view as our target star.

We utilize the custom near-infrared Engineered Diffuser described in Stefansson et al. (2017) to mitigate time-correlated noise from PSF variations and improve our observing efficiency. The diffuser delivers a top-hat PSF with an FWHM of 3''. We also minimize the time-correlated noise contribution from flat-fielding errors by utilizing precision guiding software (Zhao et al. 2014). WIRC does not have a separate guide camera, but instead guides on science images by fitting 2D Gaussian profiles to comparison stars and determining guiding offsets on each image. For these observations, we find that the position of the star typically varies by less than 2–3 pixels over the course of the night, with the largest position drift occurring at high airmass where accurate centroid measurements become more challenging.

Dates, times, and airmasses for each observation are reported in Table 1. For Kepler-29, Kepler-36, and Kepler-177, we observed continuously during the observation windows. During our observation of KOI-1783, there were three breaks in data acquisition, due to a malfunctioning torque motor causing a temporary loss of telescope pointing.

Exposure times are also reported in Table 1 and were chosen to keep the detector in the linear regime. WIRC commissioning tests have shown the detector to be linear to ∼0.5% at 22,000 ADU (Tinyanont et al. 2019). When choosing exposure times, we aimed to keep the maximum count level at or below 20,000 ADU in order to accommodate potential changes in airmass and sky background. In some cases, frames were coadded during the night to increase observing efficiency as noted in Table 1.

For each night, we construct a median dark frame and a flat field. During the construction of the dark and flat, we also construct a global bad pixel map with the procedure described by Tinyanont et al. (2019). Each image is dark-subtracted and flat-fielded, and each bad pixel is replaced with the median of the 5 pixel × 5 pixel box surrounding the errant value. The total number of bad pixels is approximately 0.6% of the full array (Tinyanont et al. 2019). During the calibration sequence, mid-exposure times are converted to Barycentric Julian Date in the Barycentric Dynamical Time standard (BJD_TDB), following the recommendation of Eastman et al. (2010). All of the above steps are performed by the WIRC Data Reduction Pipeline, which was originally developed to automatically handle large sets of polarimetric data (Tinyanont et al. 2019).

We perform aperture photometry using the photutils package (Bradley et al. 2016). We begin by using the first science image as a "finding frame" and detect sources using the DAOStarFinder function (based on Stetson 1987). Sources that are close to the detector edge and those with overlapping apertures are removed automatically. The target star is registered by comparison to an Aladin Lite finding chart (Bonnarel et al. 2000; Boch & Fernique 2014). We then perform the photometry using a range of circular apertures with radii ranging between 6 and 18 pixels in one-pixel steps, using the same aperture for all stars in each image. With WIRC's $\sim 0\buildrel{\prime\prime}\over{.} 25$/pixel scale, the diffuser is expected to deliver stellar PSFs with an FWHM of 12 pixels, but the actual FWHM changes with stellar brightness. For each image, we calculate and subtract the median background via iterative 3σ clipping with the sigma_clipped_stats function in astropy with a five-iteration maximum specified (Astropy Collaboration et al. 2013; Collaboration et al. 2018). After this, we recalculate the source centroids via iterative flux-weighted centroiding and shift apertures accordingly for each individual image. The local sky background is then estimated using an annular region around each source with inner radius of 20 pixels and outer radius of 50 pixels. We find that iterative sigma-clipping of this background region (this time with a 2σ threshold) is sufficient to reconstruct the mean local background, even though the fields are fairly crowded.

After raw light curves are obtained for each aperture size, we choose the 10 comparison stars that best track the time-varying flux of the target star (i.e., those that have minimal variance from the target star). We clean the target and comparison light curves by applying a moving median filter (of width 10 data points) to the target star data set and removing 3σ outliers. We then select the optimal aperture by minimizing the root mean square (rms) scatter after the light-curve fitting described in the next section. Our optimal aperture radii were 8 pixels for Kepler-29b, 14 pixels for Kepler-36c, 10 pixels for KOI-1783.01, and 10 pixels for Kepler-177c. We find that our preferred apertures for each target increase in size with increasing stellar brightness, and all preferred apertures are comparable in size to the aforementioned 12 pixel FWHM expected for the diffuser.

Of the four planets in our sample, only one (Kepler-29b) had a transit duration short enough to allow us to observe a full transit; for the other three planets, our observations spanned ingress or egress, but not both. This introduces a degeneracy between the midtransit time and transit duration (parameterized here by the inclination and semimajor axis) in our fits to these four transits. We resolve this degeneracy by carrying out joint fits with the original Kepler photometry, where we assume common values for the transit depth ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$, the inclination i, and the scaled semimajor axis a/R_⋆. Although we would expect the transit depth to vary as a function of wavelength if any of these planets have atmospheres, the maximum predicted magnitude for this variation (corresponding to a cloud-free, hydrogen-rich atmosphere) is much smaller than our expected measurement uncertainty for the change in transit depth ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$ between the optical Kepler band and our J-band photometry. This effect would be strongest for the low-density planet Kepler-177c, but even then, the maximal variation is of order 200 ppm versus our WIRC J-band precision of roughly 1300 ppm. We found that constraining the transit depth to the Kepler value resulted in smaller transit timing uncertainties for our partial transit observations, which otherwise exhibited correlations between the transit depth, the transit time, and the linear trend in time.

We processed the Kepler long-cadence simple aperture photometry light curves for each star in our sample using the kepcotrend function in the PyKE package (Still & Barclay 2012). To avoid errors in light-curve shape introduced by assuming a linear ephemeris, we cut out individual light curves from the cotrended Kepler data using lists of individual transit times from Holczer et al. (2016) when possible and otherwise using Rowe & Thompson (2015). We selected our trim window to provide two transit durations of both pre-ingress and post-egress baseline. After dividing out a linear trend fit to the out-of-transit baseline for each light curve, we combined all transits into a single transit light curve with flux as a function of time from transit center.

This process assumes that TDVs do not strongly bias our retrieved transit shapes. For systems with large-amplitude TDVs, it may become necessary to perform photodynamical modeling in order to properly treat the time-varying transit shape (e.g., Freudenthal et al. 2018). However, Holczer et al. (2016) examined data spanning the full length of the Kepler mission and did not detect TDVs for any of the targets in our sample. To further justify our assumption that TDVs have a negligible impact on the measured signals, we calculated the expected TDV amplitude for Kepler-177c (a planet with long period and large impact parameter that is more prone to nodal precession). The maximum TDV amplitude is of order 0.1 hr over the 10 yr baseline. The WIRC data alone are not sensitive to transit duration changes on this timescale, because we only detect ingress or egress for most transits. Additionally, the precision on the transit timing in the joint fits tend to be much more uncertain than 0.1 hr, meaning that TDV effects will not compromise our final TTV constraints. We conclude that we can safely ignore TDVs in our treatment of these data.

To fit the Kepler and WIRC light curves, we first constructed light-curve models defined by observed quantities and fit parameters. We then constructed appropriate likelihood and prior functions and sampled the resultant posterior probability numerically to obtain estimates of the best-fit parameters and their associated uncertainties. The outputs of the WIRC photometry pipeline are an array of times ${\boldsymbol{t}}=({t}_{1},{t}_{2},\,\ldots ,\,{t}_{n})$, the target data array ${\boldsymbol{y}}=({y}_{1},{y}_{2},\,\ldots ,\,{y}_{n})$ (with y_i referring to the measurement at time t_i), and comparison star arrays ${{\boldsymbol{x}}}_{j}\,=({x}_{1},{x}_{2},\,\ldots ,\,{x}_{n})$. Collectively, the comparison stars define a matrix X, with one comparison star ${{\boldsymbol{x}}}_{j}$ in each row of the matrix.

We aim to fit the target ${\boldsymbol{y}}$ with a model ${\boldsymbol{M}}$ that depends on the depth of the transit ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$, the transit center time t₀, the inclination i, the ratio of the semimajor axis to stellar radius a/R_⋆, and a linear trend in time α. That model can be written as follows (loosely following the notation of Diamond-Lowe et al. 2018):

$$\begin{eqnarray}&&{\boldsymbol{M}}=[\alpha {\boldsymbol{t}}+{\boldsymbol{S}}]\times {{\boldsymbol{T}}}_{\mathrm{WIRC}}({({R}_{{\rm{p}}}/{R}_{\star })}^{2},{t}_{0},i,a/{R}_{\star }),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{S}}$ is the systematics model, ${{\boldsymbol{T}}}_{\mathrm{WIRC}}$ is the transit model, and the multiplication is meant to denote a pointwise product. We use the batman code to construct the transit model (Kreidberg 2015) and fix the planet eccentricities to zero. The eccentricities of multiplanet Kepler systems are typically small, with a population mean of $\bar{e}={0.04}_{-0.04}^{+0.03}$ (Xie et al. 2016), and the effect of these eccentricities on the shape of the transit light curve is negligible for these data. We use four-parameter nonlinear limb-darkening coefficients from Claret & Bloemen (2011), assuming stellar parameter values from Petigura et al. (2017a) that are reproduced in Table 2.

Table 2.  Stellar Parameters for the Stars in Our Sample

Target	Teff	[Fe/H]	log g	${M}_{\star }$	${R}_{\star }$
	(K)	(dex)	(log(cm s−2))	(${M}_{\odot }$)	(${R}_{\odot }$)
Kepler-29	${5378}_{-60}^{+60}$	$-{0.44}_{-0.04}^{+0.04}$	${4.6}_{-0.1}^{+0.1}$	${0.761}_{-0.028}^{+0.024}$	${0.732}_{-0.031}^{+0.033}$
Kepler-36	${5979}_{-60}^{+60}$	$-{0.18}_{-0.04}^{+0.04}$	${4.1}_{-0.1}^{+0.1}$	${1.034}_{-0.022}^{+0.022}$	${1.634}_{-0.040}^{+0.042}$
KOI-1783	${5922}_{-60}^{+60}$	${0.11}_{-0.04}^{+0.04}$	${4.3}_{-0.1}^{+0.1}$	${1.076}_{-0.032}^{+0.036}$	${1.143}_{-0.030}^{+0.031}$
Kepler-177	${5732}_{-60}^{+60}$	$-{0.11}_{-0.04}^{+0.04}$	${4.1}_{-0.1}^{+0.1}$	${0.921}_{-0.023}^{+0.025}$	${1.324}_{-0.051}^{+0.053}$

Note. Spectroscopic parameters (T_eff, [Fe/H], and log g) are taken from Fulton et al. (2017) and physical parameters (${M}_{\star }$ and ${{\rm{R}}}_{\star }$) are from Fulton & Petigura (2018).

Download table as:  ASCIITypeset image

For ground-based observations, we expect the measured flux from each star to vary as a function of the airmass, centroid drift, seeing changes, transparency variations, and other relevant parameters. However, all of the stars on our wide-field detector should respond similarly to changes in the observing conditions. In particular, we expect that stars of approximately the same J magnitude and color will track closely with the light curve of our target star. We therefore define our systematics model as a linear combination of comparison star light curves. This allows us to empirically model these effects without explicitly relating them to the relevant atmospheric and telescope state parameters via a parametric model. We determine the coefficients for the linear combination via a linear regression fit to the target light curve after dividing out the transit light-curve model (which we call the "target systematics" ${{\boldsymbol{S}}}_{\mathrm{target}}$). We calculate new linear coefficients every time the transit light curve is modified. Mathematically, the target systematics can be written as

$$\begin{eqnarray}&&{{\boldsymbol{S}}}_{\mathrm{target}}=\displaystyle \frac{{\boldsymbol{y}}}{{{\boldsymbol{T}}}_{\mathrm{WIRC}}({({R}_{{\rm{p}}}/{R}_{\star })}^{2},{t}_{0},i,a/{R}_{\star })}-\alpha {\boldsymbol{t}},\end{eqnarray} \tag{ 2 }$$

where division is meant to be pointwise, and the linear regression defining the systematics model can be written as

$$\begin{eqnarray}&&{\boldsymbol{S}}={{\boldsymbol{PS}}}_{\mathrm{target}},\end{eqnarray} \tag{ 3 }$$

where the projection matrix P comes from the comparison stars and can be written as

$$\begin{eqnarray}&&{\boldsymbol{P}}={{\boldsymbol{X}}}^{T}{\left({{\boldsymbol{XX}}}^{T}\right)}^{-1}{\boldsymbol{X}}\end{eqnarray} \tag{ 4 }$$

Equations (1)–(4) thus define the model ${\boldsymbol{M}}$ solely as a function of the observed quantities {${\boldsymbol{t}},{\boldsymbol{y}},{\boldsymbol{X}}$} and the fit parameters {${({R}_{{\rm{p}}}/{R}_{\star })}^{2},{t}_{0},\alpha ,i,a/{R}_{\star }$}. To give a sense for how our systematics removal looks in practice, in Figure 2 we show the raw and detrended light curves for KOI-1783.01 along with the best systematics and transit models.

Zoom In Zoom Out Reset image size

As discussed in Section 3.2, we fit the WIRC photometry jointly with the Kepler photometry in order to avoid a strong degeneracy between midtransit time and transit duration. The Kepler photometry consists of an array of times ${{\boldsymbol{t}}}_{{Kep}}\,=({t}_{1},{t}_{2},\,\ldots ,\,{t}_{n})$ and the corresponding detrended target data array ${{\boldsymbol{y}}}_{{Kep}}=({y}_{1},{y}_{2},\,\ldots ,\,{y}_{n})$. Because these data are already detrended and phased together, the model ${{\boldsymbol{M}}}_{{Kep}}$ for the Kepler data is simply a batman transit model:

$$\begin{eqnarray}&&{{\boldsymbol{M}}}_{{Kep}}={{\boldsymbol{T}}}_{{Kep}}({({R}_{{\rm{p}}}/{R}_{\star })}^{2},i,a/{R}_{\star })\end{eqnarray} \tag{ 5 }$$

We supersampled the Kepler light curves to 1 minute cadence and used four-parameter nonlinear limb-darkening coefficients from Sing (2010) calculated specifically for the Kepler bandpass.

Having defined our models, we can now define our likelihood function. We assume measurements to be Gaussian-distributed and uncorrelated (correlated noise is considered briefly in Section 4.1) such that the likelihood takes the form

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({ \mathcal L }) & = & -\displaystyle \frac{1}{2}\sum _{i}\mathrm{log}(2\pi {\sigma }_{i}^{2})-\displaystyle \frac{1}{2}\sum _{i}{\left(\displaystyle \frac{{y}_{i}-{M}_{i}}{{\sigma }_{i}}\right)}^{2}\\ & & -\displaystyle \frac{1}{2}\sum _{i}\mathrm{log}(2\pi {\sigma }_{{Kep},i}^{2})\\ & & -\displaystyle \frac{1}{2}\sum _{i}{\left(\displaystyle \frac{{y}_{{Kep},i}-{M}_{{Kep},i}}{{\sigma }_{{Kep},i}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 6 }$$

where the uncertainties σ_i and ${\sigma }_{{Kep},i}$ are quadrature sums of the Poisson noise from the target star and extra noise terms that can be fitted:

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}=\sqrt{{{\boldsymbol{\sigma }}}_{\mathrm{phot},\mathrm{WIRC}}^{2}+{\sigma }_{\mathrm{extra},\mathrm{WIRC}}^{2}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\sigma }}}_{{Kep}}=\sqrt{{{\boldsymbol{\sigma }}}_{\mathrm{phot},{Kep}}^{2}+{\sigma }_{\mathrm{extra},{Kep}}^{2}}.\end{eqnarray} \tag{ 8 }$$

Because the extra noise terms are always positive, we fit for $\mathrm{log}({\sigma }_{\mathrm{extra},\mathrm{WIRC}})$ and $\mathrm{log}({\sigma }_{\mathrm{extra},{Kep}})$ as a numerical convenience. Also, rather than fitting for t₀ itself, we define all times relative to the predicted transit times in Table 1 and fit for the offset from that time Δt₀.

We impose priors on all parameters. They are either Gaussian, taking the functional form

$$\begin{eqnarray}&&\mathrm{log}({{ \mathcal P }}_{k})=-\displaystyle \frac{1}{2}\mathrm{log}(2\pi {\sigma }_{k}^{2})-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{k-{\mu }_{k}}{{\sigma }_{k}}\right)}^{2},\end{eqnarray} \tag{ 9 }$$

or uniform, taking the functional form

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({{ \mathcal P }}_{k}) & = & \mathrm{log}\left(\displaystyle \frac{1}{{k}_{\max }-{k}_{\min }}\right),{k}_{\min }\lt k\lt {k}_{\max };\\ & & -\infty \,\mathrm{otherwise}.\end{array}\end{eqnarray} \tag{ 10 }$$

We placed physically motivated Gaussian priors on a/R_⋆ calculated from the stellar parameters reported by Fulton & Petigura (2018) and used uniform priors for all other variables. We list our priors for the physical fit parameters in Table 4.

Table 4.  System Parameters for the Joint Photometric Fits

Parameter	Symbol	Values				Units	Source
		Kepler-29b	Kepler-36c	KOI-1783.01	Kepler-177c
Fixed Parameters
Orbital period	P	10.3392924	16.23192004	134.4786723	49.41117582	days	(1, 2)
Predicted transit time	${t}_{0}$	2457990.852	2458023.9115	2458229.7971125	2458242.93807	BJD	⋯
Eccentricity	e	0.	0.	0.	0.	⋯	⋯
Kepler limb-darkening coefficients	${a}_{1}$	0.4959	0.4639	0.6034	0.5716	⋯	(3)
	${a}_{2}$	0.0222	0.3045	−0.1382	−0.1145	⋯	(3)
	${a}_{3}$	0.5708	0.0751	0.6330	0.6579	⋯	(3)
	${a}_{4}$	−0.3485	−0.1251	−0.3506	−0.3667	⋯	(3)
WIRC limb-darkening coefficients	${b}_{1}$	0.3634	0.3982	0.4832	0.4421	⋯	(4)
	${b}_{2}$	0.5846	0.5452	0.2998	0.3993	⋯	(4)
	${b}_{3}$	−0.6152	−0.6817	−0.3634	−0.4523	⋯	(4)
	${b}_{4}$	0.1997	0.2508	0.1152	0.1474	⋯	(4)
Fit Priors
Transit depth prior	${{ \mathcal P }}_{{({R}_{{\rm{p}}}/{R}_{\star })}^{2}}$	${ \mathcal U }(0,2000)$	${ \mathcal U }(0,1000)$	${ \mathcal U }(0,10000)$	${ \mathcal U }(0,8000)$	ppm	⋯
Transit timing offset prior	${{ \mathcal P }}_{{\rm{\Delta }}{t}_{0}}$	${ \mathcal U }(-100,100)$	${ \mathcal U }(-100,100)$	${ \mathcal U }(-100,100)$	${ \mathcal U }(-100,100)$	min	⋯
Inclination prior	${{ \mathcal P }}_{i}$	${ \mathcal U }(85,90)$	${ \mathcal U }(85,90)$	${ \mathcal U }(85,90)$	${ \mathcal U }(85,90)$	°	⋯
Scaled semimajor axis prior	${{ \mathcal P }}_{a/{R}_{\star }}$	${ \mathcal N }(24.906,1.125)$	${ \mathcal N }(16.696,0.436)$	${ \mathcal N }(99.030,2.840)$	${ \mathcal N }(41.649,1.674)$	⋯	(5)
Fit Posteriors
Transit depth	${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$	${1020}_{-34}^{+31}$	${425.3}_{-3.5}^{+3.8}$	${5044}_{-64}^{+87}$	${3643}_{-57}^{+55}$	ppm	⋯
Transit timing offset	${\rm{\Delta }}{t}_{0}$	−14.3${}_{-2.7}^{+16.7}$	−17.9${}_{-4.7}^{+11.8}$	${16}_{-11}^{+10}$	${45.2}_{-7.1}^{+8.7}$	min	⋯
Inclination	i	${89.13}_{-0.23}^{+0.45}$	${89.36}_{-0.29}^{+0.45}$	${89.4413}_{-0.0082}^{+0.0076}$	${88.795}_{-0.035}^{+0.037}$	°	⋯
Scaled semimajor axis	$a/{R}_{\star }$	${24.95}_{-0.91}^{+1.34}$	${16.69}_{-0.31}^{+0.26}$	${94.8}_{-1.1}^{+1.1}$	${42.08}_{-0.94}^{+1.04}$	⋯	⋯
Derived Parameters
Planet–star radius ratio	${R}_{{\rm{p}}}/{R}_{\star }$	${0.03194}_{-0.00054}^{+0.00048}$	${0.02062}_{-0.00009}^{+0.00009}$	${0.07102}_{-0.00045}^{+0.00061}$	${0.06036}_{-0.00047}^{+0.00045}$	⋯	⋯
Impact Parameter	b	${0.379}_{-0.185}^{+0.083}$	${0.186}_{-0.131}^{+0.080}$	${0.9239}_{-0.0023}^{+0.0026}$	${0.8848}_{-0.0065}^{+0.0056}$	⋯	⋯
Transit duration	${T}_{14}$	${3.041}_{-0.052}^{+0.045}$	${7.46}_{-0.017}^{+0.021}$	${5.874}_{-0.040}^{+0.039}$	${5.243}_{-0.054}^{+0.054}$	hr	⋯

Note. (1) Morton et al. (2016), (2) Thompson et al. (2018), (3) Sing (2010), (4) Claret & Bloemen (2011), and (5) Fulton & Petigura (2018). Also, ${ \mathcal N }(a,b)$ indicates a normal (Gaussian) prior with mean a and standard deviation b described by Equation (9), whereas ${ \mathcal U }(a,b)$ indicates a uniform prior with lower bound a and upper bound b described by Equation (10).

Download table as:  ASCIITypeset image

With the likelihood and priors defined, we can finally write the posterior probability with Bayes' Theorem (up to a constant proportional to the evidence):

$$\begin{eqnarray}&&\mathrm{log}(\mathrm{Prob})=\mathrm{log}({ \mathcal L })+\sum _{k}\mathrm{log}({{ \mathcal P }}_{k})\end{eqnarray} \tag{ 11 }$$

Then, we seek a solution for the fit parameters ${({R}_{{\rm{p}}}/{R}_{\star })}^{2},{\rm{\Delta }}{t}_{0},i,a/{R}_{\star },\alpha ,\mathrm{log}({\sigma }_{\mathrm{extra},\mathrm{WIRC}})$, and $\mathrm{log}({\sigma }_{\mathrm{extra},{Kep}})$ that maximizes $\mathrm{log}(\mathrm{Prob})$. We carry out an initial fit using scipy's Powell minimizer (Jones et al. 2001) and use this solution as a starting point for the affine-invariant ensemble Markov Chain Monte Carlo sampler emcee (Foreman-Mackey et al. 2013). We burn the chains in for 2 × 10³ steps and then run for 10⁵ steps. This corresponds to at least 500 integrated autocorrelation times for each parameter. The maximum a posteriori parameter estimates with associated 68% confidence intervals for all model parameters aside from $\alpha ,\mathrm{log}({\sigma }_{\mathrm{extra},\mathrm{WIRC}})$, and $\mathrm{log}({\sigma }_{\mathrm{extra},{Kep}})$ are given in Table 4. The best-fit light curves are shown in Appendix A. Additionally, we plot the posterior distributions for these parameters in Appendix B.

Our fits to the ground-based WIRC photometry typically resulted in a non-Gaussian posterior for the midtransit time. We accounted for these skewed distributions in our dynamical fits by dividing the posteriors into 20 bins and normalized the probability density to give a likelihood for each bin, as illustrated in the marginalized timing distributions from Appendix B. We then ran two sets of dynamical fits for each system using either these skewed timing posteriors or a symmetric Gaussian distribution with a width equal to the average of our positive and negative uncertainties.

We fitted dynamical models to the transit timing data using a Differential Evolution Markov Chain Monte Carlo algorithm (Ter Braak 2006; Nelson et al. 2014; Jontof-Hutter et al. 2015, 2016). We used uniform priors for the orbital period and phase, and uniform positive definite priors for the dynamical masses. For each eccentricity vector component, we assumed a Gaussian distribution centered on 0 with a width of 0.1 for the prior. This is wider than the inferred eccentricity distribution among Kepler's multiplanet systems (Fabrycky et al. 2014; Hadden & Lithwick 2014), but TTV modeling is subject to an eccentricity–eccentricity degeneracy whereby aligned orbits can have larger eccentricities than allowed by our prior with little effect on the relative eccentricity (Jontof-Hutter et al. 2016). The results of our dynamical modeling are given in Table 5. This table includes orbital periods (solved at our chosen epoch of BJD = 2,455,680), masses, and eccentricity vectors for retrievals with only the Kepler data, retrievals including the new WIRC transit time with a Gaussian uncertainty distribution, and retrievals using the skewed WIRC timing posterior. We find that our fits using Gaussian posteriors are generally in good agreement with results from fits utilizing the skewed transit timing posteriors.

Table 5.  Results from Our Dynamical Analysis

Planet	Data Set	P [days]	$\left(\tfrac{{M}_{{\rm{p}}}}{{M}_{\oplus }}\right)\left(\tfrac{{M}_{\odot }}{{M}_{\star }}\right)$	$e\cos (\omega )$	$e\sin (\omega )$
Kepler-29b	Kep LC	${10.33838}_{-0.00027}^{+0.00030}$	4.6${}_{-1.5}^{+1.4}$	−0.060${}_{-0.071}^{+0.072}$	−0.030${}_{-0.072}^{+0.072}$
	Kep LC + WIRC (G)	${10.33974}_{-0.00015}^{+0.00014}$	3.7${}_{-1.3}^{+1.3}$	0.013${}_{-0.071}^{+0.071}$	−0.016${}_{-0.063}^{+0.056}$
	Kep LC + WIRC (S)	${10.33966}_{-0.00017}^{+0.00015}$	3.8${}_{-1.0}^{+1.1}$	0.003${}_{-0.070}^{+0.068}$	−0.088${}_{-0.058}^{+0.059}$
Kepler-29c	Kep LC	${13.28843}_{-0.00053}^{+0.00048}$	4.07${}_{-2.29}^{+2.87}$	0.007${}_{-0.062}^{+0.063}$	−0.022${}_{-0.063}^{+0.063}$
	Kep LC + WIRC (G)	${13.28613}_{-0.00021}^{+0.00026}$	3.28${}_{-1.08}^{+1.06}$	−0.023${}_{-0.062}^{+0.061}$	−0.022${}_{-0.055}^{+0.045}$
	Kep LC + WIRC (S)	${13.28633}_{-0.00027}^{+0.00031}$	3.39${}_{-0.84}^{+0.86}$	−0.007${}_{-0.061}^{+0.059}$	−0.085${}_{-0.051}^{+0.051}$
Kepler-36b	Kep LC	${13.86834}_{-0.00051}^{+0.00050}$	3.990${}_{-0.092}^{+0.093}$	0.050${}_{-0.025}^{+0.023}$	−0.026${}_{-0.033}^{+0.034}$
	Kep LC + WIRC (G)	${13.86825}_{-0.00050}^{+0.00050}$	3.972${}_{-0.074}^{+0.078}$	0.041${}_{-0.020}^{+0.019}$	−0.011${}_{-0.018}^{+0.018}$
	Kep LC + WIRC (S)	${13.86821}_{-0.00049}^{+0.00049}$	3.964${}_{-0.068}^{+0.077}$	0.037${}_{-0.018}^{+0.019}$	−0.004${}_{-0.015}^{+0.012}$
Kepler-36c	Kep LC	${16.21867}_{-0.00010}^{+0.00010}$	7.456${}_{-0.168}^{+0.167}$	0.053${}_{-0.023}^{+0.021}$	−0.039${}_{-0.031}^{+0.031}$
	Kep LC + WIRC (G)	${16.21865}_{-0.00010}^{+0.00010}$	7.397${}_{-0.107}^{+0.104}$	0.046${}_{-0.018}^{+0.017}$	−0.026${}_{-0.017}^{+0.017}$
	Kep LC + WIRC (S)	${16.21865}_{-0.00010}^{+0.00010}$	7.371${}_{-0.093}^{+0.092}$	0.042${}_{-0.016}^{+0.017}$	−0.019${}_{-0.014}^{+0.012}$
KOI-1783.01	Kep LC	${134.4622}_{-0.0038}^{+0.0035}$	${90.2}_{-23.2}^{+30.3}$	0.0079${}_{-0.0050}^{+0.0080}$	−0.039${}_{-0.021}^{+0.012}$
	Kep LC + WIRC (G)	${134.4628}_{-0.0035}^{+0.0033}$	${78.1}_{-12.9}^{+15.1}$	0.0073${}_{-0.0046}^{+0.0067}$	−0.048${}_{-0.015}^{+0.014}$
	Kep LC + WIRC (S)	${134.4629}_{-0.0036}^{+0.0033}$	${76.4}_{-9.6}^{+11.8}$	0.0072${}_{-0.0045}^{+0.0067}$	−0.049${}_{-0.012}^{+0.014}$
KOI-1783.02	Kep LC	${284.230}_{-0.031}^{+0.044}$	${17.1}_{-4.3}^{+5.1}$	0.018${}_{-0.015}^{+0.018}$	−0.011${}_{-0.032}^{+0.027}$
	Kep LC + WIRC (G)	${284.215}_{-0.021}^{+0.026}$	${16.2}_{-3.8}^{+4.7}$	0.017${}_{-0.015}^{+0.015}$	−0.020${}_{-0.028}^{+0.034}$
	Kep LC + WIRC (S)	${284.212}_{-0.018}^{+0.024}$	${16.1}_{-3.8}^{+4.6}$	0.017${}_{-0.014}^{+0.015}$	−0.020${}_{-0.026}^{+0.034}$
Kepler-177b	Kep LC	${35.8591}_{-0.0017}^{+0.0019}$	5.76${}_{-0.81}^{+0.84}$	−0.026${}_{-0.075}^{+0.074}$	−0.014${}_{-0.068}^{+0.065}$
	Kep LC + WIRC (G)	${35.8601}_{-0.0014}^{+0.0015}$	5.44${}_{-0.75}^{+0.78}$	0.017${}_{-0.054}^{+0.052}$	−0.001${}_{-0.063}^{+0.062}$
	Kep LC + WIRC (S)	${35.8601}_{-0.0012}^{+0.0013}$	5.38${}_{-0.74}^{+0.78}$	0.020${}_{-0.048}^{+0.047}$	0.005${}_{-0.061}^{+0.061}$
Kepler-177c	Kep LC	${49.40964}_{-0.00097}^{+0.00097}$	${14.6}_{-2.5}^{+2.7}$	−0.027${}_{-0.065}^{+0.064}$	−0.014${}_{-0.059}^{+0.056}$
	Kep LC + WIRC (G)	${49.40926}_{-0.00077}^{+0.00078}$	${13.9}_{-2.5}^{+2.7}$	0.010${}_{-0.046}^{+0.045}$	−0.003${}_{-0.054}^{+0.053}$
	Kep LC + WIRC (S)	${49.40921}_{-0.00074}^{+0.00072}$	${13.5}_{-2.3}^{+2.5}$	0.013${}_{-0.041}^{+0.040}$	0.003${}_{-0.053}^{+0.052}$

Note. In the Data Set column, "Kep LC" refers to the transit timings from the Kepler long-cadence light curves, "WIRC (G)" refers to the transit timing from our observations when assumed to have Gaussian uncertainties, and "WIRC (S)" refers to the transit timing from our observations taking into account the skewed shape of our timing posteriors. Also, the orbital period P is solved for at our chosen epoch of BJD = 2,455,680.

Download table as:  ASCIITypeset image

Our best photometric performance is for Kepler-177c, where we were only ∼20% above the shot noise. We also investigate how well WIRC mitigates time-correlated noise, which can lead to underestimated uncertainties in reported transit times. We calculate the rms versus bin size for each observation and show the corresponding plots in the bottom-right panels of Figures 7–10. We find that Kepler-29b and KOI-1783.01 appear to have minimal time-correlated noise (see the bottom-right panels in Figures 7 and 9, respectively). Kepler-36c has some time-correlated trends on longer timescales, and for Kepler-177c, quasi-periodic noise is readily visible in both the best-fit residual plot and in the rms versus bin size plot (see also the bottom-right panel in Figures 8 and 10, respectively). We tried adding sinusoids to our fits for these planets, but found that this had a negligible effect on the overall quality of the fits and the resulting transit timing posteriors.

To derive a representative noise statistic for WIRC, we first calculated the scatter in 10 minute bins for each of our observations. These statistics were then scaled to the equivalent values for observations of a 14th magnitude star. In some of our earliest observations, we used a suboptimal coaddition strategy, resulting in relatively inefficient observations (for Kepler-36c, this increased the noise by 31.1% relative to a more optimal strategy). We therefore applied an additional correction factor to rescale the noise for these inefficient observations to the expected value for better-optimized observations. Averaging these corrected noise statistics together, we find that WIRC can deliver 1613 ppm photometry per 10 minute bin on a J = 14 magnitude star. If we assume that we are able to collect two hours of data in transit and two hours out of transit, this equates to a precision of 659 ppm on the transit depth measurement for planets around a J = 14 magnitude star. To highlight the range of parameter space that this precision opens up, we plot transit depths for all confirmed transiting exoplanets against the host star J magnitude in Figure 3 along with the 3σ detection thresholds of WIRC and Spitzer. While Spitzer performs better for brighter stars, WIRC begins to outperform Spitzer for stars fainter than ∼10 magnitude, doing a factor of 1.6 better at J = 14. In practice, the achieved photometric precision will also depend on factors such as atmospheric background, amount of baseline obtained, diurnal constraints, and the number of available comparison stars of comparable magnitude, but the first-order considerations in Figure 3 suggest that ground-based, diffuser-assisted infrared photometry can indeed outperform some current space-based facilities for typical Kepler transiting planet systems.

Zoom In Zoom Out Reset image size

Kepler-29b is a sub-Neptune near the 5:4 and 9:7 mean-motion resonances with the sub-Neptune Kepler-29c. Both low-density planets were originally confirmed by Fabrycky et al. (2012) using TTVs; subsequent dynamical analyses have shown that the pair may actually be in the second-order 9:7 resonance (Migaszewski et al. 2017), but the TTV curve is likely also affected by proximity to the first-order 5:4 resonance (Jontof-Hutter et al. 2016). We detect a transit of Kepler-29b at 3.5σ confidence in the WIRC data. The final detrended Kepler and WIRC light curves, models, residuals, and rms binning plots for Kepler-29b are shown in Figure 7, and the corresponding posterior probability distributions are shown in Figure 11. Although the transit shape is poorly constrained by the WIRC data alone, both ingress and egress are visible by eye in the WIRC light curve, and the relative timing of these two events provides a solid estimate of the transit time when we constrain the transit shape using the Kepler photometry. We find that the resulting posterior distribution for our new WIRC transit time is fairly asymmetric, with the final timing offset determined to $-{14}_{-3}^{+17}$ minutes.

Our new observation was obtained in an epoch where the Kepler-only dynamical fits yield substantially divergent transit times, and as a result, our new transit time provides an improved constraint on the planet masses and eccentricities as shown in Figure 15. We find that the dynamical mass estimate for Kepler-29c has improved by almost a factor of 3 in our updated fits. Our new results favor dynamical masses on the low side of (but not incompatible with) the mass distributions inferred by Jontof-Hutter et al. (2016) for Kepler-29b and c.

Despite these decreased masses, our updated densities for these planets (1.7 ± 0.5 and 1.9 ± 0.5 g cm⁻³, respectively) are larger than the densities reported by Jontof-Hutter et al. (2016). This is because we utilize updated stellar parameters of $M\,={0.761}_{-0.028}^{+0.024}\,{M}_{\odot }$ and $R={0.732}_{-0.031}^{+0.033}\,{R}_{\odot }$ from Fulton & Petigura (2018), which are smaller than the values of $M=0.979\,\pm 0.052\,{M}_{\odot }$ and $R=0.932\pm 0.060\,{M}_{\odot }$ adopted by Jontof-Hutter et al. (2016). For a fixed planet–star radius ratio, a smaller stellar radius implies a correspondingly smaller planet radius. Similarly, a smaller stellar mass implies a larger planet mass for the same best-fit dynamical mass ratio. Both changes therefore act to increase the measured planetary density. Even with these increased density estimates, it is likely that both of these planets have retained a modest hydrogen-rich atmosphere (see Section 5.2.1). The masses and radii of both planets also remain quite similar, in good agreement with the "peas in a pod" trend wherein multiplanet Kepler systems tend to host planets that are similar in both size and bulk density (Millholland et al. 2017; Weiss et al. 2018).

The Kepler-36 system includes two planets with strikingly dissimilar densities: Kepler-36b is a rocky super-Earth close to 7:6 mean-motion resonance with the low-density sub-Neptune Kepler-36c (Carter et al. 2012). The latter planet was included in our sample, and we detect it with a significance of 5.3σ. We present the final light curves and associated statistics for our new transit observation of Kepler-36c in Figure 8 and plot the corresponding posteriors in Figure 12. The posterior distribution on the WIRC transit time is again fairly asymmetric, with the offset constrained to $-{18}_{-5}^{+12}$ minutes. We obtain masses and densities for both planets consistent with previous investigations (though on the low side for Kepler-36b; Carter et al. 2012; Hadden & Lithwick 2017). In Figure 16, we provide updated dynamical masses, eccentricity vectors, and transit timing for this system. Future constraints from TESS should allow for improved mass estimates in this system, especially for Kepler-36c (Goldberg et al. 2019).

The rms scatter achieved for this measurement was 2× the photon noise limit (see bottom-right panel of Figure 8), which is higher than any of the other observations presented in this work. This is due in part to scintillation noise (Stefansson et al. 2017), as Kepler-36 was our brightest target, and we used correspondingly short integration times. For this star, the scintillation noise at an airmass of 1.5 is ∼650 ppm, which is comparable to the shot noise. Our use of short integration times also limited our observing efficiency, resulting in higher photometric scatter than might otherwise have been expected for this relatively bright star. Both problems could be mitigated by increasing the number of coadds, resulting in a longer effective integration time and higher overall observing efficiency.

As we will discuss in Section 5.1, there is already compelling evidence in the literature establishing the planetary nature of this system, which contains two long-period (134 and 284 days, respectively) gas-giant planet candidates located near a 2:1 period commensurability. We present the final light curves and associated statistics for our new transit observation of KOI-1783.01 in Figure 9 and plot the corresponding posteriors in Figure 13. This planet is detected with a significance of 5.9σ in the WIRC data, and we achieve a timing precision of about 10 minutes. These results are in good agreement with a model of the KOI-1783 system that assumes the source of TTVs to be near-resonant planet–planet perturbations. In Figure 17, we present updated constraints on dynamical masses, eccentricities, and transit timing for KOI-1783. Our new transit observation reduces the uncertainty on the dynamical mass of KOI-1783.01 by approximately a factor of 2. When combined with the stellar parameters from Fulton & Petigura (2018), these new constraints provide the most detailed picture of this system to date. We find that KOI-1783.01 is slightly smaller than Saturn, with ${R}_{{\rm{p}}}={8.9}_{-0.2}^{+0.3}{R}_{\oplus }$ and ${M}_{{\rm{p}}}={71}_{-9}^{+11}{M}_{\oplus }$. This corresponds to a density of $\rho \,={0.56}_{-0.09}^{+0.10}$ g cm⁻³, consistent with the presence of a substantial gaseous envelope; we discuss the corresponding implications for this planet's bulk composition in more detail in Section 5.2.2. KOI-1783.02 has a mass of ${M}_{{\rm{p}}}={15}_{-4}^{+4}{M}_{\oplus }$, a radius of ${R}_{{\rm{p}}}={5.4}_{-0.3}^{+0.5}{R}_{\oplus }$, and a density of $\rho ={0.5}_{-0.2}^{+0.2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, again indicative of a substantial gaseous envelope. Both planets appear to have low orbital eccentricities ($e\lesssim 0.05$), in agreement with the overall Kepler TTV sample (Fabrycky et al. 2014; Hadden & Lithwick 2014; Xie et al. 2016). Additionally, we note that the uncertainty on $e\cos (\omega )$ for KOI-1783.01 is an order of magnitude lower than for the other planets in this study, corresponding to a $\pm 1\sigma$ uncertainty of approximately 13 hr in the secondary eclipse phase. Although this is quite good for a planet on a 134 day orbit, the star's faintness and the planet's low equilibrium temperature make this a challenging target for secondary eclipse observations.

The Kepler-177 system contains a low-density sub-Neptune (Kepler-177b) and a very low-density sub-Neptune (Kepler-177c) located near the 4:3 mean-motion resonance. This system was initially confirmed via TTVs by Xie (2014) and subsequently reanalyzed by Jontof-Hutter et al. (2016) and Hadden & Lithwick (2017). Our final light curves and associated statistics for Kepler-177c are given in Figure 10, and the posteriors are given in Figure 14. We detect the transit at 5.5σ significance and measure the corresponding transit time with a 1σ uncertainty of approximately 10 minutes. Although our new dynamical fits for this system (shown in Figure 18) result in modestly lower mass uncertainties, our transit observation was taken close to one TTV super-period away from the Kepler data, where diverging solutions reconverge and thus our new observations provided limited leverage to constrain these dynamical models. If the TESS mission is extended, it should provide additional transit observations that would further reduce the mass uncertainties in this system (Goldberg et al. 2019), but our observations demonstrate that this system is also accessible to ground-based follow-up at a more favorable epoch.

As the only unverified planet candidate in our sample, KOI-1783.01 represents a special case for this program. A transiting planet candidate around KOI-1783 (KIC 10005758) was first reported by Batalha et al. (2013), and a second candidate in the system was identified by the Planet Hunters citizen science collaboration (Lintott et al. 2013). While the a priori probability of both transit signals being false positives is quite low (Lissauer et al. 2011, 2012, 2014; Lintott et al. 2013; Rowe et al. 2014), a few characteristics of this system precluded a quick confirmation. First, the transit signals for both candidates are near-grazing (the grazing parameter $X=b+{R}_{{\rm{p}}}/{R}_{\star }$ is ${0.9949}_{-0.0027}^{+0.0032}$ for KOI-1783.01 from our posteriors and ${0.932}_{-0.015}^{+0.065}$ for KOI-1783.02 from the Thompson et al. 2018 catalog), with "V"-shaped morphologies that Batalha et al. (2013) noted as being potentially diagnostic of an eclipsing binary. Additionally, the Kepler Data Validation reports show a fairly large offset ($\sim 0\buildrel{\prime\prime}\over{.} 25$) of the stellar centroid during the transit relative to the KIC position, which is also typical of stellar blends.

The two transit candidates in this system have a period ratio of 2.11, near the 2:1 commensurability. Such an architecture can generate detectable TTVs, which previous studies have used to confirm the planetary nature of transit candidates (Nesvorný et al. 2013; Steffen et al. 2013). Early analyses of the transit times of KOI-1783.01 (Ford et al. 2012; Mazeh et al. 2013) noted the potential presence of TTVs, but concluded that the significance of the deviation from a linear ephemeris was too low to be conclusive. As Kepler continued to observe this target, evidence for TTVs of both planet candidates in this system grew stronger (Rowe et al. 2014; Holczer et al. 2016). An independent analysis of this system by the Hunt for Exomoons with Kepler Project found evidence for dynamical interactions (Kipping et al. 2015), selecting a TTV model over a linear ephemeris model by 17.2σ for KOI-1783.02. The spectral TTV analysis of Ofir et al. (2018) also found evidence of dynamical interactions, yielding Δχ² values for the TTV signals over a linear model of 49 and 264 for KOI-1783.01 and KOI-1783.02, respectively (the authors note that Δχ² ≳ 20 is a reliable detection threshold).

For non-dynamically interacting systems, it is common to use statistical arguments to establish that the planetary hypothesis is the most likely explanation for a given transit signal using codes such as the publicly available false-positive probability (FPP) calculator vespa (Morton 2012, 2015). The vespa package has been used to statistically validate more than a thousand exoplanet candidates from Kepler and K2 thus far (Crossfield et al. 2016; Morton et al. 2016; Livingston et al. 2018b, 2018a; Mayo et al. 2018), although refutation of some previously validated planets suggests that caution is necessary when validating with limited follow-up data (Santerne et al. 2016; Cabrera et al. 2017; Shporer et al. 2017). Morton et al. (2016) obtained FPPs for all KOIs, including KOI-1783.01 (FPP = 0.680 ± 0.014) and KOI-1783.02 (FPP = 0.200 ± 0.012). However, TTVs were not considered in the construction of the light curves for these planets, which can inflate the FPP by making the transits look more "V"-shaped. Additionally, Morton et al. (2016) found four confirmed planets with anomalously high FPPs: three exhibited TTVs, and the other had grazing transits. Our analysis suggests that KOI-1783 system is a near-grazing TTV system, making it very likely to have an overestimated FPP.

In a six-year campaign, Santerne et al. (2016) performed RV observations of a sample of 125 KOI stars, including KOI-1783. They observed KOI-1783 two times with SOPHIE and detected no RV variation. Additionally, they establish 99% upper limits on the RV semiamplitude (K < 81.3 m s⁻¹) and corresponding mass (M < 2.83M_J). While these upper limits were derived by fitting a circular orbit with no TTVs, the lack of detected RV variations rule out the eclipsing binary false-positive mode to very high confidence.

In addition to high-resolution spectroscopic follow-up, three ground-based adaptive optics (AO) follow-up observations of KOI-1783 have been performed to date, as listed by Furlan et al. (2017) and the Exoplanet Follow-up Observing Program. The Robo-AO team observed this star in their LP600 filter with the Palomar 60 inch telescope, achieving a contrast of ΔM = 4.00 mag at $0\buildrel{\prime\prime}\over{.} 30$ (Law et al. 2014). Additionally, Wang et al. (2015) observed KOI-1783 in the K_s band with PHARO on the Hale 200 inch telescope at Palomar Observatory, achieving a contrast of ΔM = 4.33 mag at $0\buildrel{\prime\prime}\over{.} 50$. More stringent contrast constraints of ΔM = 7.96 mag at $0\buildrel{\prime\prime}\over{.} 50$ were obtained with NIRC2 on the Keck II Telescope using the Brγ filter (Furlan et al. 2017). These observations demonstrate that there are no nearby stars that might explain the $0\buildrel{\prime\prime}\over{.} 25$ offset noted in the Data Validation Report.

Published RV data rule out the existence of an eclipsing binary, and AO imaging data rule out the existence of companions. Combined with the aforementioned multiple independent analyses all supporting dynamical interactions between the bodies in the system, these follow-up constraints lead us to conclude that the two transit candidates in the KOI-1783 system should be confirmed as bona fide planets.

5.2.1. TTVs Probe Warm Sub-Neptune-sized Planets

There are currently very few sub-Neptune-sized transiting planets with well-measured masses at large orbital distances (P > 100 days); these systems are quite rare to begin with, and most are too small and faint to be amenable to RV follow-up (Jontof-Hutter 2019). TTV studies that probe this regime are thus quite valuable, as planets that receive low incident fluxes are much more likely to retain their primordial atmospheres than their more highly irradiated counterparts (e.g., Owen & Wu 2013; Mazeh et al. 2016). Even if mass loss is common for these longer-period planets, the mechanism by which it occurs may be quite different. For highly irradiated exoplanets, atmospheric mass loss is primarily driven by thermal escape processes as the intense XUV flux heats the upper atmospheres (e.g., Owen 2019). However, for planets on more distant orbits, nonthermal processes are competitive with or dominant over photoevaporative escape; this is, for instance, the present case for terrestrial planets like Mars (Tian et al. 2013; Tian 2015). Density constraints for this population of long-period extrasolar planets at low ($\lesssim 100{F}_{\oplus }$) incident fluxes are therefore critical for building a holistic understanding of atmospheric mass loss in the regime relevant for potentially habitable terrestrial planets.

In Figure 4, we plot the masses and radii of our sub-Neptune-sized sample ($M\lt 17{M}_{\oplus }$) along with those from the NASA Exoplanet Archive and compare their radii to their incident fluxes. Other than the rocky super-Earth Kepler-36b (Carter et al. 2012), all of the planets in our sample are more inflated than they would be if they were purely composed of silicate rock (Fortney et al. 2007), implying that they possess at least modest volatile-rich envelopes. Even after allowing for water-rich compositions, our bulk density estimates for the planets in Table 6 are still too low and likely require a modest hydrogen-rich atmosphere. For Kepler-29b, Kepler-29c, Kepler-36c, and Kepler-177b, the grids of Lopez & Fortney (2014) suggest hydrogen-helium envelope fractions of 2%–5% in mass. For the more massive sub-Neptunes KOI-1783.02 and Kepler-177c, these grids suggest hydrogen-helium envelope fractions greater than 10% in mass. In the following section, we explore the bulk composition of KOI-1783.01, KOI-1783.02, and Kepler-177c in more detail.

Zoom In Zoom Out Reset image size

5.2.2. Bulk Metallicities of the Giant Planets KOI-1783.01, KOI-1783.02, and Kepler-177c

TTVs can also deliver masses and radii for giant planets in the low-insolation regime. This is crucial for estimates of bulk metallicity, as gas giants hotter than approximately 1000 K appear to have inflated radii that are inconsistent with predictions from standard interior models (e.g., Laughlin et al. 2011; Thorngren et al. 2016; Thorngren & Fortney 2018). Relatively cool, dynamically interacting planets such as KOI-1783.01 are not expected to be affected by this inflation mechanism and are therefore ideal candidates for these studies.

We measure the mass of the gas giant KOI-1783.01 to ∼15% precision and its radius to ∼3%, as this star has relatively accurate stellar parameters from Fulton & Petigura (2018). When combined with our incident flux constraints and stellar age estimates from Fulton & Petigura (2018), these parameters yield a bulk metallicity of Z_p = 0.30 ± 0.03 for KOI-1783.01 using the statistical model of Thorngren & Fortney (2019). Using the stellar metallicity from Table 2 and the ${Z}_{\mathrm{star}}=0.014\times {10}^{[\mathrm{Fe}/{\rm{H}}]}$ prescription from Thorngren et al. (2016), this corresponds to ${Z}_{{\rm{p}}}/{Z}_{\mathrm{star}}={16.6}_{-2.2}^{+2.4}$. We note that when masses and radii are constrained to this level of precision, we should also consider the additional uncertainties introduced by the choice of models, which are not accounted for in these error bars (Thorngren et al. 2016; Thorngren & Fortney 2019). This bulk metallicity value is nevertheless in excellent agreement with the mass–metallicity relation previously inferred for gas-giant planets at higher incident fluxes (Thorngren et al. 2016; Thorngren & Fortney 2019), as shown in Figure 5.

Zoom In Zoom Out Reset image size

This bulk metallicity also yields an upper limit on the atmospheric metallicity, as the metallicity observable in a planetary atmosphere will always be less than the total metal content of the planet (Thorngren & Fortney 2019). For KOI-1783.01, this (95th percentile) upper limit is Z_atm ≤ 79 × solar, where "solar" refers to the Asplund et al. (2009) photospheric metal fraction of 1.04 × 10⁻³. This calculation assumes an average mean molecular mass of 18 (that of water) for this heavy element component; if this is not the case, then the true upper limit on the atmospheric metallicity should be scaled by 18/μ_Z (Thorngren & Fortney 2019).

We calculate comparable bulk composition estimates for the two sub-Neptunes in our sample, KOI-1783.02 and Kepler-177c. In this mass regime, differences in equation of state between rock and water ice become important, adding another degree of freedom to the calculation. We construct models composed of a rock layer, a water layer, and low-density H/He layer enriched to Neptune's metallicity (90× solar) by borrowing water from the water layer. We do not include mass loss in our simulation, and we assume negligible amounts of iron in the calculation. We use constraints on the mass, radius, host star age, and incident flux to retrieve the composition, including the relative amounts of rock, water, and H/He. Although we are not able to place strong constraints on the relative amounts of rock versus water as the radius is still fairly insensitive to the core composition details (Lopez & Fortney 2014; Petigura et al. 2017b), we are able to place a strong constraint on the total bulk metallicity Z_p and the corresponding the H/He fraction ${f}_{{\rm{H}}/\mathrm{He}}=1-{Z}_{{\rm{p}}}$.

As hinted at by their low bulk densities, these two planets have large H/He mass fractions: ${f}_{{\rm{H}}/\mathrm{He}}=0.31\pm 0.08$ for KOI-1783.02 and ${f}_{{\rm{H}}/\mathrm{He}}=0.74\pm 0.04$ for Kepler-177c. The value for Kepler-177c is somewhat problematic from a planet formation perspective, as it implies a maximum core mass of just 4 ${M}_{\oplus }$. Depending on the planet's formation location, it may be difficult to explain how such a small core could have accreted such a massive gas envelope. One explanation is that the core formed outside 1 au and experienced relatively dust-free accretion, as is typically invoked for super-puffs (Lee & Chiang 2016). We note, however, that super-puffs are a few times less massive than Kepler-177c despite having similar inferred core masses, implying that the gas-to-core mass ratio (GCR) of Kepler-177c exceeds that of a typical super-puff. Although it is possible that our estimate of this maximum core mass might have been biased by assumptions made in our models, accounting for atmospheric mass loss would have preferentially removed hydrogen and helium, and including iron in the model would have increased the ${f}_{{\rm{H}}/\mathrm{He}}$. We conclude that these assumptions are unlikely to explain the large inferred H/He mass fraction for this planet. The MIST isochrone-derived age estimate for this planet from Fulton & Petigura (2018) appears to be quite secure, with $\mathrm{log}(\mathrm{age})=10.07\pm 0.04$, so it is unlikely that this planet's radius is inflated by residual heat from formation.

Can Kepler-177c be inflated by internal heating mechanisms such as Ohmic dissipation (Pu & Valencia 2017) or obliquity tides (Millholland 2019)? Its large total mass and low insolation makes this scenario unlikely. We assess the scenario of Kepler-177c having a core mass of $14.5{M}_{\oplus }$ and an envelope mass of $0.2{M}_{\oplus }$ (envelope mass fraction of 1%). Its estimated equilibrium temperature is ∼800K, too low for Ohmic dissipation to puff up Kepler-177c to ≳8${R}_{\oplus }$ (see Figures 8 and 9 of Pu & Valencia 2017). Next, we assess heating by obliquity tides. Even if we assume maximal obliquity, the expected thickness of the envelope is ∼0.48${R}_{\oplus }$ (see Equation (13) of Millholland 2019). If the composition of the Kepler-177c core is similar to that of Earth, we expect its core size to be ∼1.95${R}_{\oplus }$ (assuming $R\propto {M}^{1/4}$), so that the expected total radii of the planet is only ∼2.43${R}_{\oplus }$, far too small to explain the measured $8.73{R}_{\oplus }$. Even at a GCR of 10%, the expected total radii is just $3.74{R}_{\oplus }$.

We conclude that Kepler-177c rightfully belongs in the small sample of $\sim 15{M}_{\oplus }$ planets with extremely low bulk densities (and thus extremely large envelope fractions). This sample also includes Kepler-18d (Cochran et al. 2011; Petigura et al. 2017b) and K2-24c (Petigura et al. 2018). Petigura et al. (2018) suggest a formation scenario for the latter planet wherein the disk dissipates just as the planet begins to enter runaway accretion. Lee (2019) showed that the sub-Saturn population can indeed be explained by the timing of disk dispersal, but they note as a prerequisite that their cores must be massive enough to trigger runaway accretion during the disk lifetime, $\gtrsim 10{M}_{\oplus }$. For cores less massive than this, the maximum GCR is set by the amount of gas that can be accreted by cooling. In Figure 6, we reproduce the Lee (2019) GCR plot as a function of core mass and accretion time, which highlights the different regimes dictating the maximum envelope fraction for a given core mass. While KOI-1783.01 and KOI-1783.02 can largely be explained within the framework of disk dispersal timing relative to the onset of runaway accretion, Kepler-177c cannot, nor can K2-24c or Kepler-18d. These low-density $15{M}_{\oplus }$ planets are outliers, lying above their theoretical maximum GCRs, as are the super-puffs Kepler-51b (Masuda 2014), Kepler-223e (Mills et al. 2016), Kepler-87c (Ofir et al. 2014), and Kepler-79d (Jontof-Hutter et al. 2014).

Zoom In Zoom Out Reset image size

As a result, Lee (2019) suggests that these more massive low-density planets may share a formation pathway with the less-massive super-puffs. Super-puffs likely accreted their envelopes farther from their star and then migrated inwards (Ikoma & Hori 2012; Lee et al. 2014; Ginzburg et al. 2016; Lee & Chiang 2016; Schlichting 2018), and additionally should have experienced "dust-free" accretion, meaning that dust did not contribute much to the overall opacity due to, e.g., grain growth or sedimentation (Lee & Chiang 2015, 2016). To test the feasibility of this hypothesis, we can estimate the amount of time that Kepler-177c must have spent undergoing dust-free accretion and compare to typical disk lifetimes. If this timescale is longer than the typical disk dispersal timescale, then a mechanism other than dust-free accretion is necessary; if it is comparable or shorter, then dust-free accretion may be feasible. For Kepler-177c (${M}_{\mathrm{core}}\approx 3.8{M}_{\oplus },\mathrm{GCR}\approx 2.8$), we can approximate the dust-free accretion time necessary to achieve the observed GCR beyond 1 au in a gas-rich disk using the analytic scaling relation of Lee & Chiang (2015, see their Equation (24)):

$$\begin{eqnarray}&&t\sim 1\ \mathrm{kyr}{\left[\left(\displaystyle \frac{\mathrm{GCR}}{0.1}\right)\left(\displaystyle \frac{5{M}_{\oplus }}{{M}_{\mathrm{core}}}\right)\right]}^{2.5}\approx 8.2\ \mathrm{Myr},\end{eqnarray} \tag{ 12 }$$

where for simplicity we have assumed their nominal values for the f factor, the nebular gas metallicity Z, the adiabatic gradient ∇_ad, and the temperature and mean molecular weight at the radiative–convective boundary T_rcb = 200 K and μ_rcb. The outer layers of dust-free envelopes are largely isothermal so the adopted temperature corresponds to the nebular temperature at the formation location. The estimated accretion timescale required to build Kepler-177c is comparable to typical disk lifetimes (∼5 Myr; see, e.g., Alexander et al. 2014 and references therein). We note that Equation (12) is derived assuming the self-gravity of the envelope is negligible compared to the gravity of the core. The rate of accretion starts to accelerate once GCR ≳ 0.5, so a more careful calculation would provide an even shorter timescale. We suggest that $15{M}_{\oplus }$ planets with large GCRs may indeed share a dust-free accretion history with their lower-mass super-puff counterparts. As such, detailed characterization of Neptune-mass planets with low (ρ ≲ 0.3 g cm⁻³) bulk densities may provide invaluable insights into super-puff formation processes.