K2-146: Discovery of Planet c, Precise Masses from Transit Timing, and Observed Precession

TL;DR

This paper reports the discovery and detailed characterization of a two-planet system around a mid-M dwarf, including precise mass measurements, transit timing variations, and observed precession, revealing insights into their composition and orbital dynamics.

Contribution

It provides the first precise mass measurements of small exoplanets in a resonant system around an M dwarf, utilizing transit timing variations and dynamical modeling.

Findings

Both planets are sub-Neptunes with low densities.

The planets are in a 3:2 resonance with significant transit timing variations.

The planets may have retained substantial gaseous envelopes.

Abstract

K2-146 is a mid-M dwarf ($M_\star = 0.331 \pm 0.009 M_\odot$; $R_\star = 0.330 \pm 0.010 R_\odot$), observed in Campaigns 5, 16, and 18 of the K2 mission. In Campaign 5 data, a single planet was discovered with an orbital period of $2.6$~days and large transit timing variations due to an unknown perturber. Here we analyze data from Campaigns 16 and 18, detecting the transits of a second planet, c, with an orbital period of $4.0$~days, librating in a 3:2 resonance with planet b. Large, anti-correlated timing variations of both planets exist due to their resonant perturbations. The planets have a mutual inclination of $2.40^\circ\pm0.25^\circ$, which torqued planet c more closely into our line-of-sight. Planet c was grazing in Campaign 5 and thus missed in previous searches; in Campaigns 16 and 18 it is fully transiting, and its transit depth is three times larger. We improve the stellar properties using data from Gaia DR2, and using dynamical fits find that both planets are sub-Neptunes: their masses are $5.77\pm0.18$ and $7.50\pm0.23 M_{\oplus}$ and their radii are $2.04\pm0.06$ and $2.19\pm0.07$ R$_\oplus$, respectively. These mass constraints set the precision record for small exoplanets (a few gas giants have comparable relative precision). These planets lie in the photoevaporation valley when viewed in Radius-Period space, but due to the low-luminosity M-dwarf host star, they lie among the atmosphere-bearing planets when viewed in Radius-Irradiation space. This, along with their densities being 60%-80% that of Earth, suggests that they may both have retained a substantial gaseous envelope.

Tables

Table 1. Table 1: Stellar information relevant to follow-up efforts. Stellar mass, radius, and density are derived as described in Section  2.1 .

Parameter	Value	Source
RA	08:40:06.424218	(Gaia Collaboration et al., 2018)
Dec	+19:05:34.429575	(Gaia Collaboration et al., 2018)
$m_V$	$16.25\pm 0.01$	(Zacharias et al., 2013)
$m_J$	$12.183\pm 0.021$	(Cutri et al., 2003)
$m_H$	$11.605\pm 0.018$	(Cutri et al., 2003)
$m_K$	$11.370\pm 0.023$	(Cutri et al., 2003)
$T_{\mathrm{eff}}$ (K)	$3385\pm 70$	(Hirano et al., 2018)
$M_{\star }$ (M⊙)	$0.331\pm 0.009$	This work
$R_{\star }$ (R⊙)	$0.330\pm 0.010$	This work
${\rho }_{\star }$ (${\rho }_{\odot }$)	$9.3\pm 0.8$	This work
Distance (pc)	$79.31\pm 0.45$	(Gaia Collaboration et al., 2018)

Table 2. Table 2: Estimates of selected planetary, stellar, and general parameters from AIMCMC fitting of lightcurves.

\topruleparameter	planet b	planet c
$b_{C5}\left[R_{\star }\right]$	${0.12}_{-0.08}^{+0.12}$	${0.994}_{-0.008}^{+0.007}$
$b_{C16}\left[R_{\star }\right]$	${0.56}_{-0.06}^{+0.04}$	${0.900}_{-0.014}^{+0.011}$
$b_{C18}\left[R_{\star }\right]$	${0.53}_{-0.09}^{+0.06}$	${0.887}_{-0.022}^{+0.017}$
${\tau }_{C5}\left[hr\right]$	$1.27\pm 0.02$	$0.51\pm 0.07$
${\tau }_{C16}\left[hr\right]$	$1.366\pm 0.014$	$0.74\pm 0.02$
${\tau }_{C18}\left[hr\right]$	$1.352\pm 0.016$	$0.73\pm 0.03$
$r/R_{\star }$	$0.0582\pm 0.0006$	$0.0606\pm 0.0011$
$a/R_{\star }$	$16.84\pm 0.4$	$22.1\pm 0.6$
${\left(e\sin \omega \right)}_{C5}$	$-0.02\pm 0.04$	$-0.015\pm 0.18$
${\left(e\sin \omega \right)}_{C16}$	$-0.23\pm 0.04$	$0.04\pm 0.06$
${\left(e\sin \omega \right)}_{C18}$	$-0.21\pm 0.05$	$0.08\pm 0.08$
${\rho }_{\star }/{\rho }_{\odot }$	$9.1\pm 0.7$
$u_1$	$0.36\pm 0.15$
$u_2$	$0.25\pm 0.23$
$u_1+u_2$	$0.61\pm 0.10$
$\log \left({\sigma }_{C5}\right)$	$-3.443\pm 0.017$
$\log \left({\sigma }_{C16}\right)$	$-2.776\pm 0.004$
$\log \left({\sigma }_{C18}\right)$	$-2.741\pm 0.005$

Table 3. Table 3: Dynamical properties, constrained by transit shapes and times. The dynamical epoch on the below Jacobian coordinates is BKJD= 3370 3370 3370 .

\topruleparameter	planet b	planet c
$P$ [days]	$2.64460\pm 0.00006$	$4.00498\pm 0.00011$
$T_0\left[-2454833\right]$	$3371.5902\pm 0.0005$	$3371.1316\pm 0.0007$
$\sqrt{e}\cos \omega$	$0.165\pm 0.054$	$-0.177\pm 0.034$
$\sqrt{e}\sin \omega$	$-0.318\pm 0.021$	$0.200\pm 0.036$
$i$ [deg]	$88.93\pm 0.11$	$87.54\pm 0.04$
$\Omega$ [deg]	$0.0\pm 0.0$	$1.9\pm 0.3$
$m_p\left[M_{Jup}\right]$	$0.01816\pm 0.00024$	$0.02358\pm 0.00031$
$R_p/R_{\star }$	$0.0572\pm 0.0005$	$0.0612\pm 0.0009$
derived:
$e$	$0.129\pm 0.019$	$0.075\pm 0.016$
$\omega$ [deg]	$-64\pm 6$	$130\pm 5$
$m_p\left[M_E\right]$	$5.77\pm 0.18$	$7.49\pm 0.24$
$R_p\left[R_E\right]$	$2.05\pm 0.06$	$2.19\pm 0.07$

Table 4. Table 4: Planetary, stellar, and noise parameters for best-fitting lightcurve model.

pl.	$r_p/R_{\star }$	$b_{C5}/R_{\star }$	${\left(e\sin \omega \right)}_{C5}$	${\left(e\cos \omega \right)}_{C5}$	$b_{C16}/R_{\star }$	${\left(e\sin \omega \right)}_{C16}$	${\left(e\cos \omega \right)}_{C16}$	$b_{C18}/R_{\star }$	${\left(e\sin \omega \right)}_{C18}$	${\left(e\cos \omega \right)}_{C18}$
b	0.0584872	0.0829449	-0.0063295	-0.0013517	0.5971109	-0.2887000	0.0980958	0.5717864	-0.2867856	-0.2544630
c	0.0609491	1.0010763	-0.2172455	0.1146076	0.9038561	0.0052956	-0.0004055	0.9023093	0.0134860	0.0159588
	${\rho }_{\star }/{\rho }_{\odot }$	$u_1$	$u_2$	$\log \left({\sigma }_{C5}\right)$	$\log \left({\sigma }_{C18}\right)$	$\log \left({\sigma }_{C18}\right)$
	9.8060629	0.3381495	0.2524213	-3.4276080	-2.7734142	-2.7472132

Table 5. Table 5: Best-fitting dynamical model, Jacobian coordinates at dynamical epoch BJD 2458203.000 (BKJD + 3370 3370 +3370 ). Assumed M ⋆ = 0.331 ​ M ⊙ subscript 𝑀 ⋆ 0.331 subscript 𝑀 direct-product M_{\star}=0.331M_{\odot} . The fitting parameter minimized is ( 2 ​ ln ⁡ Δ ​ ℒ ) = 163.3337 2 Δ ℒ 163.3337 (2\ln\Delta{\cal L})=163.3337 . Planet masses given in units of M Jup = 9.545 × 10 − 4 ​ M ⊙ subscript 𝑀 Jup 9.545 superscript 10 4 subscript 𝑀 direct-product M_{\rm Jup}=9.545\times 10^{-4}M_{\odot} .

planet	$P$ (days)	$T_0$ (BKJD)	$\sqrt{e}\cos \omega$	$\sqrt{e}\sin \omega$	$i$ (deg)	$\Omega$ (deg)	$M_p$ ($M_{\mathrm{Jup}}$)
b	2.64461662195	3371.59032337250	0.14179457798	-0.31149532326	88.94655839504	0 (definition)	0.01825083297
c	4.00491983123	3371.13139767470	-0.19543722063	0.21749371377	87.51200896614	2.01729810288	0.02371931622

Table 6. Table 6: A preview of samples described in Section  2.4 . The full table is available as a machine-readable table.

${\rho }_{\star }\left({\rho }_{\odot }\right)$	$P_b$ (days)	$T_{ob}$ (BKJD)	$\sqrt{e_b}\cos {\omega }_b$	$\sqrt{e_b}\sin {\omega }_b$	$i_b$ (deg)	${\Omega }_b$ (deg)	$M_b\left(M_{Jup}\right)$	$R_b\left(R_{\star }\right)$
10.2549796	2.64459256	3371.59078	0.13014153	-0.3336921	88.7788353	0.0	0.01821541	0.0581626
8.63854893	2.64469711	3371.59074	0.14716385	-0.3012654	88.9552587	0.0	0.01833648	0.0576053
9.50392865	2.64457815	3371.59015	0.16020293	-0.3429077	88.8031510	0.0	0.01813466	0.0575424
…	…	…	…	…	…	…	…	…
	$P_c$ (days)	$T_{oc}$ (BKJD)	$\sqrt{e_c}\cos {\omega }_c$	$\sqrt{e_c}\sin {\omega }_c$	$i_c$ (deg)	${\Omega }_c$ (deg)	$M_c\left(M_{Jup}\right)$	$R_c\left(R_{\star }\right)$
	4.00502692	3371.13169	-0.1995415	0.19396097	87.5118793	1.48106165	0.02360255	0.0599834
	4.00483110	3371.13112	-0.2054754	0.22162647	87.4948920	2.44051956	0.02378850	0.061674
	4.00498325	3371.13186	-0.1783291	0.17854464	87.5673562	1.99859017	0.02360026	0.0608582
	…	…	…	…	…	…	…	…

Table 7. Table 7: Expected future transit times, for planet b.

\toprulePlanet	Transit Number	Transit time	Uncertainty
		(BKJD)	(minutes)
b	38	3472.7760	0.5
	39	3475.4451	0.5
	40	3478.1139	0.5
	41	3480.7802	0.5
	…	…	…
	568	4880.9795	16.8
	569	4883.6502	16.4
	570	4886.3208	16.2
	571	4888.9916	15.9
	572	4891.6610	15.4
	573	4894.3304	14.9
	574	4897.0000	14.4
	575	4899.6673	13.5
	576	4902.3347	12.8
	577	4905.0023	12.1
	578	4907.6667	11.0
	579	4910.3313	10.2
	580	4912.9961	9.3
	581	4915.6567	8.4
	582	4918.3177	7.7
	583	4920.9788	7.2
	584	4923.6353	7.1
	585	4926.2921	7.4
	586	4928.9490	8.1
	587	4931.6015	8.9
	588	4934.2542	9.9
	589	4936.9070	11.1
	590	4939.5559	12.0
	591	4942.2052	13.1
	592	4944.8544	14.3
	593	4947.5008	15.0
	594	4950.1473	15.9
	595	4952.7937	16.8
	596	4955.4384	17.4
	597	4958.0829	18.0
	598	4960.7276	18.6
	599	4963.3711	18.9
	600	4966.0145	19.3
	601	4968.6580	19.7
	602	4971.3011	19.8
	603	4973.9439	20.0
	604	4976.5871	20.1
	605	4979.2302	19.9
	606	4981.8730	19.9
	607	4984.5165	19.7
	608	4987.1603	19.4
	609	4989.8035	19.1
	610	4992.4481	18.6

Table 8. Table 8: Expected future transit times, for planet c.

\toprulePlanet	Transit Number	Transit time	Uncertainty
		(BKJD)	(minutes)
c	25	3470.5692	0.6
	26	3474.5371	0.6
	27	3478.5058	0.6
	28	3482.4772	0.6
	29	3486.4493	0.6
	…	…	…
	374	4861.7518	16.8
	375	4865.7170	17.2
	376	4869.6814	17.5
	377	4873.6461	17.6
	378	4877.6103	17.5
	379	4881.5752	17.4
	380	4885.5402	17.0
	381	4889.5062	16.5
	382	4893.4730	15.7
	383	4897.4409	14.9
	384	4901.4109	13.7
	385	4905.3817	12.5
	386	4909.3560	11.0
	387	4913.3311	9.5
	388	4917.3111	8.0
	389	4921.2916	7.0
	390	4925.2779	6.7
	391	4929.2648	7.3
	392	4933.2571	8.5
	393	4937.2503	10.2
	394	4941.2479	11.7
	395	4945.2465	13.4
	396	4949.2481	14.6
	397	4953.2508	15.9
	398	4957.2552	16.8
	399	4961.2606	17.7
	400	4965.2668	18.3
	401	4969.2738	18.8
	402	4973.2808	19.0
	403	4977.2884	19.1
	404	4981.2954	19.0
	405	4985.3023	18.6

Equations

$$\frac{b}{R_{\star }}=\frac{a\cos i}{R_{\star }}\frac{1-e^2}{1+e\sin \omega }.$$

$$\frac{\tau }{T}=\frac{\sqrt{(r_p+R_{\star }{)}^2-b^2}}{\pi a}\frac{\sqrt{1-e^2}}{1+e\sin \omega }.$$

$$\Delta \mathcal{L}=\exp ({\chi }^2/2)+\exp (-\text{text@frac}12{\vec{x}}^t{\mathbf{A}}^{-1}\vec{x})$$

$$f=\frac{\pi }{2}-\omega$$

$$f=\frac{\pi }{2}-\omega -\frac{e\cos \omega {\cos }^2(i)}{1+e\sin \omega }$$

$$\begin{array}{rl}\Delta t& =\Delta f\frac{(1-e^2{)}^{3/2}}{(1+e\sin \omega {)}^2}\frac{T}{2\pi }\\ & =\frac{-e\cos \omega {\cos }^2(i)(1-e^2{)}^{3/2}}{(1+e\sin \omega {)}^3}\frac{T}{2\pi }\end{array}$$

Full text

Planets and Satellites: Detection – Planets and Satellites: Dynamical Evolution and Stability – Planets and Satellites: Individual (EPIC 211924657, K2-146)

††software: numpy (Van Der Walt et al., 2011), matplotlib (Hunter et al., 2007), lightkurve (Barentsen et al., 2019), isochrones (Morton, 2015), astropy (Astropy Collaboration et al., 2018), scipy (Jones et al., 2001–), kadenza (Barentsen & Cardoso, 2018), everest (Luger et al., 2018b), emcee (Foreman-Mackey et al., 2013), corner (Foreman-Mackey, 2016), batman (Kreidberg, 2015)

1 Introduction

The K2 mission (Howell et al., 2014) observed 20 fields across the ecliptic plane. Because of the limitations of the spacecraft, with the loss of two reaction wheels, the telescope could only stably point in limited directions. With much of the ecliptic plane already observed by the telescope, later campaign fields overlapped with previous fields and some stars were observed in two or three campaigns.

Campaign 5 was observed as a part of the K2 mission from 2015 Apr 27 to 2015 Jul 10. Campaign 16 revisited a fraction of this field from 2017 Dec 07 to 2018 Feb 25, while Campaign 18 again revisited this field from 2018 May 12 to 2018 Jul 02. These observations combine to provide 205 days of observations over the “K2 Legacy Field” spread over 3.2 years. By combining these data for stars observed by multiple campaigns, planetary systems can be characterized in higher detail than for each individual campaign alone (e.g. Chakraborty et al., 2018). Moreover, when planets orbiting the same star are dynamically interacting, the long observing baseline provides an opportunity to observe and understand transit timing variations (TTVs) to measure planet masses and orbital parameters, as was commonly done with data from the original Kepler mission (e.g. Huber et al., 2013; Hadden & Lithwick, 2014; Jontof-Hutter et al., 2016).

1.1 K2-146 system

One of the stars observed in Campaigns 5, 16, and 18 was EPIC 211924657, also named K2-146 after the validation of its first planet.111This star was proposed as a target in many programs. For Campaign 5: GO5020, GO5097, GO5054, GO5011, and GO5006. For Campaign 16: GO16005, GO16011, GO16052, GO16077, GO16101, and GO16015, the final two of which requested short cadence observations. For Campaign 18: GO18068, GO18027, GO18048, and GO18063, of which all but the first requested short cadence observations. When Campaign 5 data was available, Pope et al. (2016) quickly identified a planet candidate with a 2.6 day period; the planet was also detected by the Petigura et al. (2018) search, and the host star and planet better characterized with follow-up spectroscopy by Dressing et al. (2017) — who also noted that the planet showed TTVs. Hirano et al. (2018) performed a detailed analysis of the system, validating the planet by combining spectroscopy and adaptive optics imaging with the K2 data. These authors also noticed the presence of TTVs, finding a rapid change in the observed orbital period of $\sim 10$ minutes, at time BKJD222BKJD = BJD - 2454833.$=2360$. Livingston et al. (2018) similarly analyzed and validated planet b while noting the existence of TTVs. Each of these teams only reported the presence of a single transiting planet in Campaign 5, and hence were not able to interpret the cause of these TTVs.

1.2 Precessing Planets

It was recognized early in the history of transiting exoplanets that a planet with a large impact parameter, nearly grazing its star, is exquisitely sensitive to additional planets via the out-of-plane torque exerted between the orbits (Miralda-Escudé, 2002). A changing inclination would generate large transit duration/depth variations (TDVs). In a first example, tentative upper limits on the transit depth of GJ436b seemed to conflict with later measurements of its transits, indicating that its inclination was changing on yearly timescales Ribas et al. (2008). Later, a series of transits of GJ436b were shown to have a constant shape (Ballard et al., 2010), providing upper limits on inclination change. Thus a limit was placed on the product of the mass, $m_{c}$, and rotation around the line-of-sight, $\Delta\Omega$, of secular perturbers of various periods. Our detection, to be described, allows us to constrain planetary masses by TTVs and then determine $\Delta\Omega$ via interpreting the depth change, much as envisioned by Ribas et al. (2008).

Since then, several more robust examples of orbital precession in exoplanets have been discovered.

First, Kepler-13 was shown to have both a non-symmetric transit shape and a changing transit duration (Szabó et al., 2012). Both are attributed to a misaligned orbit around a rapidly spinning host star: the former is due to the temperature decrease from pole to equator, and the latter is due to the torque from the stellar oblateness.

Second, several circumbinary planets have orbits misaligned from the orbital plane of the host binary; the torque from the binary causes the planets’ orbits to precess on decade-long timescales (Welsh et al., 2012; Kostov et al., 2014).

Finally, three planetary systems have shown planet-planet interactions that precess one another’s orbits, Kepler-117 (Almenara et al., 2015), Kepler-9 (Freudenthal et al., 2018; Borsato et al., 2019) and Kepler-119 (Mills & Fabrycky, 2017). The latter requires a large mutual inclination of $24^{+11}_{-8}$ degrees to fit the TTVs and TDVs without invoking additional planets.

Another effect that can be probed by TDVs is in-plane precession of eccentricity. The resonant term of periastron torque (on the libration/TTV period) has been inferred by Nesvorný et al. (2013), but so far no claim of a secular term has been made, to our knowledge. As described by Pál & Kocsis (2008), (a) at low impact parameter, precession of the periastron affects TDVs by changing the velocity of the planet at transit: transit near periastron yields a shorter than average duration; (b) at high impact parameters, precession changes the distance between the planet and the star at transit, which affects the transit chord length: transit near periastron yields a longer transit chord and a longer than average duration. For small orbital eccentricities, the impact parameter separating these two regimes is $\sqrt{1/2}$, and we shall show that the two-planet system K2-146 furnishes one example in each of these regimes. These sources of TDVs, and others, are described in Agol & Fabrycky (2018).

1.3 Verification of TTV perturbers

TTVs were proposed as a method of detecting companions of transiting exoplanets (Agol et al., 2005; Holman & Murray, 2005). A clear precedent for our discovery has been made in Kepler-36 (Carter et al., 2012), where a planet with a 16-day orbital period was seen to exhibit strong and sharp (non-sinuosidal) TTVs. It took about two years of data to recognize that its perturber was also transiting, after which it was immediately clear that the new, 13-day planet could explain the TTV pattern.

Recently there have been a few other examples of systems where a single planet has TTVs, and the perturber was predicted and later verified by radial velocity. Barros et al. (2014) were the first to achieve this, finding the perturbing companion predicted by Nesvorný et al. (2013). A planetary perturber was inferred in the Kepler-19 system by Ballard et al. (2011), albeit with incomplete information, and a radial-velocity planet with the right properties was later detected by Malavolta et al. (2017). In the Kepler-419 system, a highly eccentric planet’s TTVs predicted an additional planet (Dawson et al., 2012) which was then confirmed and further characterized by radial velocity measurements (Dawson et al., 2014). A perturbing planet’s mass and period was determined via highly precise TTV data by Nesvorný et al. (2012); the second planet has not yet been detected by other means. In none of these systems is the perturbing planet known to transit, but with enough mutual inclination, it might be torqued by the currently-transiting planet into a transiting orbit of its own.

In Section 2 we discuss how we update the understanding of K2-146 using Gaia data, how we treat the light curve, and how we model the transit timing. Section 3 gives the numerical results. In Section 4 we discuss the dynamical evolution of the system. Finally, in Section 5 we conclude with inferences we may draw and further steps worth taking in the study of this exciting system.

2 Data Analysis

2.1 Stellar Parameters enabled by Gaia

Hirano et al. (2018) obtained a spectrum of K2-146 with the High Dispersion Spectrograph on the Subaru 8.2 m telescope (Noguchi et al., 2002). Those authors inferred a stellar effective temperature of $3385\pm 70$ K, a metallicity [Fe/H] of $-0.02\pm 0.12$, and a surface gravity log(g) of $4.906\pm 0.041$. From these stellar parameters, the authors infer a stellar radius of $0.350\pm 0.035$ R⊙. This star was also observed by Rodríguez Martínez et al. (2018), who use near-IR spectroscopy from Palomar/TripleSpec to infer a temperature of $3766\pm 195$ K, a log luminosity of $-1.91\pm 0.16$ L⊙, and a radius of $0.68\pm 0.16$ R⊙.

We combine the spectroscopic stellar parameters with data from Gaia DR2 (Gaia Collaboration et al., 2016, 2018) in order to improve the precision on the inferred radius. The published DR2 parallax for this star is $12.5821\pm 0.0750$ milliarcsec. We apply a correction of $+0.029$ mas to this value to account for the zeropoint offset inherent in the DR2 results (Lindegren et al., 2018). We use this parallax, the inferred stellar parameters from Hirano et al. (2018), and JHK photometry from the 2-Micron All Sky Survey (2MASS, Cutri et al., 2003), as listed in Table 1, as inputs to the isochrones package, a Python interface to fit derived stellar parameters, photometry, and parallax information simultaneously to stellar models (Morton, 2015).

Fitting these data to the MESA Isochrones and Stellar Tracks (MIST) models (Choi et al., 2016), we find an inferred radius of $0.323\pm 0.008$ R⊙ and mass of $0.323\pm 0.011$ M⊙. Alternatively, fitting to the Dartmouth Stellar Evolution Database (Dotter et al., 2008) result in a radius of $0.330\pm 0.005$ R⊙ and mass of $0.331\pm 0.007$ M⊙. We choose to use the results of the fit to the Dartmouth models, which provide a higher maximum log-likelihood fit to our input data $\Delta\mathcal{\log L}=2.92$, although we note that both sets of stellar models may suffer from systematic uncertainties at the few percent level in their inference of stellar radii. Mann et al. (2015) fit an empirical relation between the $K_{s}$ absolute magnitude, metallicity, and radius of M dwarfs. Using this relation, we predict a radius of 0.326 R⊙, in line with our model predictions. These authors find a $2.7\%$ scatter in M dwarf radii around this relation, and a $1.8\%$ scatter in M dwarf masses. We consider these values to be extra systematic uncertainties above the model predictions and add them in quadrature to our derived parameters. From this, we find a radius for K2-146 of $0.330\pm 0.010$ R⊙ and a mass of $0.331\pm 0.009$ M⊙, which are the values we employ through the remainder of this work. We also calculate the density of K2-146 to be $9.3\pm 0.8$ $\rho_{\odot}$, where the uncertainty mostly arises from the aforementioned scatter in M dwarf radii and masses.

Recently, Parsons et al. (2018) performed a uniform analysis of low-mass stars, noting that slowly rotating M dwarfs are consistent with models, albeit with a large amount of scatter. As there is no obvious rotation period signal in the light curve, nor are any flares observed during the campaigns, the star is consistent with being a slowly rotating M dwarf, so we apply the results of the Dartmouth models without including any corrections to its output.

Using the less precise derived stellar parameters from Rodríguez Martínez et al. (2018) and repeating this analysis, we find an increase in the derived mass and radius of less than 1%, smaller than our quoted uncertainties. In this case, the Gaia DR2 parallax combined with the broadband photometry enable a precise measurement of the stellar mass and radius, overwhelming the spectroscopically-derived temperature and log luminosity, which isochrones accepts as priors. The stellar parameters and their uncertainties used in this work are listed in Table 1.

2.2 Location of transits using Kepler

2.2.1 Data acquisition and search for periodic signals

When first examining the system, our goal is to determine the number of significant transit signals and obtain initial estimates of transit times. For this precursory work, we use long cadence data from campaigns 5, 16, and 18 .

We obtain C5 data from the EVEREST database of detrended light curves at the Mikulski Archive for Space Telescopes (MAST) (Luger et al., 2018a). The C16 and C18 light curves we download from MAST unprocessed, as those campaigns have not yet been fully processed within the EVEREST pipeline. We use the standalone.py module within the EVEREST Python package to detrend these light curves locally, selecting apertures that encompass the star throughout their respective campaigns and minimize the CDPP (Combined Differential Photometric Precision, a measure of the light curve’s noisiness described in Christiansen et al. (2012)). The chosen apertures contained 16 pixels for C16 and 14 pixels for C18, and are not altered during their campaigns. We further detrend and normalize the light curves by dividing each by a running median spanning approximately 2 days. At this point, we note that no data points fall more than 5ppt below unity, while some outliers remain as high as 80 ppt above unity. Thus, we mask out data points more than 2ppt, corresponding to 4-5 $\sigma$.

For each campaign, we use the box least squares (BLS) algorithm to search for periodic signals. Initially, we consider 10 transit durations from 1.2 to 4.8 hours and a range of periods from 0.4 d to 30 d. However, when it becomes clear that no significant signals are present at the low or high ends of this period range, we restrict our search to between 1 d and 10 d to save computation time. This process reveals a significant signal near 2.6d and a questionable signal near 4.0d in C5, and significant signals at both stated periods in C16 and C18. While the 4.0d signal in C5 is far from significant, we examine it further because it coincides with significant signals present in C16 and C18. We use these periods to fold the data separately for each campaign, resulting in 6 folded light curves.

2.2.2 Determination of transit midpoints

We perform the process that follows for each of these six curves individually. We use the batman Python package (Kreidberg, 2015), which analytically calculates transit models using a method detailed in (Mandel & Agol, 2002), to plot over the top of the folded light curve. To develop an initial estimate for the transit parameters, we adjust five of the model’s parameters manually until model and folded lightcurve appear to be in reasonable agreement. The parameters we adjust here are: the ratio between the orbit’s semimajor axis and the star’s radius ($a/R_{\star}$), the ratio between the planet’s radius and the star’s radius ($r_{p}/R_{\star}$), the inclination angle of the orbit ($i$), and two limb-darkening parameters (LDPs) for a quadratic limb-darkening model ($u_{1}$ and $u_{2}$). At this point, the orbit is taken to be circular, although this constraint is removed later in our analysis. Using the period estimate provided by the BLS search and phase information from the folded curve, we select a window of data points around each transit. For each window, we compute 1000 batman models with midpoint times spaced evenly throughout the window. For each midpoint time, we compare the model and data within the window and calculate the resultant likelihood. We take our prior distribution for each transit’s midpoint time to be flat, so our posterior distribution is directly proportional to the calculated likelihoods. We take the midpoint time with maximum likelihood to be our point estimate of a given transit’s true midpoint time. We assume the likelihoods to be Gaussian, and thus divide their full widths at half maximum by 2.355 to estimate the uncertainties in the transit times. In cases of multimodal posteriors, we use the widest pair of points at half maximum in order to avoid underestimating transit time uncertainties. We repeat this process—using the new transit times to fold the data, allowing for better guesses of batman parameters, leading to more accurate transit times—and find that after three iterations we cannot noticeably improve the fit simply adjusting the transit parameters manually.

To improve the fit from here, we perform a more careful analysis while switching to short cadence data for C16 and C18 for greater temporal resolution. Because the short cadence data is prohibitively large for detrending in full, we split it into pieces roughly 2.7d long that consist of a similar number of data points as a full campaign of long cadence data. We detrend these pieces individually with EVEREST while masking data around the transit times obtained previously. Because further detrending is heavily impacted by outliers, we mask out points more than about 7 ppt above the highest point or below the lowest point of a running median spanning a campaign with a width of 1.5 hours. 7 ppt is about $4\sigma$ for the short cadence data and large enough that we don’t risk eliminating transits. We calculate the best-fit 4th-order polynomial for each piece while masking out the transits, and use these polynomials to detrend and normalize the fluxes. We also improve our C5 lightcurve by performing the same piecewise polynomial detrending on the masked post-EVEREST C5 long cadence data. We calculate the residuals between the folded, normalized light curves and the transit models and then repeatedly mask out points with residuals more than 3 standard deviations from zero until no 3-$\sigma$ outliers remain.

Next, we use scipy to minimize $\chi^{2}$ between data and transit model by optimizing twelve parameters: $a/R_{\star}$ for both planets, $r_{p}/R_{\star}$ for both planets, inclination angle for three campaigns each for both planets, and two LDPs. Finally, we improve our measurements of the transit midtimes using the updated light curves and transit parameters, following the same procedure described above.

The transit times we obtain from this analysis are not constrained to physically feasible orbits. Thus, we perform a preliminary TTV analysis similar to the one described in Section 2.4 to obtain a set of physically sensible transit times. For this preliminary analysis, we fit to the photometrically-derived transit times and rate of change of the outer planet’s impact parameter. The times obtained in this way are more precise because they align with smooth ephemerides, while the times obtained by individually fitting transits to a transit shape have freedom to align with photometric noise.

2.3 Lightcurve Fitting

We use an affine-invariant Markov chain Monte Carlo algorithm (AIMCMC, Goodman & Weare (2010)) via emcee (Foreman-Mackey et al., 2013) to obtain posteriors on 26 parameters that influence the shapes of the planets’ transits. We do this by computing batman models for each planet and campaign and comparing them with lightcurves folded on the dynamically-constrained transit times. The use of an MCMC method allows us to explore 26 degrees of freedom and form quantitative posteriors near the best-fit in a computationally feasible manner.

The stellar fit parameters are $\rho_{\star}/\rho_{\odot}$ and two quadratic LDPs (3 parameters). Rather than fitting eccentricities $e$ and arguments of periastron $\omega$ themselves, we fit $e\sin{\omega}$, and $e\cos{\omega}$. We fit impact parameter $b/R_{\star}$, $e\sin{\omega}$, and $e\cos{\omega}$ separately for both planets and all three campaigns (18 parameters). Each planet’s radius is fit in the form of $r/R_{\star}$ (2 parameters). Finally, the relative flux uncertainty $\sigma$ is fit separately in each campaign to ensure we neither overestimate nor underestimate them (3 parameters). Because the planets’ periods are known precisely, we fit $\rho_{\star}/\rho_{\odot}$ and calculate both $a/R_{\star}$ values from it. $b/R_{\star}$ can be constrained by photometry because of its close relation to transit depth. To generate a transit model, batman requires $i$, $e$, and $\omega$. While we don’t fit for those parameters explicitly, we can calculate $e$ and $\omega$ from $e\cos{\omega}$ and $e\sin{\omega}$; we calculate $i$ using the relationship

[TABLE]

where $\omega=90^{\circ}$ when the transit occurs at periastron. We choose to fit $e\sin{\omega}$ and $e\cos{\omega}$ instead of $e$ and $\omega$ because the transit duration is primarily affected by $e\sin{\omega}$. This can be seen from the following expression for the transit duration as a fraction of orbital period:

[TABLE]

(In this paper, transit duration always refers to the time between beginning of ingress and end of egress, since other definitions can become undefined for grazing orbits.) The second eccentricity component, $e\cos{\omega}$, plays only a minor role, and is largely unconstrained photometrically when an artificial timing offset is corrected for (see the Appendix for a description of this effect).

We utilize a Gaussian prior on $\rho_{\star}/\rho_{\odot}$ based on the results of Section 2.1. We find that using priors on the LDPs based on the results of Claret & Bloemen (2011) does not significantly alter the best-fit values for the other parameters, so we use flat priors for limb-darkening as well (constraining to physical situations, e.g. negative limb brightnesses are disallowed). We place flat priors on $b/R_{\star}$ and restrict them to be less than 1.1 to ensure the planets transit. We restrict $e\sin{\omega}$ and $e\cos{\omega}$ to be between -0.3 and 0.3, as we don’t expect larger eccentricities to be dynamically feasible. Further, we impose a log-likelihood penalty for each planet and campaign equal to the logarithm of eccentricity, which effectively converts our priors from flat in $e\sin{\omega}$ and $e\cos{\omega}$ to priors flat in eccentricity and $\omega$. Lastly, we use priors flat in $r/R_{\star}$ and $\log_{10}{\sigma}.$

For the AIMCMC analysis, we create 400 walkers with starting parameters randomly and independently chosen from uniform distributions spanning the values that seem most likely. These walkers are evolved for 10500 steps of burn-in, allowing them to explore the region of phase space that yields reasonable fits. We verify convergence of the walkers to the posterior distribution for all 26 parameters using the test of Geweke (1992). We evolve them for a further 10000 steps to collect samples. For each parameter, the autocorrelation function of each walker’s value is averaged together to obtain a total autocorrelation function for that parameter. The number of steps required for it to decrease to $1/e$ is taken as the autocorrelation length. The autocorrelation lengths of the 26 parameters are found to be between 650 and 900, giving us an effective minimum of 11 independent samples per walker, for a total of 4400 samples.. We present corner plots of these posteriors in the Appendix (figs. 17 and 18), and highlight selected results in Section 3.

2.4 Transit Timing Analysis

We perform N-body simulations to fit the transit midtime data and impact parameter shifts. These integrations use the GNU Scientific Library’s implementation of the Prince-Dormand method (gsl_odeiv2_step_rk8pd) to integrate Newton’s equations of motion for 3 bodies. The stellar mass is taken at the best-fit value of $0.331M_{\odot}$ and not varied, but for the purpose of computing durations the radius $R_{\star}$ is taken as a fitting parameter and is constrained by the value of $\rho_{\star}$ derived from the lightcurve shape. The orbital parameters ($P$, $T_{0}$, $\sqrt{e}\cos\omega$, $\sqrt{e}\sin\omega$, $i$, and $\Omega$) — a priori taken as uniformly distributed — are defined at a time in the center of the short cadence data, namely $BKJD=3370$. Each planets’ mass ($m_{p}$) and radius ($R_{p}/R_{\star}$) are also fitting parameters. The planets’ orbits are evolved backwards in time to find mid-transit times during C16 and C5 and forwards to find mid-times in C18. We hone in on mid-times by a Newton-Raphson technique described by Fabrycky (2010). When the sky-projection of the separation vector and the velocity vector are perpendicular, we record them to yield the impact parameter (in $AU$) and the transit durations; these are averaged over campaign to compare with the lightcurve fits.

These theory values are compared with the transit timing data via a $\chi^{2}$ statistic. An additional likelihood value is taken from the correlation matrix of the shape parameters. The full likelihood constraining the solutions is:

[TABLE]

Here $\vec{x}=\vec{x}_{\rm model}-\vec{x}_{\rm data}$ is the difference between the parameters predicted by a model to the parameters that best-fit the data; these are given in Table 2. The covariance matrix ${\bf A}$ we used to constrain the model was computed from 363638 samples of the shape model. Due to the low impact parameter allowed for planet b in campaign 5, with a probability indifferent to different models near $b=0$, half of the values of that parameter were flipped negative, so that its best value is taken to be [math] with Gaussian profile that better matches the actual samples.

The median values and uncertainties on parameters are determined through Differential-Evolution Markov-Chain Monte Carlo (DEMCMC, Ter Braak 2006). Forty walkers are initialized at the best-fit that was found by a Levenberg-Marquardt algorithm, with separations in the high-dimensional space by small amounts. At the heart of the algorithm is choosing what step to try for each parameter, to produce the next generation of parameter values. To make that choice for each walker, two other walkers are chosen randomly, and their separation in each parameter differenced. That difference vector is multiplied by a number $\gamma r$, that makes the step more conservative. The value $r$ is chosen randomly each draw. The value $\gamma$ is common for all walkers, but it is a function of generation number. It starts at $2.38/\sqrt{2N}=0.4$, where $N=17$ is the number of dimensions, but it is multiplied each generation by either $0.9$ or $1.1$ to make the steps more or less conservative, based on whether fewer or more than $23\%$ of the past generation’s proposed steps were accepted. In this way, we observed $\gamma$ to settle to, then vary around, $\sim 0.1$. With each walker having been assigned a trial parameter set, the model is computed and then the step is either taken if the likelihood either increased or it decreased with the ratio of the new likelihood to the old likelihood that was greater than a draw from a uniformly-distributed random number. For jumps that are not accepted, that walker gives another copy of the current state in the chain.

Given the small differences in initialization of parameters, the first steps are very conservative and most proposed jumps are accepted, leading to the separations of the walkers growing exponentially until they span the region allowed by the data, including its covariances. The values of all parameters are plotted versus generation, with initial ones discarded as a burn-in, and the latter ones validated as fair samples via measuring the autocorrelation. In different parameters, the autocorrelation length was between 450 and 1650 steps, and we have completed 98000 steps after a burn-in of 2000 steps (for each of the 40 walkers), giving at least 2376 effectively independent samples of the posterior. Using the relationship between effective sample size and tolerance level discussed by Vats et al. (2019), we expect this procedure to sufficiently sample the parameters’ posteriors at tolerance levels between 0.04 and 0.08, corresponding to credibility intervals between 0.96 and 0.92.

3 Results

3.1 Lightcurve results

Light curves folded on the transit midpoints obtained via orbit modeling are plotted alongside best-fit transit models in Figures 2 and 3. The two figures correspond to the two planets, and the subplots to the three campaigns of data. In both cases, the colored lines represent unbinned data, while the black lines represent data binned in 30-minute intervals. Since the duration of a short cadence is much shorter than the duration of a transit, the unbinned models are suitable for comparison to short cadence data. The binned models are plotted only for C5, where they provide a more realistic comparison to the long cadence data than the unbinned models do.

We list our photometry-based estimates of the transit parameters in Table 2. Note that these are not our final parameter estimates; those are updated with the dynamical simulations and instead presented in Section 2.4. Our estimates of the LDPs are consistent with the values we derive from Claret & Bloemen (2011) of $u_{1}=0.36\pm 0.08$ and $u_{2}=0.35\pm 0.06$. We find that the sum of the LDPs is more well-determined than either of their values individually. This is likely because their sum directly determines the limb darkening at an impact parameter of $R_{\star}$, and thus is constrained by the depths of planet c’s transits. Our fitting procedure leads to $\sigma$ values 15-30% lower than the median relative flux errors listed in the EVEREST-generated lightcurves. This corresponds to a difference in $\log(\sigma)$ of 0.07-0.14, a highly significant difference.

Our photometric posteriors for impact parameters are shown in more detail in Figures 4 and 5. These distributions provide significant evidence that the impact parameter of both planets have changed over time, and this is reflected in the results of the dynamical simulations of Section 2.4. The use of these impact parameter distributions as priors in our fitting of the dynamics constrains the system to configurations that yield large ($\sim 0.07-0.1R_{\star}$) changes in impact parameters over the 3-year baseline between C5 and C18. Impact parameters near $R_{\star}$ are much easier to constrain via depth and duration of transit, which is why the impact parameter of planet c is far less uncertain than that of planet b. We note that for planet c to not be grazing during C5, our estimated impact parameter would need to be incorrect by nearly 7$\sigma$.

Planet b’s impact parameter is constrained to grow because its transits are deeper during C5 than during C16 and C18. Notably, its transits are briefer during C5 despite its lower impact parameter causing it to cross a longer chord of the star. This constrained $e\sin{\omega}$ to be negative for planet b in C16 and C18 because transit durations are inflated near apastron for eccentric orbits, as can be seen from equation 2.

Figure 6 shows observed and modeled long cadence light curves for planet c. The models use the same transit parameters as Figure 3, but account for the blurring effect of 0.02-day long cadences. This blurring serves to reduce the depth of the observed signal by approximately 30% in C5 relative to what would be observed in short cadence data, as can be seen by comparing the solid and dashed green curves. Even when comparing to C16/18 long cadence data, planet c’s transits are 3 times shallower in C5 due to the planet’s area partially missing the star during the grazing transit and also to the stronger effect of limb darkening at the edge of the star. These depth-reducing effects contributed significantly to the original nondetection of planet c in the C5 data.

Planet c’s detection is also made more difficult by its large-amplitude TTVs and short transit duration. For planets which display a sinusoidally-varying ephemeris with an amplitude which exceeds the transit duration, as is the case for K2-146 b/c, the BLS algorithm can degrade in sensitivity due to the transit-timing variations causing a “smearing” of the transit signal (García-Melendo & López-Morales, 2011). In practice this loss of sensitivity is primarily an issue for very shallow transits for which each transit is at or below the noise level, such as planet c during C5. For deeper transits, like planet b, the BLS algorithm can easily identify the linear portions of the ephemeris, with a reduction in the signal-to-noise arising from areas outside of the linear region which only contribute noise.

Finally, planet c’s transit duration in C5 of just 30 minutes is comparable to a single long cadence. Most planet searches use grids of transit durations with a minimum of 1.5-2 hours. These searches will therefore have reduced sensitivity to any planets with transit durations much shorter than that minimum search duration. K2’s frequent thruster fires caused many one cadence outliers, which likely also made other groups less willing to search for bona fide candidates of such short durations because of the increased rate of false positives. Thus, the nondetection of planet c in all previous searches of C5 can be attributed to its shallow depth, duration shorter than the minimum used in the searches, and TTVs smearing out any remaining signal.

Despite these challenges, a dedicated search would have been able to detect planet c with the appropriate search parameters, most notably using grids of shorter transit durations. As discussed in Section 2.2.1, our own BLS search uncovered a tentative signal of planet c in C5 alone using a minimum duration of 1.2 hours — more than twice the transit’s actual duration. Kruse et al. (submitted ) searched the star with an updated version of QATS designed specifically to detect small planets with TTVs. While they also found planet b, they missed planet c: again because their minimum search duration of 2 hours was too long. Re-searching the star with their new version of QATS, but with transit durations down to half an hour, recovers planet c in C5 alone, as well as capturing its anti-correlated TTVs to planet b.

3.2 TTV results

Here we report results of the transit timing analysis by an N-body code driven by DEMCMC, as described in Section 2.4, in table 3 and Figures 7 and 8. In the electronic journal article, readers may download chains that are thinned to every 200 steps, and a brief sample of these is shown in the Appendix in Table 6.

Correlation plots of each parameter versus the others are given in the Appendix, in Figure 19.

We fixed $M_{\star}$ in the dynamical analysis, described in Section 2.4. In order to fit the shape constraint on $\rho_{\star}$, we solved for the parameter $R_{\star}$, as well as $R_{p}/R_{\star}$ for each planet. We obtain $\sim 1.0\%$ and $\sim 1.5\%$ precision on $R_{p}/R_{\star}$ values of planets b and c. Lightcurves alone can only yield densities, as the equations of motion connect distances with the $1/3$ power of mass. So, after the transit timing fitting, we rejection-sample the joint distribution of $M_{\star}$ and $R_{\star}$ obtained in Section 2.1 based on the updated posterior for $\rho_{\star}$. Finally, we derive $M_{p}$ and $R_{p}$ values from the ratios relative to the star multiplied by the new stellar mass and radius, and propagate the uncertainty from both sources of uncertainty. These planetary masses and radii are given as the last two entries to table 3.

4 Dynamical Evolution

In this section we explore the dynamical behavior of the system.

The aspect of the dynamics which is linked most closely to the large TTVs is the resonant libration. The planets have resonant arguments composed of their mean longitudes and longitudes of periastron. These oscillate around a mean value, which itself secularly changes (see Figure 9). The long axes of the two orbits librate around an antiparallel configuration.

The dynamical results require that the two planets are librating in the 3:2 resonance.

Figure 10 shows the eccentricity components as a function of campaign for each planet. The increase in transit duration of the outer planet is partially due to out-of-plane torque changing the inclination of the planet to the sky plane, and partially due to in-plane precession causing the periastron to swing towards the observer. This latter effect shortens the distance between the planet and the star by about 5%, and the sky-projected projected component of this separation therefore decreases by 5%. Thus, the impact parameter from 1.0 to 0.95 by that effect alone, which dramatically increases the depth of the transit.

The eccentricity and the arguments of periapsis of each planet are shown in figures 11 and 12, respectively. Over a timespan longer than the data, the eccentricities variations cause angular momenta to secularly swap between the planets, and the periastra of both planets make complete circuits around the star.

The mutual inclination of $2.40^{\circ}\pm 0.25^{\circ}$ causes a nodal precession with a period of $\sim 106$ years (shown for one model in Figure 13). This rate is much slower than the resonant precession rate (200 days), the 20-year $\Delta\varpi$ and eccentricity oscillation, and the 22-year periastron precession rate. Such a hierarchy of timescales has been noted for other resonant systems as well (Correia et al., 2010).

Figure 13 also shows impact parameter variations for the next 300 years, according to one model. It appears that planet c is more likely to be grazing (in the gray region) or missing the star (above the gray region) than planet b. Running 146 of the models described in Section 2.4 forwards 1000 years provided further evidence for this . Across the 146 models, planet c spent an average of $16\pm 2\%$ of the time grazing and $12\pm 3\%$ of the time missing; meanwhile, planet b had an average grazing fraction of $10\pm 3\%$ and missing fraction of $<5\%$ (at 95% confidence). Additionally, these results suggest that it is not rare in general for both of these planets to be observable at once.

We also confirmed that our solutions to the system are long-term stable. Some known systems have poorly characterized libration, or other sensitivity, such that not all solutions to the data are viable models. That is, since we are observing this planetary system that is likely Gyr old, we expect that it is in a stable state, and running it forward a few Myrs would show a regular pattern of motion. We selected 20 draws from the posterior and ran them forward in time with the Mercury integrator. They all survived for 10 Myr. On the one hand, this task suggests that all the models that are allowed by the data are viable; on the other hand, we cannot use a stability criterion to further constrain the system.

5 Discussion

The striking transit timing and depth variations of the K2-146 system can be compared to what was discovered in the Kepler prime mission. No duration changes as dramatic were found for planets orbiting single stars — only for circumbinary planets, whose precession and moving hosts caused impact parameters to vary from transit to transit. Regarding transit timing, Holczer et al. (2016) have measured the transit times of many of the planet candidates, computing transit timing amplitude of the dominant sinusoidal component in each planet. We compare the distributions of that amplitude and average planetary period to the K2-146 planets in Figure 14.

We see that the planets lie in an unusual spot in this space: being in resonance with large amplitude libration allows a large TTV amplitude to be visible, while having such a small orbital period allows us to see many cycles of the variation. This combination has led to the $\sim 1\%$ precision on $M_{p}/M_{\star}$ values. This privileged location is shared only with KOI-984.01, which is the only known transiting planet in its system. Thus it is very difficult to solve the KOI-984 system uniquely, and it is impossible to constrain the mass of KOI-984.01 to the level we have constrained the planets of K2-146.

5.1 Structure of the K2-146 planets

We turn to comparing the mass and radius constraints to other exoplanets and to theoretical models, depicted in Figure 15. The K2-146planets’ densities of $0.67\pm 0.04$ and $0.71\pm 0.05\rho_{\oplus}$ (or $3.69\pm 0.21$ and $3.92\pm 0.27g/cm^{3}$), respectively, are consistent with being either water worlds or having substantial atmospheres (Fortney et al., 2007; Lopez & Fortney, 2014, see also Figure 15). At first glance, this might be surprising given their short orbital periods. Owen & Wu (2013) proposed the existence of an “evaporation valley,” suggesting that planets with radii of $\sim 2$ R⊕ should be rare relative to planets with radius of 1.5 R⊕, which have had their atmospheres photoevaporated away, or 2.5 R⊕, which have strong enough gravity to retain their atmospheres. Such a valley was confirmed by Fulton et al. (2017); Van Eylen et al. (2018) expanded on this result, showing that the radius gap is a function of orbital period, in line with the predictions of Owen & Wu (2013).

Figure 16 shows the orbital period and radius of the K2-146 planets relative to the evaporation valley as inferred by Van Eylen et al. (2018). From those two parameters alone, it would appear that K2-146 b is below the gap and therefore should have lost its primordial atmosphere due to photoevaporation. However, these analyses used samples of FGK stars observed spectroscopically and asteroseismic targets, which are mostly more massive than the Sun. These planets, orbiting a mid-M dwarf, likely faced a very different UV environment over their first 100 Myr (e.g. Stelzer et al., 2013; Ansdell et al., 2015). There is a large scatter in the details of UV emission from M dwarf to M dwarf (Shkolnik & Barman, 2014), making it challenging to interpret the early UV environment for these planets. As a simple proxy, we can consider the bolometric incident flux received by the planet at the present time. If we reframe the evaporation valley in that context, as shown in Fig. 16, we see the planets are both comfortably above the observed gap, suggesting both planets could retain their atmospheres to the present day. The existence of K2-146 b in its present form strongly suggests that the specific stellar environment in which a planet resides should be considered when trying to interpret the history of its atmosphere and its susceptibility to photoevaporation.

5.2 Future Observations

K2-146 is faint ($V=16.2$) in the optical (Zacharias et al., 2012), limiting the possibility of RV follow-up with many facilities. However, the star is significantly brighter at redder wavelengths ($i=14.2,J=12.2$; Cutri et al., 2003; Zacharias et al., 2012), meaning current and planned RV facilities optimized to observe M dwarfs at redder wavelengths could be used to obtain additional observations of this system. RVs could also detect the presence of non-transiting planets in this system. This was the case for K2-18, where follow-up RV observations identified an inner non-transiting companion to a transiting planet (Montet et al., 2015; Cloutier et al., 2017). Additionally, while giant planets only orbit a few percent of mid-M dwarfs (Montet et al., 2014), the presence of a massive, outer planet could provide clues about the formation and dynamical evolution of these transiting planets (e.g. Otor et al., 2016; Gratia & Fabrycky, 2017).

While the star is faint, the observed transits of the inner planet are 0.4% in depth, making these transits visible from the ground with modest telescopes. RVs and additional transit photometry both provide future opportunities on reasonable timescales to determine which plausible solution accurately describes the K2-146 system. The errors on planet properties are already driven by the uncertainties on stellar characterization. The main motivation to continue measuring the planetary periods is to refine the ephemerides of when future transits will occur. This would be important if a campaign is envisioned for probing the thin atmospheres on these planets. Predicted transit times and uncertainties until September 2022 are available for digital download, and a preview is shown in Tables 7 and 8 in the Appendix; ephemerides can be improved most efficiently by observing future transits with the highest timing uncertainties. We expect that both planets would have JWST S/N below the cutoff recommended by Kempton et al. (2018) for atmospheric characterization. However, because their masses are known more precisely than is usually possible with RV measurements, they may be valuable candidates for atmospheric study nonetheless.

5.3 The K2 Legacy Field as a TESS Testbed

K2-146 demonstrates the power of the K2 mission, and space-based transit missions in general, for detailed characterization of planetary systems through repeated observations of fields. Without any ground-based follow up observations, we have succeeded in measuring the mass of both planets to a precision of 3% despite the Kepler spacecraft observing this system for less than 20% of the observational baseline. In this case, the Campaign 5 observations were fortuitously timed so that the presence of TTVs could be inferred from that data set alone, but longer-baseline TTVs in other systems will only be detectable by combining data from Campaign 5, 16, and 18 together. Moreover, continued observations of previously observed fields enable the detection of new, previously undiscovered planets. This is both true in the case of K2-146 c, where the planet’s transit depth increased between campaigns, and in the likely more typical case where additional observations raise the S/N of a phase-folded transit signal above the criterion for significance. This system, and the K2 legacy field in general, thus serves as a dress rehearsal for the possibilities that a TESS extended mission will enable, when transits are observed in sectors separated by long data gaps. In future work, we will investigate other systems in this field, including a search for previously undetected planets and TTVs.

We thank Leslie Rogers for discussing the nature of gaseous envelopes on small planets and Rodrigo Luger for help with the EVEREST pipeline. A.C.H. was supported by a K2GO5 grant from NASA, number 80NSSC18K0364. Work by B.T.M. was performed under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute.

This paper includes data collected by the K2 mission. Funding for the K2 mission is provided by the NASA Science Mission directorate. We are grateful to the entire Kepler and K2 team, past and present. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Some/all of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

Kepler, Gaia, Exoplanet Archive