Diffuser-Assisted Infrared Transit Photometry for Four Dynamically Interacting \textit{Kepler} Systems

TL;DR

This study demonstrates that ground-based infrared transit observations, aided by an engineered diffuser, can achieve space-like photometric precision, enabling improved mass and radius measurements for exoplanets in dynamically interacting systems.

Contribution

The paper introduces a diffuser-assisted ground-based infrared photometry technique that rivals space-based precision, reducing uncertainties in exoplanet mass and radius measurements.

Findings

Achieved infrared photometric precision comparable to space telescopes.

Detected transits with >3σ significance for all four planets.

Reduced dynamical mass uncertainties by up to a factor of three.

Abstract

We present ground-based infrared transit observations for four dynamically interacting \textit{Kepler} planets, including Kepler-29b, Kepler-36c, KOI-1783.01, and Kepler-177c, obtained using the Wide-field Infrared Camera on the Hale 200" telescope at Palomar Observatory. By utilizing an engineered diffuser and custom guiding software, we mitigate time-correlated telluric and instrumental noise sources in these observations. We achieve an infrared photometric precision comparable to or better than that of space-based observatories such as the \textit{Spitzer Space Telescope}, and detect transits with greater than 3$\sigma$ significance for all planets. For Kepler-177c ($J=13.9$) our measurement uncertainties are only $1.2\times$ the photon noise limit and 1.9 times better than the predicted photometric precision for \textit{Spitzer} IRAC photometry of this same target. We find that a single transit observation obtained $4-5$ years after the end of the original \textit{Kepler} mission can reduce dynamical mass uncertainties by as much as a factor of three for these systems. Additionally, we combine our new observations of KOI-1783.01 with information from the literature to confirm the planetary nature of this system. We discuss the implications of our new mass and radius constraints in the context of known exoplanets with low incident fluxes, and we note that Kepler-177c may be a more massive analog to the currently known super-puffs given its core mass (3.8$\pm0.9M_\Earth$) and large gas-to-core ratio (2.8$\pm0.7$). Our demonstrated infrared photometric performance opens up new avenues for ground-based observations of transiting exoplanets previously thought to be restricted to space-based investigation.

Tables

Table 1. Table 1: Observational parameters for our four nights of data collection.

Star	J magaaJ band magnitudes from the 2MASS catalogue (Cutri et al., 2003).	Date	Start Time	End Time	Event TimebbPredicted mid-transit time.	Event Duration	Start/Min/End Airmass	Exposure Time
		(UTC)	(UTC)	(UTC)	(UTC)	(hr)		(s)
Kepler-29	14.13	2017 August 25	05:35:24	11:57:00	08:26:53	3.046	1.03/1.03/3.01	25string
Kepler-36	11.12	2017 September 27	03:06:20	08:55:42	09:52:34	7.461	1.04/1.04/2.50	16dd4 co-adds of 4 second exposures.
KOI-1783	12.92	2018 April 21	08:19:42	12:04:05	07:07:51	5.871	1.73/1.05/1.05	20string
Kepler-177	13.86	2018 May 4	07:17:36	12:09:04	10:30:49	5.245	1.73/1.02/1.02	75ee3 co-adds of 25 second exposures.

Table 2. Table 2: Stellar parameters for the stars in our sample.

Target	$T_{\mathrm{eff}}$	[Fe/H]	$\log \left(g\right)$	$M_{\star }$	$R_{\star }$
	(K)	(dex)	(log(cm/s2))	($M_{☉}$)	($R_{☉}$)
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}$

Table 3. Table 3: Photometric quality statistics for the observations presented in this work.

Planet	WIRC Transit Coverage	Kepler RMS	WIRC RMS	WIRC RMS	WIRC Binned RMS	$\log \left({\sigma }_{\mathrm{extra},\mathrm{WIRC}}\right)$
	(%)	(ppm)	(ppm)	($\times$ photon noise)	($\times$ photon noise)
Kepler-29b	100	504	4222	1.20	1.27	-2.627
Kepler-36c	41.8	75	1305	2.10	2.46	-2.943
KOI-1783.01	33.7	157	2862	1.48	1.29	-2.680
Kepler-177c	66.9	320	2403	1.22	1.46	-2.851

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	d	(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	string0.4959	string0.4639	string0.6034	string0.5716	—	(3)
	a_2	string0.0222	string0.3045	-0.1382	-0.1145	—	(3)
	a_3	string0.5708	string0.0751	string0.6330	string0.6579	—	(3)
	a_4	-0.3485	-0.1251	-0.3506	-0.3667	—	(3)
WIRC limb darkening coefficients	b_1	string0.3634	string0.3982	string0.4832	string0.4421	—	(4)
	b_2	string0.5846	string0.5452	string0.2998	string0.3993	—	(4)
	b_3	-0.6152	-0.6817	-0.3634	-0.4523	—	(4)
	b_4	string0.1997	string0.2508	string0.1152	string0.1474	—	(4)
Fit Priors
Transit depth prior	P_(R_p/R_⋆)^2	$𝒰(0,2000)$	$𝒰(0,1000)$	$𝒰(0,10000)$	$𝒰(0,8000)$	ppm	—
Transit timing offset prior	P_Δt_0	$𝒰(-100,100)$	$𝒰(-100,100)$	$𝒰(-100,100)$	$𝒰(-100,100)$	min	—
Inclination prior	P_i	$𝒰(85,90)$	$𝒰(85,90)$	$𝒰(85,90)$	$𝒰(85,90)$	°	—
Scaled semi-major axis prior	P_a/R_⋆	$𝒩(24.906,1.125)$	$𝒩(16.696,0.436)$	$𝒩(99.030,2.840)$	$𝒩(41.649,1.674)$	—	(5)
Fit Posteriors
Transit depth	(R_p/R_⋆)^2	1020${}_{-34}{}^{}{}_{}{}^{+31}$	425.3${}_{-3.5}{}^{}{}_{}{}^{+3.8}$	5044${}_{-64}{}^{}{}_{}{}^{+87}$	3643${}_{-57}{}^{}{}_{}{}^{+55}$	ppm	—
Transit timing offset	Δ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 semi-major axis	a/R_⋆	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_p/R_⋆	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	–

Table 5. Table 5: Results from our dynamical analysis.

Planet	Dataset	$P$ [days]	$\left(\frac{M_p}{M_{\oplus }}\right)\left(\frac{M_{☉}}{M_{\star }}\right)$	$e\cos \left(\omega \right)$	$e\sin \left(\omega \right)$
Kepler-29b	Kep LC	${10.33838}_{-0.00027}^{+0.00030}$	string4.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}$	string3.7${}_{-1.3}{}^{}{}_{}{}^{+1.3}$	string0.013${}_{-0.071}{}^{}{}_{}{}^{+0.071}$	-0.016${}_{-0.063}{}^{}{}_{}{}^{+0.056}$
	Kep LC + WIRC (S)	${10.33966}_{-0.00017}^{+0.00015}$	string3.8${}_{-1.0}{}^{}{}_{}{}^{+1.1}$	string0.003${}_{-0.070}{}^{}{}_{}{}^{+0.068}$	-0.088${}_{-0.058}{}^{}{}_{}{}^{+0.059}$
Kepler-29c	Kep LC	${13.28843}_{-0.00053}^{+0.00048}$	string4.07${}_{-2.29}{}^{}{}_{}{}^{+2.87}$	string0.007${}_{-0.062}{}^{}{}_{}{}^{+0.063}$	-0.022${}_{-0.063}{}^{}{}_{}{}^{+0.063}$
	Kep LC + WIRC (G)	${13.28613}_{-0.00021}^{+0.00026}$	string3.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}$	string3.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}$	string3.990${}_{-0.092}{}^{}{}_{}{}^{+0.093}$	string0.050${}_{-0.025}{}^{}{}_{}{}^{+0.023}$	-0.026${}_{-0.033}{}^{}{}_{}{}^{+0.034}$
	Kep LC + WIRC (G)	${13.86825}_{-0.00050}^{+0.00050}$	string3.972${}_{-0.074}{}^{}{}_{}{}^{+0.078}$	string0.041${}_{-0.020}{}^{}{}_{}{}^{+0.019}$	-0.011${}_{-0.018}{}^{}{}_{}{}^{+0.018}$
	Kep LC + WIRC (S)	${13.86821}_{-0.00049}^{+0.00049}$	string3.964${}_{-0.068}{}^{}{}_{}{}^{+0.077}$	string0.037${}_{-0.018}{}^{}{}_{}{}^{+0.019}$	-0.004${}_{-0.015}{}^{}{}_{}{}^{+0.012}$
Kepler-36c	Kep LC	${16.21867}_{-0.00010}^{+0.00010}$	string7.456${}_{-0.168}{}^{}{}_{}{}^{+0.167}$	string0.053${}_{-0.023}{}^{}{}_{}{}^{+0.021}$	-0.039${}_{-0.031}{}^{}{}_{}{}^{+0.031}$
	Kep LC + WIRC (G)	${16.21865}_{-0.00010}^{+0.00010}$	string7.397${}_{-0.107}{}^{}{}_{}{}^{+0.104}$	string0.046${}_{-0.018}{}^{}{}_{}{}^{+0.017}$	-0.026${}_{-0.017}{}^{}{}_{}{}^{+0.017}$
	Kep LC + WIRC (S)	${16.21865}_{-0.00010}^{+0.00010}$	string7.371${}_{-0.093}{}^{}{}_{}{}^{+0.092}$	string0.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}$	string0.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}$	string0.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}$	string0.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}$	string0.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}$	string0.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}$	string0.017${}_{-0.014}{}^{}{}_{}{}^{+0.015}$	-0.020${}_{-0.026}{}^{}{}_{}{}^{+0.034}$
Kepler-177b	Kep LC	${35.8591}_{-0.0017}^{+0.0019}$	string5.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}$	string5.44${}_{-0.75}{}^{}{}_{}{}^{+0.78}$	string0.017${}_{-0.054}{}^{}{}_{}{}^{+0.052}$	-0.001${}_{-0.063}{}^{}{}_{}{}^{+0.062}$
	Kep LC + WIRC (S)	${35.8601}_{-0.0012}^{+0.0013}$	string5.38${}_{-0.74}{}^{}{}_{}{}^{+0.78}$	string0.020${}_{-0.048}{}^{}{}_{}{}^{+0.047}$	string0.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}$	string0.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}$	string0.013${}_{-0.041}{}^{}{}_{}{}^{+0.040}$	string0.003${}_{-0.053}{}^{}{}_{}{}^{+0.052}$

Table 6. Table 6: Physical parameters for the planets in this study.

Planet	$M_p$ [$M_{\oplus }$]aaCalculated from our dynamical masses and the stellar masses of Fulton & Petigura (2018).	$R_p$ [$R_{\oplus }$]bbCalculated from either our measured $R_p/R_{\star }$ or that from Thompson et al. (2018) and stellar radii from Fulton & Petigura (2018).	${\rho }_p$ [g/cm3]	$F_{\mathrm{in}}\left[F_{\oplus }\right]$ccCalculated in the low-eccentricity () approximation via $F_{\mathrm{in}}=4.62\times {10}^4{F}_{\oplus }{\left(\frac{T_{\mathrm{eff}}}{T_{☉}}\right)}^4{\left(\frac{a}{R_{\star }}\right)}^{-2}$ (Jontof-Hutter et al., 2016), with effective temperatures from Fulton et al. (2017) and scaled semi-major axes from our measurements or Thompson et al. (2018).
Kepler-29b	string5.0${}_{-1.3}{}^{}{}_{}{}^{+1.5}$	2.55${}_{-0.12}{}^{}{}_{}{}^{+0.12}$	1.65${}_{-0.49}{}^{}{}_{}{}^{+0.53}$	55.9${}_{-4.8}{}^{}{}_{}{}^{+6.5}$
Kepler-29cddRadius ratio and scaled semi-major axis taken from Thompson et al. (2018).	string4.5${}_{-1.1}{}^{}{}_{}{}^{+1.1}$	2.34${}_{-0.11}{}^{}{}_{}{}^{+0.12}$	1.91${}_{-0.54}{}^{}{}_{}{}^{+0.57}$	34.4${}_{-3.8}{}^{}{}_{}{}^{+3.8}$
Kepler-36bddRadius ratio and scaled semi-major axis taken from Thompson et al. (2018).	string3.83${}_{-0.10}{}^{}{}_{}{}^{+0.11}$	1.498${}_{-0.049}{}^{}{}_{}{}^{+0.061}$	6.26${}_{-0.64}{}^{}{}_{}{}^{+0.79}$	247${}_{-32}{}^{}{}_{}{}^{+32}$
Kepler-36c	string7.13${}_{-0.18}{}^{}{}_{}{}^{+0.18}$	3.679${}_{-0.091}{}^{}{}_{}{}^{+0.096}$	0.787${}_{-0.062}{}^{}{}_{}{}^{+0.065}$	191.0${}_{-10.4}{}^{}{}_{}{}^{+9.7}$
KOI-1783.01	71.0${}_{-9.2}{}^{}{}_{}{}^{+11.2}$	8.86${}_{-0.24}{}^{}{}_{}{}^{+0.25}$	0.560${}_{-0.085}{}^{}{}_{}{}^{+0.101}$	string5.70${}_{-0.27}{}^{}{}_{}{}^{+0.27}$
KOI-1783.02ddRadius ratio and scaled semi-major axis taken from Thompson et al. (2018).	15.0${}_{-3.6}{}^{}{}_{}{}^{+4.3}$	5.44${}_{-0.30}{}^{}{}_{}{}^{+0.52}$	0.51${}_{-0.15}{}^{}{}_{}{}^{+0.21}$	string2.49${}_{-0.35}{}^{}{}_{}{}^{+0.35}$
Kepler-177bddRadius ratio and scaled semi-major axis taken from Thompson et al. (2018).	string5.84${}_{-0.82}{}^{}{}_{}{}^{+0.86}$	3.50${}_{-0.15}{}^{}{}_{}{}^{+0.19}$	0.75${}_{-0.14}{}^{}{}_{}{}^{+0.16}$	30.4${}_{-4.0}{}^{}{}_{}{}^{+4.0}$
Kepler-177c	14.7${}_{-2.5}{}^{}{}_{}{}^{+2.7}$	8.73${}_{-0.34}{}^{}{}_{}{}^{+0.36}$	0.121${}_{-0.025}{}^{}{}_{}{}^{+0.027}$	25.4${}_{-1.6}{}^{}{}_{}{}^{+1.6}$

Equations

$$\vec{M}=[\alpha \vec{t}+\vec{S}]\times {\vec{T}}_{\mathrm{WIRC}}((R_{\mathrm{p}}/R_{\star }{)}^2,t_0,i,a/R_{\star }),$$

$${\vec{S}}_{\mathrm{target}}=\frac{\vec{y}}{{\vec{T}}_{\mathrm{WIRC}}((R_{\mathrm{p}}/R_{\star }{)}^2,t_0,i,a/R_{\star })}-\alpha \vec{t},$$

$$\vec{S}=\mathbf{P}{\vec{S}}_{\mathrm{target}},$$

$$\mathbf{P}={\mathbf{X}}^T(\mathbf{X}{\mathbf{X}}^T{)}^{-1}\mathbf{X}$$

$${\vec{M}}_{Kep}={\vec{T}}_{Kep}((R_{\mathrm{p}}/R_{\star }{)}^2,i,a/R_{\star })$$

$$\log (\mathcal{L})=$$

$$-\frac{1}{2}\sum _i\log (2\pi {\sigma }_{Kep,i}^2)$$

$$\displaystyle-\frac{1}{2}\sum_{i}\Big{(}\frac{y_{Kep,i}-M_{Kep,i}}{\sigma_{Kep,i}}\Big{)}^{2},$$

$$\vec{\sigma }$$

$${\vec{\sigma }}_{Kep}$$

$$\log(\mathcal{P}_{k})=-\frac{1}{2}\log(2\pi\sigma_{k}^{2})-\frac{1}{2}\Big{(}\frac{k-\mu_{k}}{\sigma_{k}}\Big{)}^{2},$$

$$\log ({\mathcal{P}}_k)=$$

$$-\infty \mathrm{otherwise}.$$

$$\log (\mathrm{Prob})=\log (\mathcal{L})+\sum _k\log ({\mathcal{P}}_k)$$

$$t\sim\mathrm{1\ kyr}\Bigg{[}\Big{(}\frac{\mathrm{GCR}}{0.1}\Big{)}\Big{(}\frac{5M_{\earth}}{M_{\mathrm{core}}}\Big{)}\Bigg{]}^{2.5}\approx 8.2\ \mathrm{Myr},$$

Full text

techniques: photometric – planets and satellites: fundamental parameters – planets and satellites: individual (Kepler-29b, Kepler-36c, KOI-1783.01, Kepler-177c)

††facilities: Hale (WIRC, PHARO), Kepler, OHP:1.52m (SOPHIE), PO:1.5m (ROBO-AO), Keck:II (NIRC2), ADS, Exoplanet Archive††software: photutils (Bradley et al., 2016), numpy (van der Walt et al., 2011), astropy (Astropy Collaboration et al., 2013, 2018), scipy (Jones et al., 2001), matplotlib (Hunter, 2007), batman (Kreidberg, 2015), emcee (Foreman-Mackey et al., 2013), corner (Foreman-Mackey, 2016), PyKE (Still & Barclay, 2012), Aladin Lite (Bonnarel et al., 2000; Boch & Fernique, 2014)

1 Introduction

The Kepler mission (Borucki et al., 2010) has revealed thousands of transiting exoplanets and exoplanet candidates over the past decade, many of which reside in multi-planet systems. Dynamical interactions between planets in these systems cause deviations from the expected Keplerian behavior that can change both the timing and duration of transits (Agol et al., 2005; Holman & Murray, 2005; Agol & Fabrycky, 2017). In systems where planetary periods are close to integer multiples of each other – in other words, for planets close to or occupying mean motion resonances – the amplitude of transit timing variations (TTVs) and transit duration variations (TDVs) may become observable and reveal the dynamical architecture of the system. Approximately 10% of Kepler Objects of Interest (KOIs) exhibit significant long-term TTVs (Holczer et al., 2016). Most of these planets are on $\lesssim 100$ day orbits, with eccentricities of a few percent and sizes ranging from 1-10 $R_{\earth}$ (Holczer et al., 2016; Hadden & Lithwick, 2017).

TTV analyses have yielded a wealth of information about the properties of Kepler multi-planet systems, but arguably their most valuable contribution to date has been estimates of planet masses and densities for systems that are not amenable to characterization using the radial velocity (RV) technique (e.g. Wu & Lithwick, 2013; Jontof-Hutter et al., 2016; Hadden & Lithwick, 2017). These density constraints are especially critical for interpreting the bimodal radius distribution observed for close-in planets, which peaks at approximately 1.3 and 2.5 $R_{\earth}$ (Fulton et al., 2017; Fulton & Petigura, 2018). It has been suggested that this distribution is well-matched by models in which a subset of highly irradiated rocky planets have lost their primordial atmospheres while more distant planets retain modest (few percent in mass) hydrogen-rich atmosphere that inflate their observed radii (Owen & Wu, 2013; Lopez & Fortney, 2013, 2014; Fulton et al., 2017; Owen & Wu, 2017; Fulton & Petigura, 2018). Measuring the bulk density of planets in this size regime is thus a direct test of these photoevaporative models.

In Figure 1, we plot all confirmed Kepler planets (with those exhibiting TTVs specially marked) on the radius-period plane, following Fulton et al. (2017). In general, the TTV sample allows for characterization of planets that are $1.75R_{\earth}$ and larger (on the sub-Neptune side of the bimodal radius distribution), with periods longer than a week. While the radial velocity technique is most sensitive to short-period planets with relatively high densities, TTV observations are well-suited to characterizing long period and/or low-density planets, making it an important tool for probing this region of parameter space (Steffen, 2016; Mills & Mazeh, 2017). Indeed, this technique has already revealed the existence of a separate sub-population of “super-puffs,” a rare class of super-Earths with very low bulk densities and relatively long orbital periods (Masuda, 2014; Jontof-Hutter et al., 2014). Unlike the broader super-Earth population, which some studies argue could have formed in situ, it is thought that these planets may have accreted their envelopes at large stellocentric distances and then migrated inward to their current locations in resonant chains (Ikoma & Hori, 2012; Lee et al., 2014; Ginzburg et al., 2016; Lee & Chiang, 2016; Schlichting, 2018).

These previous studies showcase the crucial role of Kepler TTVs in testing theories of planet formation and evolution. The failure of Kepler’s second reaction wheel in 2013, however, effectively limited the baseline of these TTV analyses to four years. This makes it particularly challenging to constrain masses and bulk densities for long-period planets with a relatively small set of measured transits during this four-year period. In addition, uncertainties in the orbital solutions grow over time, making future in-transit observations (for instance, those aimed at atmospheric characterization) increasingly difficult to schedule with confidence.

These problems can be ameliorated with ground- or space-based follow-up observations (Petigura et al., 2018; Wang et al., 2018). However, many of the Kepler planets exhibiting TTVs orbit faint ($V>12$) stars, making it difficult to achieve the required photometric precision using existing space-based facilities with small apertures, such as the Spitzer Space Telescope. Additionally, Spitzer will be decommissioned in January 2020, necessitating an alternative approach to follow-up observations. Although ongoing observations by the Transiting Exoplanet Survey Satellite (TESS; Ricker et al., 2015) are expected to recover a few hundred Kepler planets (Christ et al., 2018), short-cadence data from the nominal mission will only improve the mass uncertainties for 6-14 of the $\sim$150 currently known Kepler TTV planets (Goldberg et al., 2019). This is due to the limited photometric precision and relatively short baseline of TESS relative to Kepler. While TESS is expected to recover additional transits in an extended mission scenario, these detections will still constitute less than 20% of the overall Kepler TTV sample (Goldberg et al., 2019).

Ground-based observatories can in principle recover transits for faint Kepler stars with long period planets, and coordinated multi-observatory campaigns have shown promise in achieving the requisite phase coverage (Freudenthal et al., 2018; von Essen et al., 2018; Wang et al., 2018). However, their photometric precisions are typically limited by low observing efficiencies and the presence of time-correlated noise due to imperfect guiding and point-spread function (PSF) variations (Zhao et al., 2014; Croll et al., 2015; Stefansson et al., 2017). These difficulties can be mitigated by using diffusers to control the shape of the point spread function (PSF) and spread out light from the star over a larger area. Diffusers have already been installed on several ground-based telescopes and have been shown to achieve significantly better photometric precision than more traditional observing techniques (Stefansson et al., 2017, 2018; von Essen et al., 2019).

Here, we present diffuser-assisted TTV follow-up observations of four Kepler planets in dynamically interacting systems. We discuss our sample selection methodology and our observations of the four-planet sample with the Wide-field InfraRed Camera (WIRC; Wilson et al., 2003) in Section 2. In Section 3, we describe our image calibration, data reduction, light curve modeling, and dynamical modeling methods. We then present our results for each system in Section 4, along with some brief comments on the general performance of our instrument. In Section 5, we discuss some of the scientific implications of our new dynamical mass constraints within the broader exoplanet population, and we conclude with a summary of our results and a look towards future possibilities in Section 6.

2 Observations

2.1 Sample Selection

In this study we focused on the set of multi-planet 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 transit timing variations. 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 (SNR) in the limit of circular orbits. For systems exhibiting TTVs with high SNRs, 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). Since 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\sigma$ 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\sigma$ range in predicted transit times.

Within this sample of systems, we searched for targets with an ingress and/or egress visible from Palomar between August 2017 and May 2018. We then ranked the targets in our sample based on predicted signal-to-noise ratio (SNR) scaled from early WIRC commissioning data (Stefansson et al., 2017), and prioritized observations of the highest SNR 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 mid-transit times for these planets are shown in Table 1.

2.2 New WIRC Observations

We observed our four selected systems in J band with WIRC, which is located at the prime focus of the Hale 200” telescope at Palomar Observatory (Wilson et al., 2003). The current 2048 $\times$ 2048 pixel Hawaii-II HgCdTe detector was installed in January 2017, along with 32-channel readout electronics that allow for a read time of 0.92 s (Tinyanont et al., 2019). The instrument has an 8$\farcm$7 $\times$ 8$\farcm$7 field of view with a pixel scale of 0$\farcs$2487, ensuring that (at least for the magnitude range in our sample) there are always on the order of ten 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 a full width at half maximum (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 $\sim 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 co-added during the night to increase observing efficiency as noted in Table 1.

3 Data Reduction and Model Fits

3.1 Image Calibration and Photometry

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 $\times$ 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 (BJDTDB), 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\farcs 25$/pixel scale, the diffuser is expected to deliver stellar PSFs with a 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$\sigma$ clipping with the sigma_clipped_stats function in astropy with a five-iteration maximum specified (Astropy Collaboration et al., 2013, 2018). After this, we re-calculate 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\sigma$ 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 ten comparison stars that best track the time-varying flux of the target star (i.e. those that have the 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 dataset and removing 3$\sigma$ 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.

3.2 Kepler Light Curves

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 mid-transit time and transit duration (parameterized here by the inclination and semi-major 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_{\mathrm{p}}/R_{\star})^{2}$, the inclination $i$, and the scaled semi-major axis $a/R_{\star}$. 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_{\mathrm{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 (SAP) 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 year baseline. The WIRC data alone are not sensitive to transit duration changes on this timescale, since 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.

3.3 Light Curve Fitting

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 $\vec{t}=(t_{1},t_{2},...,t_{n})$, the target data array $\vec{y}=(y_{1},y_{2},...,y_{n})$ (with $y_{i}$ referring to the measurement at time $t_{i}$), and comparison star arrays $\vec{x}_{j}=(x_{1},x_{2},...,x_{n})$. Collectively, the comparison stars define a matrix $\mathbf{X}$, with one comparison star $\vec{x}_{j}$ in each row of the matrix.

We aim to fit the target $\vec{y}$ with a model $\vec{M}$ that depends on the depth of the transit $(R_{\mathrm{p}}/R_{\star})^{2}$, the transit center time $t_{0}$, the inclination $i$, the ratio of semi-major axis to stellar radius $a/R_{\star}$, and a linear trend in time $\alpha$. That model can be written as follows (loosely following the notation of Diamond-Lowe et al., 2018):

[TABLE]

where $\vec{S}$ is the systematics model, $\vec{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 multi-planet Kepler systems are typically small, with a population mean of $\bar{e}=0.04^{+0.03}_{-0.04}$ (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.

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” $\vec{S}_{\mathrm{target}}$). We calculate new linear coefficients every time the transit light curve is modified. Mathematically, the target systematics can be written:

[TABLE]

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

[TABLE]

where the projection matrix $\mathbf{P}$ comes from the comparison stars and can be written:

[TABLE]

Equations (1)–(4) thus define the model $\vec{M}$ solely as a function of the observed quantities {$\vec{t},\vec{y},\mathbf{X}$} and the fit parameters {$(R_{\mathrm{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.

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

[TABLE]

We supersampled the Kepler light curves to 1 min 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 §4.1) such that the likelihood takes the form:

[TABLE]

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

[TABLE]

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

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

[TABLE]

or uniform, taking the functional form:

[TABLE]

We placed physically motivated Gaussian priors on $a/R_{\star}$ 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.

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

[TABLE]

Then, we seek a solution for the fit parameters $(R_{\mathrm{p}}/R_{\star})^{2}$, $\Delta t_{0}$, $i$, $a/R_{\star}$, $\alpha$, $\log(\sigma_{\mathrm{extra,WIRC}})$, and $\log(\sigma_{\mathrm{extra},Kep})$ that maximizes $\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\times 10^{3}$ steps and then run for 105 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$, $\log(\sigma_{\mathrm{extra,WIRC}})$, and $\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.

3.4 Dynamical Modeling

Our fits to the ground-based WIRC photometry typically resulted in a non-Gaussian posterior for the mid-transit time. We accounted for these skewed distributions in our dynamical fits by dividing the posteriors into twenty 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 multi-planet 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 = 2455680), 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.

4 Results

We determine the significance of each detection in the WIRC data by re-running the joint fit and allowing the WIRC transit depth to vary independent of the Kepler transit depth. The confidence is then estimated using the width of the posterior on the WIRC transit depth. We detect transit signals for all four of our targets with $3\sigma$ or greater confidence in the WIRC data alone.

We show various quality statistics for each night of photometry in Table 3 (see Section 4.1 for additional details). Our results for the photometric fits to each observed planet are given in Table 4, and the resulting orbital periods, masses, and eccentricity vectors are presented in Table 5. We combine our photometric and dynamical results with previously computed stellar parameters to yield the physical planet parameters we report in Table 6. Below we discuss WIRC’s overall photometric performance as well as results for each individual system.

4.1 Instrument Performance

Our best photometric performance is for Kepler-177c, where we were only $\sim 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 sub-optimal co-addition 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 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 host star $J$ magnitude in Figure 3 along with the $3\sigma$ detection thresholds of WIRC and Spitzer. While Spitzer performs better for brighter stars, WIRC begins to out-perform Spitzer for stars fainter than $\sim 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.

4.2 Kepler-29

Kepler-29b is 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\sigma$ 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^{+17}_{-3}$ min.

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 three 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$\pm 0.5$ and 1.9$\pm 0.5$ g/cm3, 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.024}_{-0.028}~{}M_{\sun}$ and $R=0.732^{+0.033}_{-0.031}~{}R_{\sun}$ from Fulton & Petigura (2018), which are smaller than the values of $M=0.979\pm 0.052~{}M_{\sun}$ and $R=0.932\pm 0.060~{}M_{\sun}$ 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 §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 multi-planet Kepler systems tend to host planets that are similar in both size and bulk density (Millholland et al., 2017; Weiss et al., 2018).

4.3 Kepler-36

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$\sigma$. 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${}^{+12}_{-5}$ 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\times$ 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 $\sim 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 co-adds, resulting in a longer effective integration time and higher overall observing efficiency.

4.4 KOI-1783

As we will discuss in §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\sigma$ 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 two. 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_{\mathrm{p}}=8.9^{+0.3}_{-0.2}R_{\earth}$ and $M_{\mathrm{p}}=71^{+11}_{-9}M_{\earth}$. This corresponds to a density of $\rho=0.56^{+0.10}_{-0.09}$ g/cm3, consistent with the presence of a substantial gaseous envelope; we discuss the corresponding implications for this planet’s bulk composition in more detail in §5.2.2. KOI 1783.02 has a mass of $M_{\mathrm{p}}=15^{+4}_{-4}M_{\earth}$, a radius of $R_{\mathrm{p}}=5.4^{+0.5}_{-0.3}R_{\earth}$, and a density of $\rho=0.5^{+0.2}_{-0.2}$g/cm3, 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 hours 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.

4.5 Kepler-177

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 re-analyzed 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\sigma$ significance and measure the corresponding transit time with a $1\sigma$ uncertainty of approximately 10 minutes. Although our new dynamical fits for this system 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 re-converge 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.

5 Discussion

5.1 Confirmation of the KOI-1783 System

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; Lintott et al., 2013; Lissauer et al., 2014; 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_{\mathrm{p}}/R_{\star}$ is 0.9949${}^{+0.0032}_{-0.0027}$ for KOI-1783.01 from our posteriors, and 0.932${}^{+0.065}_{-0.015}$ 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\arcsec.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 (Steffen et al., 2013; Nesvorný 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$\sigma$ for KOI-1783.02. The spectral TTV analysis of Ofir et al. (2018) also found evidence of dynamical interactions, yielding $\Delta\chi^{2}$ values for the TTV signals over a linear model of 49 and 264 for KOI-1783.01 and .02, respectively (the authors note that $\Delta\chi^{2}\gtrsim 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., 2018a, b; 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\pm 0.014$) and KOI-1783.02 (FPP = $0.200\pm 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 semi-amplitude ($K<81.3$ m/s) and corresponding mass ($M<2.83\ M_{\mathrm{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” telescope, achieving a contrast of $\Delta M=4.00$ mag at $0\arcsec.30$ (Law et al., 2014). Additionally, Wang et al. (2015) observed KOI-1783 in $K_{s}$ band with PHARO on the Hale 200” telescope at Palomar Observatory, achieving a contrast of $\Delta M=4.33$ mag at $0\arcsec.50$. More stringent contrast constraints of $\Delta M=7.96$ mag at $0\arcsec.50$ were obtained with NIRC2 on the Keck II Telescope using the Br$\gamma$ filter (Furlan et al., 2017). These observations demonstrate that there are no nearby stars that might explain the $0\arcsec.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 Population-Level Trends

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, non-thermal 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 100F_{\earth}$) 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<17M_{\earth}$) 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.

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 $\sim 15\%$ precision and its radius to $\sim 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_{\mathrm{p}}=0.30\pm 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/H]}}$ prescription from Thorngren et al. (2016), this corresponds to $Z_{\mathrm{p}}/Z_{\mathrm{star}}=16.6^{+2.4}_{-2.2}$. 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.

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_{\mathrm{atm}}\leq 79\times$ solar, where “solar” refers to the Asplund et al. (2009) photospheric metal fraction of $1.04\times 10^{-3}$. 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/\mu_{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$\times$ 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_{\mathrm{p}}$ and the corresponding the H/He fraction $f_{\mathrm{H/He}}=1-Z_{\mathrm{p}}$.

As hinted at by their low bulk densities, these two planets have large H/He mass fractions: $f_{\mathrm{H/He}}=0.31\pm 0.08$ for KOI-1783.02 and $f_{\mathrm{H/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_{\earth}$. 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 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_{\mathrm{H/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 $\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 $\sim$800K, too low for Ohmic dissipation to puff up Kepler-177c to $\gtrsim$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 $\sim$0.48$R_{\oplus}$ (see equation 13 of Millholland, 2019). If the composition of Kepler-177c core is similar to that of Earth, we expect its core size to be $\sim$1.95$R_{\oplus}$ (assuming $R\propto M^{1/4}$), so that the expected total radii of the planet is only $\sim$2.43$R_{\oplus}$, far too small to explain the measured 8.73$R_{\oplus}$. Even at gas-to-core mass ratio of 10%, the expected total radii is just 3.74$R_{\oplus}$.

5.3 A Possible Formation Scenario for Kepler-177

We conclude that Kepler-177c rightfully belongs in the small sample of $\sim 15M_{\earth}$ 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) show 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 10M_{\earth}$. For cores less massive than this, the maximum gas-to-core mass ratio (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 $15M_{\earth}$ 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).

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 feasibilty 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.8M_{\earth},\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):

[TABLE]

where for simplicity we have assumed their nominal values for the $f$ factor, the nebular gas metallicity $Z$, the adiabatic gradient $\nabla_{\mathrm{ad}}$, and the temperature and mean molecular weight at the radiative-convective boundary $T_{\mathrm{rcb}}=200$ K and $\mu_{\mathrm{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 ($\sim$ 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 $\gtrsim$ 0.5, so a more careful calculation would provide an even shorter timescale. We suggest that 15$M_{\earth}$ 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 ($\rho\lesssim 0.3$ g/cm3) bulk densities may provide invaluable insights into super-puff formation processes.

6 Conclusions and Future Prospects

We presented infrared photometry for four dynamically interacting Kepler systems. With precise telescope guiding and the use of an engineered diffuser, we achieved a precision with WIRC that is comparable to or better than Spitzer for stars fainter than $J=9.5$. Most of the planets we observed have host stars that are too faint for standard Doppler-based follow-up, but their masses can be measured to a high relative precision by fitting their transit timing variations. Our new transit measurements demonstrate that a single, well-timed follow-up observation taken years after the Kepler mission’s conclusion can improve mass estimates by almost a factor of three. Perhaps unsurprisingly, we found that observing in epochs of maximally divergent transit times for differing dynamical solutions yields the largest improvements in mass estimates. The potential information gain is also larger for long-period systems with relatively few transits observed during the original Kepler mission. The systems we have studied highlight the diverse range of science cases made possible by diffuser-assisted photometry, including the confirmation of long-period planet candidates in TTV systems as well as bulk composition studies for relatively cool planets ranging in size from sub-Neptunes to gas giants.

WIRC’s demonstrated infrared photometric precision opens up multiple new opportunities for ground-based studies of transiting planets and brown dwarfs. For dynamically interacting systems bright enough for RV observations, diffuser-assisted transit observations can provide an extended TTV baseline for joint RV-TTV modeling. These kinds of studies can constrain the structures of planetary systems without reliance on stellar models (Almenara et al., 2015, 2016; Agol & Fabrycky, 2017; Weiss et al., 2017; Almenara et al., 2018; Petigura et al., 2018). For highly irradiated gas giant planets, WIRC can be used to complement existing space-based emission and transmission spectroscopy from Spitzer and the Hubble Space Telescope by observing photometric transits and secondary eclipses at wavelengths that are inaccessible to these telescopes. This extended wavelength coverage is important for reducing degeneracies in atmospheric retrievals (e.g. Benneke & Seager, 2012; Line et al., 2012, 2013, 2014). WIRC can also measure low-amplitude rotational variability in brown dwarfs at infrared wavelengths. Current ground-based infrared measurements can constrain variability at the $\sim 0.7\%$ level (Wilson et al., 2014; Radigan, 2014) in these objects; for the brighter ($J$ = 14-15) variable brown dwarfs, WIRC will be able to push these limiting amplitudes below $0.1\%$. We are only beginning to explore the parameter space made available by diffuser-assisted photometry, but the prospects for new ground-based studies of brown dwarfs and transiting planets are promising.

The authors thank the entire Palomar Observatory staff for their tireless support of our work. We additionally acknowledge Jessie Christiansen for helpful discussions on KOI-1783, B. J. Fulton for assistance with the California Kepler Survey dataset, Erik Petigura for useful comments on time-correlated noise and joint RV-TTV modeling, Nicole Wallack for discussions on light curve fitting, and Gudmundur Stefansson for conversations regarding diffuser-assisted photometry at Palomar and other observatories. Support for this program was provided by NSF Career grant 1555095 and by NASA Origins grant NNX14AD22G. This work was partially supported by funding from the Center for Exoplanets and Habitable Worlds, which is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Kepler and WIRC Light Curves

Appendix B Posterior Probability Distributions

Appendix C Dynamical Modeling Results