Discovery of a Third Transiting Planet in the Kepler-47 Circumbinary System

TL;DR

This paper reports the discovery of a third transiting planet in the Kepler-47 circumbinary system, providing improved mass constraints and insights into the system's formation and stability.

Contribution

The discovery of a third planet in Kepler-47 and the refined mass and density estimates offer new understanding of circumbinary planetary systems.

Findings

Discovery of Kepler-47 d with a 187.4-day orbit and 7 R⊕ radius

Mass constraints for all three planets are significantly improved

Outer planets are low-density and tightly packed, indicating stable, gentle migration

Abstract

Of the nine confirmed transiting circumbinary planet systems, only Kepler-47 is known to contain more than one planet. Kepler-47 b (the "inner planet") has an orbital period of 49.5 days and a radius of about $3\,R_{\oplus}$. Kepler-47 c (the "outer planet") has an orbital period of 303.2 days and a radius of about $4.7\,R_{\oplus}$. Here we report the discovery of a third planet, Kepler-47 d (the "middle planet"), which has an orbital period of 187.4 days and a radius of about $7\,R_{\oplus}$. The presence of the middle planet allows us to place much better constraints on the masses of all three planets, where the $1\sigma$ ranges are less than $26\,M_{\oplus}$, between $7-43\,M_{\oplus}$, and between $2-5\,M_{\oplus}$ for the inner, middle, and outer planets, respectively. The middle and outer planets have low bulk densities, with $\rho_{\rm middle} < 0.68$ g cm$^{-3}$ and $\rho_{\rm outer} < 0.26$ g cm$^{-3}$ at the $1\sigma$ level. The two outer planets are "tightly packed," assuming the nominal masses, meaning no other planet could stably orbit between them. All of the orbits have low eccentricities and are nearly coplanar, disfavoring violent scattering scenarios and suggesting gentle migration in the protoplanetary disk.

Tables

Table 1. Table 1: Radial velocities for Kepler-47

Date	UT Time	BJD	RVA	RVA	Telescope
YYYY-MM-DD		(2,455,000+)	(km s-1)aaUncorrected radial velocities.	(km s-1)bbCorrected radial velocities with scaled error bars.
2012 Apr 10	13:25:48.68	1028.05942	$11.442\pm 0.011$	$\cdots$	Keck
2012 Apr 23	09:11:27.36	1040.90325	$33.534\pm 0.091$	$28.934\pm 0.164$	HET
2012 May 1	09:52:55.08	1048.93237	$35.458\pm 0.171$	$30.890\pm 0.308$	HJST
2012 May 2	07:23:45.95	1049.82882	$24.430\pm 0.440$	$19.862\pm 0.792$	HJST
2012 May 4	08:34:10.40	1051.88474	$-21.957\pm 0.159$	$-26.525\pm 0.286$	HJST
2012 May 5	08:08:26.83	1052.86692	$-26.719\pm 0.178$	$-31.287\pm 0.320$	HJST
2012 May 6	08:08:55.42	1053.86729	$-9.150\pm 0.122$	$-13.718\pm 0.220$	HJST
2012 May 18	07:35:21.15	1065.83749	$-1.843\pm 0.060$	$-6.443\pm 0.108$	HET
2012 May 20	07:37:07.78	1067.83880	$-25.681\pm 0.030$	$-30.281\pm 0.054$	HET
2012 Jun 5	06:29:24.15	1083.79236	$-5.223\pm 0.080$	$-9.823\pm 0.144$	HET
2012 Jun 26	08:03:52.39	1104.85862	$-26.743\pm 0.086$	$-31.311\pm 0.155$	HJST
2012 Aug 20	02:04:42.20	1159.5866	$26.670\pm 0.080$	$22.402\pm 0.144$	1.93m OHP
2012 Aug 21	01:53:19.7	1160.5787	$34.880\pm 0.130$	$30.612\pm 0.234$	1.93m OHP
2012 Aug 22	01:10:33.6	1161.5490	$24.220\pm 0.080$	$19.952\pm 0.144$	1.93m OHP
2012 Sep 7	19:30:08.6	1178.3126	$-19.070\pm 0.160$	$-23.339\pm 0.288$	1.93m OHP

Table 2. Table 2: Times and Durations of Planetary Transits

	Measured	Uncertainty	Duration	Model	Model	Impact
cycle#	TimeaaBJD-2,455,000	(minute)	(hr)	Timea	Duration (hr)	Parameter	Note
Inner Planet
1.0	-30.80969	4.46	$4.55\pm 0.19$	-30.80876	4.76	0.29
2.0	…	…	…	16.27666	4.10	0.30	data gap
3.0	65.24832	22.29	$5.81\pm 1.40$	65.24587	6.41	0.29
4.0	112.54179	10.51	$2.94\pm 0.69$	112.54644	3.76	0.33
5.0	160.89953	32.75	$10.83\pm 2.66$	160.90621	10.45	0.28
6.0	208.84796	7.54	$3.47\pm 0.66$	208.83932	3.69	0.36
7.0	256.34222	11.56	$8.86\pm 0.69$	256.32139	7.90	0.30
8.0	305.14101	7.61	$4.97\pm 0.61$	305.12131	3.85	0.38
9.0	352.25457	12.36	$5.41\pm 1.90$	352.25805	5.10	0.35
10.0	401.35072	11.18	$4.16\pm 0.50$	401.35481	4.34	0.40
11.0	448.43586	9.83	$3.16\pm 0.60$	448.43157	4.06	0.40
12.0	497.42556	12.52	$2.96\pm 0.89$	497.46945	5.54	0.41
13.0	544.73196	17.72	$6.90\pm 0.82$	544.69458	3.61	0.44
14.0	593.25486	8.75	$8.75\pm 0.37$	593.27974	8.76	0.41
15.0	…	…	…	640.98828	3.45	0.48	data gap
16.0	688.66867	11.27	$8.47\pm 0.47$	688.65503	8.74	0.42
17.0	737.27699	37.31	$2.99\pm 1.23$	737.27888	3.51	0.51
18.0	784.46811	19.36	$4.84\pm 3.83$	784.47271	5.25	0.47
19.0	…	…	…	833.53075	3.85	0.53	data gap
20.0	880.61567	11.00	$4.90\pm 0.99$	880.60230	3.95	0.52
21.0	…	…	…	929.68734	4.72	0.54	corrupted data
22.0	976.86285	11.31	$4.08\pm 0.72$	976.84585	3.40	0.56
23.0	1025.61184	12.69	$7.25\pm 0.53$	1025.60438	7.09	0.53
24.0	…	…	…	1073.13111	3.18	0.59	blend
25.0	1120.93797	70.76	$7.42\pm 2.54$	1121.01153	9.26	0.53
26.0	…	…	…	1169.41978	3.19	0.61	data gap
27.0	1216.66103	15.68	$5.72\pm 0.84$	1216.67569	5.38	0.57
28.0	1265.68597	20.38	$2.50\pm 0.75$	1265.67744	3.47	0.62
29.0	…	…	…	1312.74767	3.85	0.60	data gap
30.0	1361.87205	11.00	$4.14\pm 1.12$	1361.85261	4.20	0.61
31.0	1408.96020	12.58	$3.00\pm 1.37$	1408.96429	3.25	0.63
Middle Planet
1.0	…	…	…	45.15444	0.90	0.99	blend
2.0	…	…	…	231.57791	1.33	0.97	data gap
3.0	418.04185	8.93	$1.71\pm 0.86$	418.00486	1.66	0.95	one cadence missing
4.0	604.43883	5.69	$1.12\pm 0.46$	604.44003	1.97	0.94
5.0	790.88576	16.08	$1.55\pm 0.71$	790.89628	2.31	0.93
6.0	977.36103	4.98	$1.86\pm 0.65$	977.35770	2.73	0.91	original orphan
7.0	1163.84405	3.35	$2.99\pm 0.21$	1163.84154	3.31	0.90
8.0	1350.36851	6.08	$4.81\pm 0.60$	1350.36740	4.23	0.88
Outer Planet
1.0	246.65182	5.95	$4.14\pm 0.64$	246.65033	4.29	0.70
2.0	550.47424	5.08	$5.96\pm 0.39$	550.47748	6.16	0.56
3.0	850.97978	11.50	$8.07\pm 1.46$	850.98582	5.80	0.72
4.0	1154.79184	6.42	$6.15\pm 0.68$	1154.78852	4.95	0.55

Table 3. Table 3: Times of stellar eclipses

cycle #	primary	corrected	uncertainty	cycle #	secondary	uncertainty
	timeaaBJD-2,455,000	timeaaBJD-2,455,000	(min)		timeaaBJD-2,455,000	(min)
0.0	…	…	…	0.487385	$-33.12242$	2.55
1.0	$-29.30551$	$-29.30557$	0.51	1.487385	$-25.67575$	2.65
2.0	$-21.85721$	$-21.85712$	0.51	2.487385	$-18.22955$	2.65
3.0	…	…	…	3.487385	$-10.77950$	2.54
4.0	$-6.96069$	$-6.96058$	0.49	4.487385	…	…
5.0	…	…	…	5.487385	$4.11925$	3.00
6.0	$7.93602$	$7.93622$	0.48	6.487385	$11.56766$	2.89
7.0	…	…	…	7.487385	$19.01206$	2.77
8.0	$22.83277$	$22.83289$	0.48	8.487385	$26.46528$	2.65
9.0	$30.28122$	$30.28145$	0.48	9.487385	$33.91418$	2.31
10.0	$37.72961$	$37.72978$	0.48	10.487385	$41.36150$	2.65
11.0	…	…	…	11.487385	$48.80864$	2.54
12.0	…	…	…	12.487385	$56.25788$	2.77
13.0	$60.07493$	$60.07506$	0.49	13.487385	…	…
14.0	$67.52346$	$67.52350$	0.50	14.487385	$71.15074$	2.42
15.0	$74.97166$	$74.97165$	0.49	15.487385	$78.60188$	2.65
16.0	$82.42033$	$82.42022$	0.49	16.487385	$86.05079$	2.54
17.0	$89.86868$	$89.86856$	0.49	17.487385	$93.49712$	3.12
18.0	$97.31711$	$97.31693$	0.47	18.487385	$100.94813$	2.54
19.0	…	…	…	19.487385	$108.39654$	2.54
20.0	$112.21389$	$112.21356$	0.48	20.487385	$115.84383$	2.42
21.0	$119.66226$	$119.66203$	0.48	21.487385	…	…
22.0	$127.11047$	$127.11021$	0.48	22.487385	$130.73682$	2.54
23.0	$134.55888$	$134.55875$	0.50	23.487385	…	…
24.0	$142.00720$	$142.00713$	0.49	24.487385	$145.63577$	2.77
25.0	$149.45537$	$149.45534$	0.48	25.487385	…	…
26.0	$156.90382$	$156.90398$	0.48	26.487385	$160.53487$	2.65
27.0	$164.35231$	$164.35226$	0.48	27.487385	$167.98122$	2.54
28.0	$171.80081$	$171.80089$	0.49	28.487385	…	…
29.0	$179.24881$	$179.24907$	0.48	29.487385	…	…
30.0	$186.69701$	$186.69743$	0.48	30.487385	$190.32692$	2.54
31.0	$194.14522$	$194.14566$	0.48	31.487385	$197.77646$	2.65
32.0	$201.59352$	$201.59381$	0.50	32.487385	$205.22630$	2.54
33.0	…	…	…	33.487385	$212.67039$	3.46
34.0	…	…	…	34.487385	$220.12045$	2.42
35.0	$223.93909$	$223.93924$	0.50	35.487385	$227.56808$	3.46
36.0	…	…	…	36.487385	$235.01796$	2.42
37.0	$238.83626$	$238.83617$	0.48	37.487385	$242.46574$	2.42
38.0	$246.28459$	$246.28468$	0.49	38.487385	$249.91273$	2.54
39.0	$253.73305$	$253.73320$	0.47	39.487385	$257.36230$	2.54
40.0	$261.18094$	$261.18106$	0.47	40.487385	$264.81024$	2.42
41.0	$268.62946$	$268.62949$	0.48	41.487385	…	…
42.0	…	…	…	42.487385	$279.70694$	2.54
43.0	$283.52628$	$283.52618$	0.49	43.487385	…	…
44.0	$290.97491$	$290.97467$	0.49	44.487385	$294.60605$	2.42
45.0	$298.42331$	$298.42300$	0.48	45.487385	…	…
46.0	$305.87186$	$305.87149$	0.49	46.487385	$309.50082$	2.65
47.0	$313.32040$	$313.31982$	0.48	47.487385	…	…
48.0	$320.76880$	$320.76819$	0.49	48.487385	$324.39990$	2.42
49.0	$328.21707$	$328.21655$	0.48	49.487385	$331.84735$	2.42
50.0	$335.66525$	$335.66498$	0.49	50.487385	$339.29431$	2.42
51.0	$343.11343$	$343.11331$	0.48	51.487385	$346.74100$	2.54
52.0	$350.56171$	$350.56155$	0.50	52.487385	$354.19067$	2.65
53.0	$358.00989$	$358.00983$	0.49	53.487385	$361.63971$	2.54
54.0	$365.45822$	$365.45825$	0.48	54.487385	$369.08972$	2.54
55.0	$372.90720$	$372.90741$	0.63	55.487385	$376.53705$	2.65
56.0	…	…	…	56.487385	$383.98596$	2.54
57.0	$387.80347$	$387.80380$	0.47	57.487385	$391.43390$	2.65
58.0	$395.25180$	$395.25220$	0.48	58.487385	$398.88248$	2.65
59.0	$402.69992$	$402.70035$	0.50	59.487385	$406.33072$	2.77
60.0	$410.14825$	$410.14865$	0.49	60.487385	…	…
61.0	$417.59674$	$417.59717$	0.49	61.487385	…	…
62.0	$425.04520$	$425.04559$	0.50	62.487385	$428.67667$	2.42
63.0	$432.49371$	$432.49371$	0.92	63.487385	$436.12863$	3.70
64.0	$439.94229$	$439.94229$	0.48	64.487389	$443.57419$	2.54
65.0	$447.39066$	$447.39035$	0.49	65.487389	$451.01895$	2.65
66.0	$454.83932$	$454.83868$	0.50	66.487389	$458.47025$	2.65
67.0	…	…	…	67.487389	…	…
68.0	$469.73666$	$469.73569$	0.51	68.487389	$473.36713$	2.42
69.0	…	…	…	69.487389	$480.81107$	2.54
70.0	$484.63339$	$484.63257$	0.50	70.487389	$488.26462$	2.54
71.0	…	…	…	71.487389	$495.71335$	2.42
72.0	…	…	…	72.487389	$503.15891$	2.65
73.0	$506.97739$	$506.97763$	0.49	73.487389	$510.60684$	2.77
74.0	$514.42554$	$514.42584$	0.49	74.487389	$518.05542$	2.77
75.0	$521.87421$	$521.87439$	0.47	75.487389	$525.50653$	2.54
76.0	$529.32257$	$529.32281$	0.48	76.487389	$532.95557$	2.65
77.0	$536.77094$	$536.77118$	0.47	77.487389	$540.40015$	2.54
78.0	…	…	…	78.487389	$547.84821$	2.77
79.0	$551.66760$	$551.66797$	0.48	79.487389	…	…
80.0	…	…	…	80.487389	…	…
81.0	…	…	…	81.487389	…	…
82.0	$574.01306$	$574.01324$	0.48	82.487389	$577.64178$	2.65
83.0	$581.46143$	$581.46149$	0.48	83.487389	…	…
84.0	…	…	…	84.487389	$592.54077$	2.54
85.0	…	…	…	85.487389	$599.98920$	2.42
86.0	$603.80609$	$603.80603$	0.54	86.487389	$607.43634$	2.54
87.0	$611.25482$	$611.25470$	0.52	87.487389	$614.88654$	2.42
88.0	$618.70325$	$618.70300$	0.52	88.487389	$622.33221$	2.65
89.0	$626.15149$	$626.15137$	0.52	89.487389	$629.78101$	2.77
90.0	$633.60010$	$633.60004$	0.51	90.487389	…	…
91.0	…	…	…	91.487389	$644.67877$	2.54
92.0	$648.49695$	$648.49652$	0.47	92.487389	$652.12604$	2.54
93.0	$655.94531$	$655.94482$	0.47	93.487389	$659.57526$	2.54
94.0	$663.39392$	$663.39337$	0.48	94.487389	$667.02484$	2.42
95.0	$670.84247$	$670.84210$	0.48	95.487389	$674.47101$	2.54
96.0	…	…	…	96.487389	$681.92090$	2.54
97.0	$685.73883$	$685.73865$	0.48	97.487389	$689.37140$	2.65
98.0	$693.18701$	$693.18695$	0.49	98.487389	$696.81665$	2.31
99.0	$700.63531$	$700.63531$	0.49	99.487389	$704.26617$	2.54
100.0	$708.08374$	$708.08374$	0.49	100.487389	$711.71442$	2.42
101.0	$715.53192$	$715.53204$	0.48	101.487389	$719.16046$	2.42
102.0	$722.98016$	$722.98029$	0.48	102.487389	$726.60980$	2.42
103.0	$730.42853$	$730.42883$	0.47	103.487389	$734.05872$	2.54
104.0	$737.87714$	$737.87732$	0.47	104.487389	$741.50854$	2.42
105.0	$745.32550$	$745.32550$	0.48	105.487389	$748.95648$	2.54
106.0	$752.77393$	$752.77380$	0.47	106.487389	$756.40527$	2.65
107.0	…	…	…	107.487389	$763.85193$	2.65
108.0	$767.67059$	$767.67065$	0.47	108.487389	$771.30164$	2.65
109.0	$775.11890$	$775.11896$	0.48	109.487389	…	…
110.0	$782.56732$	$782.56732$	0.51	110.487389	$786.19946$	2.65
111.0	…	…	…	111.487389	…	…
112.0	$797.46442$	$797.46448$	0.49	112.487389	$801.09283$	2.54
113.0	$804.91309$	$804.91296$	0.50	113.487389	$808.54272$	2.42
114.0	$812.36139$	$812.36102$	0.47	114.487389	$815.99078$	2.77
115.0	$819.80957$	$819.80927$	0.47	115.487389	$823.43982$	2.54
116.0	$827.25800$	$827.25775$	0.47	116.487389	$830.88733$	2.65
117.0	$834.70618$	$834.70618$	0.48	117.487389	$838.33765$	2.65
118.0	…	…	…	118.487389	…	…
119.0	$849.60260$	$849.60272$	0.47	119.487389	$853.23230$	2.42
120.0	$857.05133$	$857.05127$	0.49	120.487389	$860.68353$	2.42
121.0	$864.49988$	$864.49969$	0.49	121.487389	$868.13147$	2.42
122.0	$871.94800$	$871.94781$	0.47	122.487389	$875.57574$	2.54
123.0	$879.39630$	$879.39630$	0.47	123.487389	$883.02753$	2.54
124.0	$886.84467$	$886.84460$	0.47	124.487389	$890.47528$	2.31
125.0	$894.29315$	$894.29303$	0.47	125.487389	$897.92639$	2.42
126.0	$901.74152$	$901.74146$	0.48	126.487389	…	…
127.0	$909.19000$	$909.19006$	0.48	127.487389	$912.82153$	2.77
128.0	$916.63831$	$916.63855$	0.49	128.487381	$920.26917$	2.54
129.0	$924.08630$	$924.08685$	0.48	129.487381	…	…
130.0	…	…	…	130.487381	$935.16748$	2.54
131.0	$938.98212$	$938.98297$	0.53	131.487381	$942.61377$	2.42
132.0	…	…	…	132.487381	…	…
133.0	$953.87927$	$953.87964$	0.55	133.487381	…	…
134.0	$961.32800$	$961.32831$	0.57	134.487381	$964.95990$	3.12
135.0	$968.77692$	$968.77704$	0.57	135.487381	$972.40747$	2.42
136.0	$976.22607$	$976.22589$	0.58	136.487381	$979.85541$	2.65
137.0	$983.67474$	$983.67450$	0.52	137.487381	…	…
138.0	$991.12311$	$991.12274$	0.51	138.487381	…	…
139.0	…	…	…	139.487381	$1002.20215$	2.31
140.0	$1006.01910$	$1006.01849$	0.51	140.487381	…	…
141.0	…	…	…	141.487381	$1017.09631$	2.42
142.0	…	…	…	142.487381	$1024.54651$	2.54
143.0	$1028.36389$	$1028.36365$	0.49	143.487381	$1031.99353$	2.42
144.0	$1035.81226$	$1035.81226$	0.48	144.487381	$1039.44019$	2.54
145.0	$1043.26074$	$1043.26074$	0.49	145.487381	$1046.89429$	3.23
146.0	…	…	…	146.487381	$1054.33911$	2.65
147.0	$1058.15735$	$1058.15723$	0.48	147.487381	…	…
148.0	$1065.60583$	$1065.60571$	0.49	148.487381	$1069.23608$	2.54
149.0	$1073.05432$	$1073.05444$	0.50	149.487381	$1076.68445$	2.31
150.0	$1080.50232$	$1080.50269$	0.50	150.487381	$1084.13232$	2.42
151.0	$1087.95056$	$1087.95093$	0.50	151.487381	$1091.58301$	2.54
152.0	$1095.39880$	$1095.39929$	0.50	152.487381	$1099.02747$	2.54
153.0	$1102.84717$	$1102.84753$	0.48	153.487381	…	…
154.0	$1110.29541$	$1110.29578$	0.47	154.487381	$1113.92444$	2.43
155.0	…	…	…	155.487381	$1121.37476$	2.32
156.0	…	…	…	156.487381	…	…
157.0	$1132.64062$	$1132.64111$	0.47	157.487381	$1136.27185$	2.43
158.0	…	…	…	158.487381	$1143.71802$	2.32
159.0	…	…	…	159.487381	$1151.16760$	2.21
160.0	…	…	…	160.487381	$1158.61731$	2.43
161.0	$1162.43567$	$1162.43506$	0.47	161.487381	$1166.06384$	2.55
162.0	…	…	…	162.487381	$1173.51453$	2.21
163.0	…	…	…	163.487381	$1180.96008$	2.32
164.0	$1184.78064$	$1184.77991$	0.47	164.487381	$1188.41089$	2.43
165.0	$1192.22864$	$1192.22791$	0.47	165.487381	$1195.85779$	2.43
166.0	$1199.67688$	$1199.67639$	0.47	166.487381	$1203.30750$	2.43
167.0	$1207.12500$	$1207.12476$	0.47	167.487381	$1210.75415$	2.32
168.0	…	…	…	168.487381	$1218.20361$	2.43
169.0	…	…	…	169.487381	$1225.65259$	2.21
170.0	…	…	…	170.487381	$1233.09839$	2.32
171.0	…	…	…	171.487381	$1240.54895$	2.55
172.0	$1244.36633$	$1244.36682$	0.47	172.487381	…	…
173.0	…	…	…	173.487381	$1255.44543$	2.32
174.0	$1259.26306$	$1259.26343$	0.47	174.487381	$1262.89319$	2.66
175.0	$1266.71155$	$1266.71179$	0.47	175.487381	$1270.33948$	2.21
176.0	$1274.15991$	$1274.16028$	0.47	176.487381	$1277.78967$	2.32
177.0	$1281.60815$	$1281.60840$	0.47	177.487381	$1285.24048$	2.55
178.0	$1289.05664$	$1289.05701$	0.47	178.487381	$1292.68555$	2.21
179.0	$1296.50476$	$1296.50513$	0.47	179.487381	$1300.13538$	2.21
180.0	…	…	…	180.487381	$1307.58350$	2.32
181.0	…	…	…	181.487381	…	…
182.0	…	…	…	182.487381	$1322.48193$	2.43
183.0	$1326.29883$	$1326.29883$	0.47	183.487381	$1329.93335$	2.43
184.0	…	…	…	184.487381	…	…
185.0	$1341.19531$	$1341.19568$	0.47	185.487381	$1344.82397$	2.43
186.0	$1348.64368$	$1348.64404$	0.47	186.487381	$1352.27368$	2.43
187.0	$1356.09229$	$1356.09241$	0.47	187.487381	$1359.72375$	2.43
188.0	…	…	…	188.487381	$1367.17090$	2.32
189.0	$1370.98950$	$1370.98914$	0.47	189.487381	$1374.61914$	2.32
190.0	$1378.43774$	$1378.43726$	0.47	190.487381	$1382.06519$	2.66
191.0	…	…	…	191.487381	$1389.51453$	2.21
192.0	$1393.33435$	$1393.33398$	0.47	192.487381	$1396.96460$	2.21
193.0	$1400.78259$	$1400.78247$	0.47	193.487381	$1404.41370$	2.21
194.0	$1408.23083$	$1408.23096$	0.47	194.487381	$1411.86157$	2.21
195.0	…	…	…	…	…	…
196.0	$1423.12757$	$1423.12767$	0.47	…	…	…

Table 4. Table 4: Priors for fitting parameters

Orbit	Parameter	Lower	Upper
Binary and Stellar	Period (days)	7.4480	7.44899
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$-29.316$	$-29.300$
	$e\cos \omega$	$-0.1$	$0.0$
	$e\sin \omega$	$-0.05$	$0.01$
	Inclination (deg)	89.00	90.25
	$M_1$ ($M_{\odot }$)	0.80	1.0
	$Q$	0.340	0.373
	$R_1/a$	0.046	0.061
	$R_2/a$	0.0186	0.0202
	$T_2/T_1$	0.39	0.83
	$q_1$ (primary)	0.0	1.0
	$q_2$ (primary)	0.0	1.0
	$q_1$ (secondary)	0.0	1.0
	$q_2$ (secondary)	0.0	1.0
Inner	Period (days)	49.41	49.52
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$-31.6$	$-31.14$
	$e\cos \omega$	$-0.08$	$0.10$
	$e\sin \omega$	$10.114$	$0.069$
	Inclination (deg)	89.03	90.32
	Nodal angle (deg)	$-0.390$	0.315
	$M_3$ ($M_{\oplus }$)	0.01	100.0
Middle	Period (days)	187.15	187.60
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$45.10$	$45.75$
	$e\cos \omega$	$-0.39$	$0.35$
	$e\sin \omega$	$-0.1$	$0.2$
	Inclination (deg)	90.18	90.62
	Nodal angle (deg)	$-2.8$	$0.5$
	$M_3$ ($M_{\oplus }$)	0.01	100.0
Outer	Period (days)	303.00	303.42
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$-55.8$	$-55.4$
	$e\cos \omega$	$-0.19$	$0.23$
	$e\sin \omega$	$-0.21$	$0.11$
	Inclination (deg)	90.10	90.28
	Nodal angle (deg)	$-3.0$	$1.2$
	$M_3$ ($M_{\oplus }$)	0.01	100.0
Kepler	First season contaminationaaKepler Quarters Q1, Q5, …	0.000	0.035
	Second season contaminationbbKepler Quarters Q2, Q6, …	0.000	0.035
	Third season contaminationccKepler Quarters Q3, Q7, …	0.000	0.035
	Fourth season contaminationddKepler Quarters Q4, Q8, …	0.000	0.035

Table 5. Table 5: Fitted Parameters from the Photometric-dynamical model a a T ref = subscript 𝑇 ref absent T_{\rm ref}= BJD 2,454,965.000

Parameter	Binary Orbit	Inner Orbit	Middle Orbit	Outer Orbit
Period (days)	${7.4483648}_{-0.0000270}^{+0.0000038}$	${49.4643}_{-0.0074}^{+0.0081}$	${187.366}_{0.051}^{+0.069}$	${303.227}_{-0.027}^{+0.062}$
$T_{\mathrm{conj}}$ (BJD + 2,455,000)	$-{29.305248}_{-0.000030}^{+0.000031}$	$-{31.3573}_{-0.0049}^{+0.0043}$	${45.428}_{-0.047}^{+0.044}$	$-{55.638}_{-0.046}^{+0.067}$
$e\cos \omega$	$-{0.019895}_{-0.000024}^{+0.000025}$	${0.0139}_{-0.0025}^{+0.0020}$	${0.024}_{-0.036}^{+0.038}$	${0.026}_{-0.034}^{+0.029}$
$e\sin \omega$	$-{0.0208}_{-0.0020}^{+0.0019}$	${0.0157}_{-0.0011}^{+0.0015}$	$-{0.004}_{-0.028}^{+0.023}$	$-{0.035}_{-0.025}^{+0.026}$
$i$ (deg)	${89.613}_{-0.040}^{+0.045}$	${89.752}_{-0.045}^{+0.063}$	${90.395}_{-0.012}^{+0.009}$	${90.1925}_{-0.0042}^{+0.0055}$
$\Omega$ (deg)	$0.0$	$-{0.088}_{-0.053}^{+0.056}$	$-{0.87}_{-0.10}^{+0.11}$	$-{1.28}_{-0.25}^{+0.23}$
$Q_{2,1}\equiv M_2/M_1$	${0.3569}_{-0.0020}^{+0.0021}$	…	…	…
$M_i$ ($M_{\oplus }$)	…	${2.07}_{-2.07}^{+23.70}$	${19.02}_{-11.67}^{+23.84}$	${3.17}_{-1.25}^{+2.18}$
parameter	primary	secondary	inner planet	middle planet	outer planet
$R/a$bbFractional stellar radii	${0.05347}_{-0.00016}^{+0.00016}$	${0.019332}_{-0.000065}^{+0.000056}$	…	…	…
$R_1/R$ccRatio of primary radius to planet radius	…	…	${33.51}_{-0.48}^{+0.54}$	${14.55}_{-1.35}^{+1.11}$	${21.97}_{-0.35}^{+0.42}$
$q_1$ddKipping’s first quad-law triangular limb darkening coefficient	${0.502}_{-0.063}^{+0.067}$		…	…	…
$q_2$eeKipping’s second quad-law triangular limb darkening coefficient	${0.188}_{-0.028}^{+0.020}$	${0.38}_{-0.12}^{+0.34}$	…	…	…
Other parameters
$T_2/T_1$	${0.613}_{-0.006}^{+0.012}$	…	…	…	…
$s_0$ (Season 0 contamination)		…	…	…	…
$s_1$ (Season 1 contamination)	${0.0219}_{-0.0037}^{+0.0045}$	…	…	…	…
$s_2$ (Season 2 contamination)	${0.0044}_{-0.0038}^{+0.0043}$	…	…	…	…
$s_3$ (Season 3 contamination)	${0.0232}_{-0.0029}^{+0.0047}$	…	…	…	…

Table 6. Table 6: Derived Parameters from the Photometric-dynamical Model

Parameter	Binary Orbit	Inner Orbit	Middle Orbit	Outer Orbit
$a$ (au)	${0.08145}_{-0.00037}^{+0.00036}$	${0.2877}_{-0.0011}^{+0.0014}$	${0.6992}_{-0.0033}^{+0.0031}$	${0.9638}_{-0.0044}^{+0.0041}$
$e$	${0.0288}_{-0.0013}^{+0.0015}$	${0.0210}_{-0.0022}^{+0.0025}$	${0.024}_{-0.017}^{+0.025}$	${0.044}_{-0.019}^{+0.029}$
$\omega$ (deg)	${226.3}_{-2.6}^{+2.8}$	${48.6}_{-2.3}^{+3.0}$	${352}_{-166}^{+96}$	${306}_{-67}^{+40}$
$I_{1,i}$ (deg)aaMutual inclination between the orbital planes of the binary and the planets, given by $\cos I=\sin i_p\sin i_b\cos \Delta \Omega +\cos i_p\cos i_b$	…	${0.166}_{-0.041}^{+0.037}$	${1.165}_{-0.093}^{+0.101}$	${1.38}_{-0.24}^{+0.20}$
$I_{2,i}$ (deg)bbMutual inclination between the orbital planes of the inner planet and the other planets	…	…	${1.006}_{-0.050}^{+0.089}$	${1.247}_{-0.018}^{+0.024}$
$I_{3,i}$ (deg)ccMutual inclination between the orbital planes of the middle planet and the outer planet	…	…	…	${0.442}_{-0.102}^{+0.028}$
Parameter	Primary	Secondary	Inner Planet	Middle Planet	Outer Planet
$M$ ($M_{\odot }$ or $M_{\oplus }$)	${0.957}_{-0.015}^{+0.013}$	${0.342}_{-0.003}^{+0.003}$	${2.07}_{-2.07}^{+23.70}$	${19.02}_{-11.67}^{+23.84}$	${3.17}_{-1.25}^{+2.18}$
$R$ ($R_{\odot }$ or $R_{\oplus }$)	${0.936}_{-0.005}^{+0.005}$	${0.338}_{-0.002}^{+0.002}$	${3.05}_{-0.04}^{+0.04}$	${7.04}_{-0.49}^{+0.66}$	${4.65}_{-0.07}^{+0.09}$
$\rho$ (g cm-3)	${1.65}_{-0.01}^{+0.02}$	${12.42}_{-0.11}^{+0.11}$	${0.40}_{-0.40}^{+4.87}$	${0.30}_{-0.18}^{+0.38}$	${0.17}_{-0.07}^{+0.09}$

Table 7. Table 7: Initial Cartesian Barycentric Conditions a for Kepler-47

Parameter	Primary	Secondary	Inner Planet	Middle Planet	Outer Planet
Mass ($M_{\odot }$)	9.573912601143495E$-$01	3.416626576742950E$-$01	6.208164703811071E$-$06	5.711706000000001E$-$05	9.530601983442541E$-$06
$x$ (au)	2.092317942847221E$-$02	$-5$.866162141774457E$-$02	1.302102743834203E$-$01	2.398289447308760E$-$01	$-3$.872087634688049E$-$01
$y$ (au)	$-1$.364205200783173E$-$05	3.791628623967441E$-$05	8.812713793817102E$-$04	8.391663271667140E$-$04	5.540095504781927E$-$03
$z$ (au)	$-1$.988193746560007E$-$03	5.649120240556342E$-$03	2.500211844794519E$-$01	$-6$.478290001132635E$-$01	9.270974333301998E$-$01
$v_x$ (au)	1.350420386068420E$-$03	$-3$.786681230811973E$-$03	$-3$.299902928667536E$-$02	2.207346827176963E$-$02	$-1$.773342551338714E$-$02
$v_y$ (au)	1.239137665559171E$-$04	$-3$.471736278879929E$-$04	1.256761297656883E$-$04	$-3$.936013341278077E$-$04	4.205823582940745E$-$04
$v_z$ (au)	1.835286133723581E$-$02	$-5$.142910738078619E$-$02	1.739447523236208E$-$02	8.701783937203954E$-$03	$-7$.191666020340975E$-$03

Table 8. Table 8: Initial Orbital Elements a for Kepler-47

parameter	binary	inner	middle	outer
	orbit	orbit	orbit	orbit
$a$ (AU)	8.144243202775878E$-$02	2.877434527504005E$-$01	6.992154891162726E$-$01	9.638218771090775E$-$01
$e$	2.877193854136767E$-$02	2.095514505946130E$-$02	2.416231575331903E$-$02	4.371716413425727E$-$02
$i$ (deg)	8.961321083430335E$+$01	8.975235828571519E$+$01	9.039477337702283E$+$01	9.019249877520613E$+$01
$\omega$ (deg)	2.262537367634179E$+$02	4.859496576626187E$+$01	3.517243342867096E$+$02	3.059357996708584E$+$02
$\Omega$ (deg)	$0.$000000000000000E$+$00	$-8.$764673739309498E$-$02	$-8.$659498733670760E$-$01	$-1.$280580223423379E$+$00
$f$ (deg)	3.082645652441405E$+$02	1.389173589110049E$+$01	2.985868391504580E$+$02	1.667337092860433E$+$02
$M$ (deg)	3.108188471954695E$+$02	1.332389759998393E$+$01	3.009969996738955E$+$02	1.655463689900479E$+$02
$\lambda$ (deg)	1.745183020075583E$+$02	6.239905491996927E$+$01	2.894452235638005E$+$02	1.113889287334783E$+$02

Table 9. Table 9: Parameter Ranges Displayed in Figures 18 - 20

Orbit	Parameter	Lower	Upper
Binary	Period (days)	7.44820	7.44838
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$-29.30541$	$-29.30515$
	$e\cos \omega$	$-0.02000$	$-0.01979$
	$e\sin \omega$	$-0.028$	$-0.012$
	Inclination (deg)	89.5	89.8
	$M_1$ ($M_{\odot }$)	0.8950	0.9905
	$M_2$ ($M_{\odot }$)	0.326	0.352
Inner	Period (days)	49.440	49.495
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$-31.375$	$-31.337$
	$e\cos \omega$	$0.006$	$0.024$
	$e\sin \omega$	$0.0122$	$0.0210$
	Inclination (deg)	89.55	90.05
	Nodal angle (deg)	$-0.32$	0.12
	$M_3$ ($M_{\oplus }$)	0.0	100.0
Middle	Period (days)	187.2	187.5
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$45.25$	$45.60$
	$e\cos \omega$	$-0.15$	$0.14$
	$e\sin \omega$	$-0.08$	$0.09$
	Inclination (deg)	90.34	90.45
	Nodal angle (deg)	$-0.4$	$-0.4$
	$M_3$ ($M_{\oplus }$)	0.0	100.0
Outer	Period (days)	303.2	303.37
	$T_{\mathrm{conj}}$ (BJD 2,455,000+)	$-55.78$	$-55.49$
	$e\cos \omega$	$-0.10$	$0.13$
	$e\sin \omega$	$-0.110$	$0.055$
	Inclination (deg)	90.17	90.21
	Nodal angle (deg)	$-2.4$	$-0.4$
	$M_3$ ($M_{\oplus }$)	0.0	12.0

Table 10. Table 10: Times, durations, and impact parameters of future planet transits

BJD–2,455,000	Year	Month	Day	UT	Impact	Duration
					Parameter	(hr)
Inner Planet
$3522.04361\pm 0.07701$	2019	Feb	7	13:02:47.6	$0.258\pm 0.033$	$4.86\pm 0.13$
$3569.12508\pm 0.06896$	2019	Mar	26	15:00:06.6	$0.269\pm 0.034$	$4.12\pm 0.07$
$3618.07766\pm 0.10808$	2019	May	14	13:51:49.8	$0.245\pm 0.027$	$6.66\pm 0.30$
$3665.39077\pm 0.06593$	2019	Jun	30	21:22:42.5	$0.286\pm 0.030$	$3.81\pm 0.04$
$3713.68886\pm 0.18417$	2019	Aug	18	04:31:57.6	$0.233\pm 0.023$	$10.88\pm 0.33$
$3761.67854\pm 0.06705$	2019	Oct	5	04:17:05.5	$0.303\pm 0.025$	$3.77\pm 0.03$
$3809.11491\pm 0.14101$	2019	Nov	21	14:45:28.6	$0.252\pm 0.023$	$7.74\pm 0.54$
$3857.95365\pm 0.07213$	2020	Jan	9	10:53:15.4	$0.317\pm 0.022$	$3.98\pm 0.06$
$3905.07251\pm 0.09442$	2020	Feb	25	13:44:24.5	$0.291\pm 0.021$	$5.10\pm 0.18$
$3954.17620\pm 0.08347$	2020	Apr	14	16:13:43.6	$0.329\pm 0.023$	$4.54\pm 0.13$
$4001.25330\pm 0.07808$	2020	May	31	18:04:45.0	$0.330\pm 0.021$	$4.13\pm 0.07$
$4050.26982\pm 0.10916$	2020	Jul	19	18:28:32.1	$0.335\pm 0.024$	$5.92\pm 0.30$
$4097.51864\pm 0.07189$	2020	Sep	5	00:26:50.1	$0.367\pm 0.022$	$3.72\pm 0.03$
$4146.02872\pm 0.17570$	2020	Oct	23	12:41:21.0	$0.333\pm 0.024$	$9.51\pm 0.68$
$4193.81227\pm 0.07070$	2020	Dec	10	07:29:40.1	$0.402\pm 0.023$	$3.60\pm 0.04$
$4241.40816\pm 0.16954$	2021	Jan	26	21:47:44.6	$0.349\pm 0.024$	$8.56\pm 0.63$
$4290.10058\pm 0.07352$	2021	Mar	16	14:24:49.9	$0.431\pm 0.022$	$3.70\pm 0.07$
$4337.27273\pm 0.10477$	2021	May	2	18:32:43.8	$0.395\pm 0.022$	$5.28\pm 0.23$
$4386.34736\pm 0.08170$	2021	Jun	20	20:20:11.5	$0.453\pm 0.024$	$4.09\pm 0.11$
$4433.42044\pm 0.08192$	2021	Aug	6	22:05:26.2	$0.445\pm 0.021$	$4.07\pm 0.08$
$4482.49229\pm 0.10084$	2021	Sep	24	23:48:53.6	$0.466\pm 0.023$	$5.07\pm 0.22$
$4529.67276\pm 0.07279$	2021	Nov	11	04:08:46.6	$0.489\pm 0.020$	$3.54\pm 0.03$
$4578.38021\pm 0.15083$	2021	Dec	29	21:07:30.2	$0.467\pm 0.023$	$7.72\pm 0.56$
$4625.96328\pm 0.06974$	2022	Feb	15	11:07:07.1	$0.527\pm 0.017$	$3.34\pm 0.04$
$4673.77625\pm 0.19135$	2022	Apr	4	06:37:48.3	$0.473\pm 0.018$	$9.29\pm 0.45$
$4722.25528\pm 0.07085$	2022	May	22	18:07:36.1	$0.555\pm 0.013$	$3.35\pm 0.05$
$4769.49523\pm 0.11299$	2022	Jul	8	23:53:07.8	$0.515\pm 0.013$	$5.45\pm 0.29$
$4818.51498\pm 0.07680$	2022	Aug	27	00:21:34.2	$0.571\pm 0.010$	$3.63\pm 0.07$
$4865.58850\pm 0.08351$	2022	Oct	13	02:07:26.6	$0.559\pm 0.010$	$3.95\pm 0.10$
$4914.69024\pm 0.09165$	2022	Dec	1	04:33:56.9	$0.576\pm 0.009$	$4.38\pm 0.15$
$4961.81627\pm 0.07225$	2023	Jan	17	07:35:25.7	$0.595\pm 0.008$	$3.34\pm 0.04$
$5010.65834\pm 0.13025$	2023	Mar	7	03:48:00.2	$0.569\pm 0.007$	$6.37\pm 0.34$
$5058.09452\pm 0.06831$	2023	Apr	23	14:16:06.2	$0.620\pm 0.009$	$3.10\pm 0.03$
$5106.12188\pm 0.19939$	2023	Jun	10	14:55:30.0	$0.562\pm 0.008$	$9.46\pm 0.27$
$5154.38086\pm 0.06892$	2023	Jul	28	21:08:26.2	$0.632\pm 0.011$	$3.10\pm 0.03$
$5201.68566\pm 0.12314$	2023	Sep	14	04:27:20.8	$0.589\pm 0.012$	$5.65\pm 0.40$
$5250.64114\pm 0.07437$	2023	Nov	2	03:23:14.8	$0.630\pm 0.021$	$3.37\pm 0.06$
$5297.71619\pm 0.08725$	2023	Dec	19	05:11:18.8	$0.618\pm 0.018$	$3.92\pm 0.17$
$5346.82635\pm 0.08815$	2024	Feb	6	07:49:56.9	$0.615\pm 0.028$	$4.06\pm 0.10$
$5393.91612\pm 0.07470$	2024	Mar	24	09:59:12.8	$0.635\pm 0.024$	$3.28\pm 0.12$
$5442.82740\pm 0.12356$	2024	May	12	07:51:27.2	$0.590\pm 0.034$	$5.84\pm 0.29$
$5490.17822\pm 0.07079$	2024	Jun	28	16:16:37.8	$0.637\pm 0.038$	$3.07\pm 0.11$
$5538.33935\pm 0.20718$	2024	Aug	15	20:08:39.8	$0.568\pm 0.039$	$9.59\pm 0.27$
$5586.45447\pm 0.07203$	2024	Oct	2	22:54:26.1	$0.626\pm 0.046$	$3.12\pm 0.11$
$5633.79496\pm 0.13951$	2024	Nov	19	07:04:44.3	$0.581\pm 0.046$	$6.05\pm 0.54$
$5682.70773\pm 0.07857$	2025	Jan	7	04:59:07.7	$0.601\pm 0.052$	$3.44\pm 0.13$
$5729.77921\pm 0.09714$	2025	Feb	23	06:42:03.4	$0.591\pm 0.053$	$4.13\pm 0.28$
$5778.88977\pm 0.09420$	2025	Apr	13	09:21:15.7	$0.564\pm 0.057$	$4.20\pm 0.15$
$5825.95882\pm 0.08342$	2025	May	30	11:00:42.0	$0.586\pm 0.060$	$3.48\pm 0.19$
$5874.89271\pm 0.13339$	2025	Jul	18	09:25:30.1	$0.521\pm 0.061$	$6.06\pm 0.32$
$5922.20903\pm 0.07985$	2025	Sep	3	17:01:00.3	$0.567\pm 0.065$	$3.29\pm 0.16$
$5970.40676\pm 0.23356$	2025	Oct	21	21:45:44.1	$0.487\pm 0.065$	$10.24\pm 0.34$
$6018.47654\pm 0.08223$	2025	Dec	8	23:26:13.1	$0.536\pm 0.055$	$3.39\pm 0.14$
$6065.81849\pm 0.16443$	2026	Jan	25	07:38:37.9	$0.491\pm 0.057$	$6.57\pm 0.67$
$6114.72278\pm 0.09079$	2026	Mar	15	05:20:47.8	$0.496\pm 0.057$	$3.75\pm 0.15$
$6161.78773\pm 0.11442$	2026	May	1	06:54:19.5	$0.488\pm 0.059$	$4.49\pm 0.31$
$6210.89773\pm 0.11007$	2026	Jun	19	09:32:44.3	$0.448\pm 0.057$	$4.56\pm 0.15$
$6257.96094\pm 0.09876$	2026	Aug	5	11:03:44.9	$0.470\pm 0.061$	$3.79\pm 0.19$
$6306.89174\pm 0.15754$	2026	Sep	23	09:24:06.5	$0.397\pm 0.060$	$6.55\pm 0.40$
$6354.20783\pm 0.09504$	2026	Nov	9	16:59:16.2	$0.443\pm 0.060$	$3.59\pm 0.14$
$6402.38804\pm 0.27883$	2026	Dec	27	21:18:46.6	$0.360\pm 0.059$	$10.99\pm 0.29$
Middle Planet
$3590.46599\pm 0.01449$	2019	Apr	16	23:11:01.3	$0.381\pm 0.031$	$5.02\pm 0.05$
$3776.89067\pm 0.01533$	2019	Oct	20	09:22:33.7	$0.362\pm 0.032$	$5.05\pm 0.04$
$3963.31897\pm 0.01657$	2020	Apr	23	19:39:19.0	$0.343\pm 0.033$	$5.15\pm 0.04$
$4149.75504\pm 0.01850$	2020	Oct	27	06:07:15.6	$0.326\pm 0.034$	$5.32\pm 0.04$
$4336.19874\pm 0.02071$	2021	May	1	16:46:11.3	$0.309\pm 0.034$	$5.60\pm 0.04$
$4522.65635\pm 0.02359$	2021	Nov	4	03:45:08.9	$0.292\pm 0.035$	$6.04\pm 0.05$
$4709.13842\pm 0.02774$	2022	May	9	15:19:19.7	$0.275\pm 0.037$	$6.72\pm 0.09$
$4895.65032\pm 0.03401$	2022	Nov	12	03:36:27.5	$0.256\pm 0.038$	$7.84\pm 0.16$
$5082.21936\pm 0.04505$	2023	May	17	17:15:52.9	$0.232\pm 0.039$	$9.97\pm 0.35$
$5268.92328\pm 0.07347$	2023	Nov	20	10:09:31.0	$0.199\pm 0.041$	$15.74\pm 1.35$
$5456.14904\pm 0.20665$	2024	May	25	15:34:37.5	$0.119\pm 0.045$	$34.28\pm 0.83$
$5643.63396\pm 0.08954$	2024	Nov	29	03:12:54.0	$0.014\pm 0.041$	$18.61\pm 1.73$
$5830.39948\pm 0.05933$	2025	Jun	3	21:35:14.8	$0.079\pm 0.044$	$11.14\pm 0.49$
$6016.99917\pm 0.04807$	2025	Dec	7	11:58:48.0	$0.127\pm 0.046$	$8.58\pm 0.24$
$6203.52776\pm 0.04248$	2026	Jun	12	00:39:58.4	$0.165\pm 0.047$	$7.26\pm 0.15$
$6390.01925\pm 0.03950$	2026	Dec	15	12:27:43.1	$0.197\pm 0.048$	$6.46\pm 0.11$
Outer Planet
$3575.38963\pm 0.03743$	2019	Apr	1	21:21:04.0	$0.285\pm 0.013$	$05.78\pm 0.04$
$3878.83019\pm 0.11337$	2020	Jan	30	07:55:28.7	$0.213\pm 0.018$	$15.62\pm 2.02$
$4179.70378\pm 0.06532$	2020	Nov	26	04:53:26.9	$0.306\pm 0.011$	$06.10\pm 0.07$
$4483.55306\pm 0.05994$	2021	Sep	26	01:16:24.1	$0.148\pm 0.014$	$06.61\pm 0.12$
$4784.28095\pm 0.14007$	2022	Jul	23	18:44:34.3	$0.316\pm 0.012$	$12.09\pm 1.76$
$5087.86638\pm 0.07416$	2023	May	23	08:47:35.0	$0.156\pm 0.014$	$05.88\pm 0.01$
$5391.64281\pm 0.10567$	2024	Mar	22	03:25:38.7	$0.027\pm 0.017$	$08.53\pm 0.44$
$5692.15323\pm 0.09213$	2025	Jan	16	15:40:39.0	$0.175\pm 0.013$	$07.24\pm 0.24$
$5996.00286\pm 0.07245$	2025	Nov	16	12:04:06.7	$0.006\pm 0.012$	$06.10\pm 0.09$
$6299.20901\pm 0.26220$	2026	Sep	15	17:00:58.4	$0.038\pm 0.040$	$33.5\pm 21.7$

Table 11. Table 11: Planetary mass and radius data

Name	Mass ($M_{\oplus }$)	Radius ($R_{\oplus }$)	Insolation ($S_{\oplus }$)	References
55Cnc e	${7.81}_{-0.53}^{+0.58}$	${2.17}_{-0.10}^{+0.10}$	${2442.2}_{-54.2}^{+54.2}$	13, 39
GJ1214 b	${6.55}_{-0.98}^{+0.98}$	${2.68}_{-0.13}^{+0.13}$	${16.3}_{-2.4}^{+2.4}$	7
GJ3470 b	${14.00}_{-1.80}^{+1.80}$	${4.20}_{-0.60}^{+0.60}$	${29.7}_{-6.8}^{+6.8}$	4
GJ436 b	${24.79}_{-2.54}^{+2.22}$	${4.05}_{-0.17}^{+0.17}$	${28.4}_{-2.0}^{+3.7}$	36
GJ1132 b	${1.63}_{-0.45}^{+0.54}$	${1.44}_{-0.21}^{+0.21}$	${10.5}_{-2.2}^{+2.2}$	40
GJ9827 K2-133	${3.74}_{-0.48}^{+0.50}$	${1.62}_{-0.16}^{+0.17}$	${311.6}_{-50.6}^{+50.6}$	41
GJ9827 K2-13	${1.47}_{-0.58}^{+0.59}$	${1.27}_{-0.13}^{+0.13}$	${71.5}_{-11.6}^{+11.6}$	41
GJ9827 K2-13	${2.38}_{-0.69}^{+0.71}$	${2.09}_{-0.21}^{+0.22}$	${35.2}_{-5.7}^{+5.7}$	41
CoRoT-7 b	${5.53}_{-0.78}^{+0.86}$	${1.58}_{-0.10}^{+0.10}$	${1572.7}_{-121.9}^{+121.9}$	16
HAT-P-11 b	${26.70}_{-2.22}^{+2.22}$	${4.36}_{-0.02}^{+0.10}$	${83.5}_{-2.1}^{+3.7}$	42
HAT-P-26 b	${18.75}_{-2.22}^{+2.22}$	${6.20}_{-0.35}^{+0.79}$	${162.9}_{-13.8}^{+30.1}$	15
HD97658b	${7.55}_{-0.79}^{+0.83}$	${2.25}_{-0.09}^{+0.10}$	${54.4}_{-3.9}^{+4.0}$	35
HD219134b	${4.46}_{-0.47}^{+0.47}$	${1.61}_{-0.09}^{+0.09}$	${181.0}_{-4.3}^{+4.3}$	28
Kepler-4 b	${24.50}_{-3.80}^{+3.80}$	${4.50}_{-0.12}^{+0.12}$	${1239.5}_{-44.0}^{+40.1}$	33
Kepler-10 b	${3.72}_{-0.42}^{+0.42}$	${1.47}_{-0.02}^{+0.03}$	${3574.9}_{-227.0}^{+227.0}$	2
Kepler-10 c	${13.98}_{-1.79}^{+1.79}$	${2.35}_{-0.04}^{+0.09}$	${17.5}_{-1.1}^{+1.1}$	10, 38
Kepler-11 b	${1.90}_{-1.00}^{+1.40}$	${1.80}_{-0.05}^{+0.03}$	${125.3}_{-5.4}^{+4.7}$	22
Kepler-11 c	${2.90}_{-1.60}^{+2.90}$	${2.87}_{-0.06}^{+0.05}$	${91.6}_{-3.9}^{+3.5}$	22
Kepler-11 d	${7.30}_{-1.50}^{+0.80}$	${3.12}_{-0.07}^{+0.06}$	${43.7}_{-1.9}^{+1.6}$	22
Kepler-11 e	${8.00}_{-2.10}^{+1.50}$	${4.19}_{-0.09}^{+0.07}$	${27.6}_{-1.2}^{+1.0}$	22
Kepler-11 f	${2.00}_{-0.90}^{+0.80}$	${2.49}_{-0.07}^{+0.04}$	${16.69}_{-0.72}^{+0.63}$	22
Kepler-18 b	${6.90}_{-3.40}^{+3.40}$	${2.00}_{-0.10}^{+0.10}$	${450.1}_{-33.3}^{+33.3}$	8
Kepler-18 c	${17.30}_{-1.90}^{+1.90}$	${5.49}_{-0.26}^{+0.26}$	${160.1}_{-11.9}^{+11.9}$	8
Kepler-18 d	${16.40}_{-1.40}^{+1.40}$	${6.98}_{-0.33}^{+0.33}$	${65.2}_{-4.8}^{+4.8}$	8
Kepler-20 b	${8.70}_{-2.20}^{+2.20}$	${1.91}_{-0.21}^{+0.12}$	${346.3}_{-51.8}^{+35.0}$	11, 12
Kepler-20 c	${16.10}_{-3.50}^{+3.50}$	${3.07}_{-0.31}^{+0.20}$	${82.9}_{-12.4}^{+8.4}$	11, 12
Kepler-23 b	${1.30}_{-0.50}^{+1.30}$	${1.80}_{-0.10}^{+0.10}$	${450.8}_{-36.3}^{+36.3}$	43
Kepler-23 c	${2.20}_{-0.90}^{+2.80}$	${3.20}_{-0.20}^{+0.20}$	${259.9}_{-21.0}^{+21.0}$	43
Kepler-24 b	${2.00}_{-0.70}^{+1.60}$	${2.70}_{-0.50}^{+0.60}$	${230.3}_{-54.4}^{+124.1}$	43
Kepler-24 c	${1.70}_{-0.70}^{+1.20}$	${3.20}_{-1.80}^{+1.70}$	${132.4}_{-31.3}^{+71.4}$	43
Kepler-25 b	${24.60}_{-5.70}^{+5.70}$	${5.20}_{-0.09}^{+0.09}$	${186.3}_{-9.9}^{+9.9}$	25
Kepler-25 c	${9.60}_{-4.20}^{+4.20}$	${2.71}_{-0.05}^{+0.05}$	${481.7}_{-25.5}^{+25.5}$	25
Kepler-26 b	${5.12}_{-0.61}^{+0.65}$	${2.78}_{-0.11}^{+0.11}$	${7.64}_{-0.65}^{+0.65}$	20
Kepler-26 c	${6.20}_{-0.65}^{+0.65}$	${2.72}_{-0.12}^{+0.12}$	${4.85}_{-0.41}^{+0.41}$	20
Kepler-27 b	${3.60}_{-1.60}^{+1.70}$	${4.90}_{-0.30}^{+0.40}$	${52.5}_{-5.0}^{+3.6}$	43
Kepler-27 c	${3.60}_{-1.60}^{+2.00}$	${7.00}_{-0.50}^{+0.50}$	${20.2}_{-1.9}^{+1.4}$	43
Kepler-28 b	${1.20}_{-0.50}^{+0.80}$	${1.80}_{-0.10}^{+0.10}$	${27.8}_{-7.4}^{+7.4}$	43
Kepler-28 c	${1.00}_{-0.30}^{+0.70}$	${1.70}_{-0.10}^{+0.20}$	${15.9}_{-4.2}^{+4.2}$	43
Kepler-29 b	${4.51}_{-1.47}^{+1.41}$	${3.35}_{-0.22}^{+0.22}$	${96.8}_{-10.4}^{+10.4}$	20
Kepler-29 c	${4.00}_{-1.29}^{+1.23}$	${3.14}_{-0.20}^{+0.20}$	${69.3}_{-7.4}^{+7.4}$	20
Kepler-30 b	${11.30}_{-1.40}^{+1.40}$	${3.90}_{-0.20}^{+0.20}$	${21.5}_{-4.1}^{+4.1}$	31
Kepler-30 d	${23.10}_{-2.70}^{+2.70}$	${8.80}_{-0.50}^{+0.50}$	${2.59}_{-0.50}^{+0.50}$	31
Kepler-31 b	${0.70}_{-0.60}^{+2.40}$	${5.60}_{-1.00}^{+1.00}$	${92.7}_{-41.4}^{+44.8}$	43
Kepler-31 c	${2.20}_{-1.20}^{+4.60}$	${5.30}_{-1.00}^{+0.90}$	${35.8}_{-16.0}^{+17.3}$	43
Kepler-31 d	${2.80}_{-1.40}^{+4.20}$	${4.20}_{-0.80}^{+0.80}$	${13.7}_{-6.1}^{+6.6}$	43
Kepler-32 b	${0.50}_{-0.10}^{+0.30}$	${2.00}_{-0.10}^{+0.20}$	${16.8}_{-2.8}^{+2.1}$	43
Kepler-32 c	${0.40}_{-0.10}^{+0.20}$	${1.90}_{-0.20}^{+0.20}$	${9.9}_{-1.6}^{+1.2}$	43
Kepler-33 c	${0.80}_{-0.70}^{+2.50}$	${3.20}_{-0.30}^{+0.30}$	${255.5}_{-38.3}^{+31.1}$	14, 22
Kepler-33 d	${4.70}_{-2.00}^{+2.00}$	${5.40}_{-0.50}^{+0.50}$	${130.8}_{-19.6}^{+15.9}$	14, 22
Kepler-33 e	${6.70}_{-1.30}^{+1.20}$	${4.00}_{-0.40}^{+0.40}$	${79.0}_{-11.8}^{+9.6}$	14, 22
Kepler-33 f	${11.50}_{-2.10}^{+1.80}$	${4.50}_{-0.40}^{+0.40}$	${56.2}_{-8.4}^{+6.8}$	14, 22
Kepler-36 b	${4.45}_{-0.27}^{+0.33}$	${1.49}_{-0.04}^{+0.03}$	${217.4}_{-9.3}^{+9.3}$	6
Kepler-36 c	${8.08}_{-0.46}^{+0.60}$	${3.68}_{-0.05}^{+0.05}$	${175.8}_{-7.6}^{+7.6}$	6
Kepler-47 b	${2.07}_{-2.07}^{+23.70}$	${3.05}_{-0.04}^{+0.04}$	${11.74}_{-0.43}^{+0.43}$	this work
Kepler-47 c	${19.02}_{-11.67}^{+23.84}$	${7.04}_{-0.49}^{+0.66}$	${1.989}_{-0.072}^{+0.072}$	this work
Kepler-47 d	${3.17}_{-1.25}^{+2.18}$	${4.65}_{-0.07}^{+0.09}$	${1.047}_{-0.038}^{+0.038}$	this work
Kepler-48 b	${1.90}_{-1.29}^{+5.34}$	${1.89}_{-0.23}^{+0.23}$	${171.4}_{-31.0}^{+31.0}$	25
Kepler-48 c	${2.81}_{-1.58}^{+3.78}$	${2.61}_{-0.32}^{+0.32}$	${66.9}_{-12.1}^{+12.1}$	25
Kepler-49 b	${5.09}_{-1.94}^{+2.11}$	${2.35}_{-0.09}^{+0.09}$	${14.7}_{-2.0}^{+2.0}$	20
Kepler-49 c	${3.28}_{-1.32}^{+1.45}$	${2.06}_{-0.09}^{+0.09}$	${8.5}_{-1.2}^{+1.2}$	20
Kepler-51 b	${2.10}_{-0.80}^{+1.50}$	${7.10}_{-0.30}^{+0.30}$	${16.5}_{-2.1}^{+2.1}$	26
Kepler-51 c	${4.00}_{-0.40}^{+0.40}$	${9.00}_{-1.70}^{+2.80}$	${7.04}_{-0.88}^{+0.88}$	26
Kepler-51 d	${7.60}_{-1.10}^{+1.10}$	${9.70}_{-0.50}^{+0.50}$	${4.01}_{-0.50}^{+0.50}$	26
Kepler-52 b	${2.30}_{-0.80}^{+1.50}$	${2.30}_{-0.10}^{+0.20}$	${23.0}_{-2.2}^{+1.9}$	43
Kepler-52 c	${6.30}_{-2.00}^{+5.00}$	${2.10}_{-0.10}^{+0.10}$	${8.67}_{-0.82}^{+0.73}$	43
Kepler-53 c	${2.80}_{-0.90}^{+1.90}$	${3.40}_{-0.40}^{+0.40}$	${57.0}_{-31.8}^{+31.8}$	43
Kepler-53 d	${3.50}_{-1.10}^{+2.30}$	${3.60}_{-0.40}^{+0.50}$	${21.7}_{-12.1}^{+12.1}$	43
Kepler-54 b	${0.80}_{-0.30}^{+1.00}$	${2.30}_{-0.30}^{+0.40}$	${10.0}_{-2.1}^{+2.1}$	43
Kepler-54 c	${0.70}_{-0.30}^{+0.90}$	${1.30}_{-0.10}^{+0.10}$	${5.8}_{-1.2}^{+1.2}$	43
Kepler-55 b	${1.30}_{-0.60}^{+2.10}$	${1.90}_{-0.10}^{+0.20}$	${5.99}_{-0.45}^{+0.65}$	43
Kepler-55 c	${1.80}_{-0.80}^{+2.80}$	${2.50}_{-0.30}^{+0.30}$	${3.47}_{-0.26}^{+0.38}$	43
Kepler-57 b	${16.85}_{-6.95}^{+7.24}$	${1.88}_{-0.14}^{+0.14}$	${73.7}_{-9.0}^{+9.0}$	20
Kepler-57 c	${4.12}_{-1.89}^{+1.97}$	${1.39}_{-0.10}^{+0.10}$	${28.9}_{-3.5}^{+3.5}$	20
Kepler-58 b	${2.90}_{-1.10}^{+2.10}$	${2.80}_{-0.60}^{+0.50}$	${171.6}_{-38.3}^{+84.0}$	43
Kepler-58 c	${2.30}_{-0.90}^{+1.80}$	${2.70}_{-0.50}^{+0.60}$	${97.8}_{-21.9}^{+47.9}$	43
Kepler-60 b	${4.18}_{-0.52}^{+0.56}$	${1.71}_{-0.13}^{+0.13}$	${319.1}_{-41.9}^{+41.9}$	20
Kepler-60 c	${3.85}_{-0.81}^{+0.81}$	${1.90}_{-0.15}^{+0.15}$	${236.9}_{-31.1}^{+31.1}$	20
Kepler-60 d	${4.16}_{-0.52}^{+0.56}$	${1.99}_{-0.16}^{+0.16}$	${161.3}_{-21.2}^{+21.2}$	20
Kepler-68 b	${5.97}_{-1.70}^{+1.70}$	${2.33}_{-0.02}^{+0.02}$	${406.9}_{-20.8}^{+20.8}$	25
Kepler-78 b	${1.69}_{-0.41}^{+0.41}$	${1.20}_{-0.09}^{+0.09}$	${3973.8}_{-432.5}^{+432.5}$	17, 30
Kepler-79 b	${10.90}_{-6.00}^{+7.40}$	${3.47}_{-0.07}^{+0.07}$	${162.5}_{-8.8}^{+8.8}$	18
Kepler-79 c	${5.90}_{-2.30}^{+1.90}$	${3.72}_{-0.08}^{+0.08}$	${63.1}_{-3.4}^{+3.4}$	18
Kepler-79 d	${6.00}_{-1.60}^{+2.10}$	${7.16}_{-0.16}^{+0.13}$	${26.8}_{-1.4}^{+1.4}$	18
Kepler-79 e	${4.10}_{-1.10}^{+1.20}$	${3.49}_{-0.14}^{+0.14}$	${14.85}_{-0.80}^{+0.80}$	18
Kepler-80 d	${6.75}_{-0.51}^{+0.69}$	${1.53}_{-0.07}^{+0.09}$	${126.3}_{-8.8}^{+8.8}$	24
Kepler-80 e	${4.92}_{-0.37}^{+0.49}$	${1.60}_{-0.07}^{+0.08}$	${72.8}_{-5.1}^{+5.1}$	24
Kepler-80 b	${5.99}_{-0.57}^{+0.49}$	${2.67}_{-0.10}^{+0.10}$	${41.7}_{-2.9}^{+2.9}$	24
Kepler-80 c	${5.03}_{-0.42}^{+0.40}$	${2.74}_{-0.10}^{+0.12}$	${28.0}_{-2.0}^{+2.0}$	24
Kepler-81 b	${0.20}_{-0.10}^{+1.70}$	${2.40}_{-0.10}^{+0.10}$	${40.3}_{-3.5}^{+3.6}$	43
Kepler-81 c	${0.10}_{-0.10}^{+1.10}$	${2.30}_{-0.10}^{+0.10}$	${15.8}_{-1.4}^{+1.4}$	43
Kepler-84 b	${5.00}_{-2.90}^{+5.70}$	${2.50}_{-0.40}^{+0.40}$	${221.4}_{-48.3}^{+52.7}$	43
Kepler-85 b	${0.70}_{-0.40}^{+1.90}$	${1.80}_{-0.20}^{+0.20}$	${99.4}_{-12.5}^{+17.5}$	43
Kepler-85 c	${1.00}_{-0.60}^{+2.20}$	${1.90}_{-0.20}^{+0.20}$	${57.5}_{-7.2}^{+10.1}$	43
Kepler-85 e	${0.60}_{-0.40}^{+0.50}$	${1.20}_{-0.10}^{+0.10}$	${22.6}_{-2.8}^{+4.0}$	43
Kepler-87 c	${6.40}_{-0.80}^{+0.80}$	${6.14}_{-0.29}^{+0.29}$	${6.51}_{-0.34}^{+0.34}$	29
Kepler-89 b	${10.50}_{-4.60}^{+4.60}$	${1.71}_{-0.16}^{+0.16}$	${1194.7}_{-161.8}^{+161.8}$	37
Kepler-89 c	${9.40}_{-2.10}^{+2.40}$	${4.32}_{-0.41}^{+0.41}$	${295.9}_{-40.1}^{+40.1}$	26, 37
Kepler-89 e	${13.00}_{-2.10}^{+2.50}$	${6.56}_{-0.62}^{+0.62}$	${32.8}_{-4.4}^{+4.4}$	26, 37
Kepler-93 b	${4.02}_{-0.48}^{+0.68}$	${1.48}_{-0.02}^{+0.02}$	${274.1}_{-11.8}^{+11.8}$	9
Kepler-94 b	${10.84}_{-1.40}^{+1.40}$	${3.51}_{-0.15}^{+0.15}$	${238.7}_{-22.0}^{+22.0}$	25
Kepler-95 b	${13.00}_{-2.90}^{+2.90}$	${3.42}_{-0.09}^{+0.09}$	${179.3}_{-13.8}^{+13.8}$	25
Kepler-98 b	${3.55}_{-1.60}^{+1.60}$	${1.99}_{-0.22}^{+0.22}$	${1530.0}_{-249.1}^{+249.1}$	25
Kepler-99 b	${6.15}_{-1.30}^{+1.30}$	${1.48}_{-0.08}^{+0.08}$	${99.7}_{-11.3}^{+11.3}$	25
Kepler-102 e	${8.93}_{-2.00}^{+2.00}$	${2.22}_{-0.07}^{+0.07}$	${21.1}_{-1.6}^{+1.6}$	25
KOI-115.03	${1.13}_{-0.83}^{+1.21}$	${0.73}_{-0.04}^{+0.04}$	${423.0}_{-36.4}^{+36.4}$	20
Kepler-105 b	${3.74}_{-1.99}^{+1.99}$	${2.22}_{-0.11}^{+0.11}$	${233.3}_{-20.1}^{+20.1}$	20
Kepler-105 c	${4.60}_{-0.88}^{+0.92}$	${1.31}_{-0.07}^{+0.07}$	${161.6}_{-13.9}^{+13.9}$	20
Kepler-106 c	${10.44}_{-3.20}^{+3.20}$	${2.50}_{-0.32}^{+0.32}$	${92.2}_{-22.1}^{+22.1}$	25
Kepler-114 c	${0.80}_{-0.30}^{+0.60}$	${1.00}_{-0.10}^{+0.10}$	${15.0}_{-1.6}^{+1.5}$	43
Kepler-122 f	${2.10}_{-1.00}^{+2.50}$	${2.30}_{-0.40}^{+0.30}$	${20.5}_{-5.0}^{+5.0}$	43
Kepler-127 b	${0.80}_{-0.40}^{+2.30}$	${1.40}_{-0.04}^{+0.10}$	${146.9}_{-7.7}^{+12.9}$	43
Kepler-127 c	${1.00}_{-0.60}^{+3.50}$	${2.40}_{-0.10}^{+0.10}$	${56.9}_{-3.0}^{+5.0}$	43
Kepler-127 d	${6.10}_{-3.50}^{+3.10}$	${2.60}_{-0.10}^{+0.10}$	${29.1}_{-1.5}^{+2.5}$	43
Kepler-128 b	${1.15}_{-0.52}^{+1.92}$	${1.43}_{-0.03}^{+0.03}$	${210.0}_{-13.7}^{+13.7}$	44
Kepler-128 c	${1.30}_{-0.59}^{+2.25}$	${1.42}_{-0.03}^{+0.03}$	${121.1}_{-7.9}^{+7.9}$	44
Kepler-131 b	${16.13}_{-3.50}^{+3.50}$	${2.41}_{-0.20}^{+0.20}$	${63.1}_{-9.3}^{+9.3}$	25
Kepler-138 b	${0.07}_{-0.04}^{+0.06}$	${0.52}_{-0.03}^{+0.03}$	${11.2}_{-1.6}^{+1.6}$	19
Kepler-138 c	${1.97}_{-1.12}^{+1.91}$	${1.20}_{-0.07}^{+0.07}$	${7.6}_{-1.1}^{+1.1}$	19
Kepler-138 d	${0.64}_{-0.39}^{+0.67}$	${1.21}_{-0.08}^{+0.08}$	${3.82}_{-0.56}^{+0.56}$	19
Kepler-176 c	${0.10}_{-0.04}^{+0.30}$	${2.50}_{-0.20}^{+0.10}$	${37.7}_{-4.2}^{+4.5}$	43
Kepler-176 d	${0.20}_{-0.10}^{+0.50}$	${2.30}_{-0.20}^{+0.10}$	${14.8}_{-1.6}^{+1.8}$	43
Kepler-223 b	${7.40}_{-1.10}^{+1.30}$	${2.99}_{-0.27}^{+0.18}$	${653.7}_{-91.4}^{+70.0}$	27
Kepler-223 c	${5.10}_{-1.10}^{+1.70}$	${3.44}_{-0.30}^{+0.20}$	${445.4}_{-62.3}^{+47.7}$	27
Kepler-223 d	${8.00}_{-1.30}^{+1.50}$	${5.24}_{-0.45}^{+0.26}$	${258.9}_{-36.2}^{+27.7}$	27
Kepler-223 e	${4.80}_{-1.20}^{+1.40}$	${4.60}_{-0.41}^{+0.27}$	${176.4}_{-24.7}^{+18.9}$	27
Kepler-231 b	${4.90}_{-1.30}^{+1.80}$	${1.82}_{-0.25}^{+0.26}$	${7.6}_{-1.5}^{+1.6}$	21
Kepler-231 c	${2.20}_{-1.10}^{+1.50}$	${1.69}_{-0.23}^{+0.24}$	${3.35}_{-0.68}^{+0.72}$	21
Kepler-238 f	${4.50}_{-2.20}^{+5.40}$	${2.80}_{-0.50}^{+0.50}$	${24.5}_{-7.0}^{+6.4}$	43
Kepler-277 b	${3.50}_{-1.60}^{+2.80}$	${2.90}_{-0.70}^{+0.50}$	${202.4}_{-71.4}^{+43.7}$	43
Kepler-277 c	${4.50}_{-2.00}^{+3.70}$	${3.00}_{-0.70}^{+0.50}$	${85.7}_{-30.2}^{+18.5}$	43
Kepler-279 c	${7.40}_{-1.50}^{+2.20}$	${5.50}_{-1.00}^{+1.10}$	${24.34}_{-13.21}^{+0.84}$	43
Kepler-279 d	${4.50}_{-0.90}^{+1.20}$	${4.30}_{-0.90}^{+1.00}$	${13.89}_{-7.54}^{+0.48}$	43
Kepler-289 b	${7.30}_{-6.80}^{+6.80}$	${2.15}_{-0.10}^{+0.10}$	${25.47}_{-0.88}^{+0.88}$	32
Kepler-289 c	${4.00}_{-0.90}^{+0.90}$	${2.68}_{-0.17}^{+0.17}$	${10.73}_{-0.37}^{+0.37}$	32
Kepler-305 b	${2.70}_{-1.10}^{+2.90}$	${3.10}_{-0.20}^{+0.20}$	${118.0}_{-11.2}^{+12.7}$	43
Kepler-305 c	${1.40}_{-0.70}^{+1.20}$	${2.70}_{-0.20}^{+0.20}$	${68.1}_{-6.5}^{+7.3}$	43
Kepler-305 d	${9.10}_{-3.80}^{+6.10}$	${2.70}_{-0.20}^{+0.20}$	${26.7}_{-2.5}^{+2.9}$	43
Kepler-307 b	${7.44}_{-0.87}^{+0.91}$	${2.35}_{-0.09}^{+0.09}$	${60.4}_{-3.8}^{+3.8}$	20
Kepler-307 c	${3.64}_{-0.58}^{+0.65}$	${2.20}_{-0.07}^{+0.07}$	${44.6}_{-2.8}^{+2.8}$	20
Kepler-310 d	${7.00}_{-4.10}^{+3.40}$	${2.10}_{-0.10}^{+0.20}$	${3.58}_{-0.38}^{+0.41}$	43
KOI1831.03	${1.00}_{-0.40}^{+1.20}$	${1.20}_{-0.10}^{+0.10}$	${12.1}_{-1.2}^{+1.5}$	43
Kepler-345 b	${0.50}_{-0.30}^{+0.30}$	${0.80}_{-0.10}^{+0.10}$	${35.2}_{-2.7}^{+4.0}$	43
Kepler-345 c	${2.20}_{-0.90}^{+0.90}$	${1.30}_{-0.10}^{+0.10}$	${25.7}_{-1.9}^{+2.9}$	43
Kepler-359 c	${2.90}_{-1.90}^{+2.40}$	${4.80}_{-0.90}^{+1.00}$	${16.3}_{-3.5}^{+7.4}$	43
Kepler-359 d	${2.70}_{-1.50}^{+2.50}$	${4.60}_{-0.90}^{+0.90}$	${11.1}_{-2.4}^{+5.0}$	43
Kepler-396 b	${1.10}_{-0.20}^{+0.50}$	${3.30}_{-0.40}^{+0.30}$	${14.6}_{-2.2}^{+3.0}$	43
Kepler-396 c	${1.10}_{-0.20}^{+0.60}$	${5.60}_{-0.60}^{+0.50}$	${5.6}_{-0.8}^{+1.2}$	43
Kepler-406 b	${6.35}_{-1.40}^{+1.40}$	${1.43}_{-0.03}^{+0.03}$	${739.6}_{-43.9}^{+43.9}$	25
Kepler-444 d	${0.04}_{-0.02}^{+0.06}$	${0.54}_{-0.02}^{+0.02}$	${89.2}_{-7.9}^{+8.0}$	45
Kepler-444 e	${0.03}_{-0.02}^{+0.06}$	${0.56}_{-0.02}^{+0.02}$	${66.2}_{-5.8}^{+5.9}$	45
Kepler-526 c	${0.50}_{-0.30}^{+0.40}$	${0.80}_{-0.10}^{+0.10}$	${250.4}_{-31.0}^{+57.9}$	43
Kepler-549 b	${11.00}_{-3.20}^{+4.20}$	${3.10}_{-0.40}^{+0.40}$	${12.9}_{-2.2}^{+2.3}$	43
K2-2 b	${11.82}_{-1.33}^{+1.33}$	${2.53}_{-0.18}^{+0.18}$	${50.1}_{-2.9}^{+2.9}$	46
K2-3 b	${6.60}_{-1.10}^{+1.10}$	${2.29}_{-0.23}^{+0.23}$	${11.6}_{-1.9}^{+1.9}$	47
K2-3 c	${3.10}_{-1.20}^{+1.30}$	${1.77}_{-0.18}^{+0.18}$	${3.50}_{-0.58}^{+0.58}$	47
K2-3 d	${1.20}_{-0.70}^{+1.40}$	${1.65}_{-0.17}^{+0.17}$	${1.59}_{-0.26}^{+0.26}$	47
K2-18 b	${7.96}_{-1.91}^{+1.91}$	${2.24}_{-0.23}^{+0.23}$	${1.06}_{-0.18}^{+0.18}$	48
K2-19 c	${15.90}_{-2.80}^{+7.70}$	${4.51}_{-0.47}^{+0.47}$	${62.7}_{-10.9}^{+10.9}$	49
K2-38 b	${12.00}_{-2.90}^{+2.90}$	${1.55}_{-0.16}^{+0.16}$	${466.2}_{-57.6}^{+57.6}$	50
K2-38 c	${9.90}_{-4.60}^{+4.60}$	${2.42}_{-0.29}^{+0.29}$	${128.4}_{-15.9}^{+15.9}$	50
K2-66 b	${21.30}_{-3.60}^{+3.60}$	${2.49}_{-0.34}^{+0.34}$	${840.4}_{-89.9}^{+89.9}$	51
K2-110 b	${16.70}_{-3.20}^{+3.20}$	${2.60}_{-0.10}^{+0.10}$	${27.6}_{-1.3}^{+1.3}$	52
K2-106 b	${8.36}_{-0.94}^{+0.96}$	${1.52}_{-0.16}^{+0.16}$	${3469.2}_{-535.4}^{+535.4}$	51
K2-106 c	${5.80}_{-3.00}^{+3.30}$	${2.50}_{-0.26}^{+0.27}$	${52.0}_{-8.0}^{+8.0}$	51
K2-131 b	${6.80}_{-1.60}^{+1.60}$	${1.81}_{-0.12}^{+0.16}$	${4763.9}_{-339.3}^{+339.3}$	53
K2-111 b	${8.60}_{-3.90}^{+3.90}$	${1.90}_{-0.20}^{+0.20}$	${528.8}_{-81.2}^{+81.2}$	54
K2-106 b	${9.00}_{-1.60}^{+1.60}$	${1.82}_{-0.14}^{+0.20}$	${4301.3}_{-347.6}^{+347.6}$	51
LHS1140 b	${6.65}_{-1.82}^{+1.82}$	${1.43}_{-0.10}^{+0.10}$	${0.390}_{-0.062}^{+0.062}$	54
Trappist1 b	${1.02}_{-0.15}^{+0.14}$	${1.12}_{-0.03}^{+0.03}$	${4.26}_{-0.40}^{+0.40}$	56
Trappist1 c	${1.16}_{-0.14}^{+0.13}$	${1.10}_{-0.03}^{+0.03}$	${2.27}_{-0.22}^{+0.22}$	56
Trappist1 d	${0.30}_{-0.04}^{+0.04}$	${0.78}_{-0.02}^{+0.02}$	${1.15}_{-0.11}^{+0.11}$	56
Trappist1 e	${0.77}_{-0.08}^{+0.08}$	${0.91}_{-0.03}^{+0.03}$	${0.663}_{-0.063}^{+0.063}$	56
Trappist1 f	${0.93}_{-0.08}^{+0.08}$	${1.05}_{-0.03}^{+0.03}$	${0.383}_{-0.036}^{+0.036}$	56
Trappist1 g	${0.33}_{-0.06}^{+0.05}$	${0.77}_{-0.03}^{+0.03}$	${0.383}_{-0.036}^{+0.036}$	56
Wasp-47 d	${15.20}_{-7.00}^{+7.00}$	${3.63}_{-0.14}^{+0.14}$	${155.1}_{-12.9}^{+12.9}$	3
K2-229 b	${2.59}_{-0.430}^{+0.430}$	${1.164}_{-0.048}^{+0.066}$	${2454.3}_{-110.3}^{+150.4}$	57

Equations

$$I(\mu )/I_0=1-u_1(1-\mu )-u_2(1-\mu {)}^2$$

$$\frac{S}{S_{\oplus }}={\left(\frac{T_{\mathrm{eff}}}{T_{\odot }}\right)}^4{\left(\frac{1\mathrm{AU}}{a}\right)}^2{\left(\frac{R_{\odot }}{R_{\ast }}\right)}^2$$

Full text

planets and satellites: detection – (stars:) binaries: eclipsing – methods: data analysis – techniques: photometric – methods: observational – techniques: spectroscopic –

1 Introduction

Binary stars are common: for example, roughly half of all Sun-like stars are found in pairs (Raghavan et al., 2010). Circumbinary planets (hereafter, CBPs), i.e., planets that orbit an entire binary star system, can reveal their presence in a variety of different ways. Known as P-type systems (Dvorak, 1982), these planets, if sufficiently large, may affect radial velocities of the stars or the timing of stellar eclipses (Schwarz et al., 2011). Several claims of the detection of CBPs using eclipse timing variations have been made (for example Qian et al., 2012a, b, 2013; Lee et al., 2014), but, the validity of these detections is currently under question (e.g., Marsh, 2018). A CBP may also transit one or both stars in the binary. When the system is viewed edge-on (i.e., as in an eclipsing binary), transits in the light curve provide the most secure detection of the planet. To date, ten transiting CBPs have been discovered with this method (Doyle et al., 2011; Orosz et al., 2012a, b; Kostov et al., 2013, 2014, 2016; Welsh et al., 2012, 2015; Schwamb et al., 2013) in the set of $\sim$2900 eclipsing binary stars111http://keplerebs.villanova.edu (Prša et al., 2011; Slawson et al., 2011; Kirk et al., 2016) observed with the Kepler telescope (Borucki et al., 2010a; Koch et al., 2010). Transits are particularly valuable because they allow the planet’s radius to be measured, and the large variations in transit durations and the deviations from periodicity unambiguously demonstrate that the body is in a circumbinary orbit (Orosz et al., 2012b).

Kepler-47, the topic of this paper, is the only known multi-planet circumbinary system, and the multi-planet nature of this system allows us to study an ensemble of star–planet and planet–planet dynamics. Its binary stars have an orbital period of $\sim 7.5$ days, the shortest known in any CBP system. The innermost planet (Kepler-47 b) has an orbital period of $\sim 49.5$ days, and is the smallest of the known CBPs ($\sim 3\,R_{\oplus}$). The Uranus-size outer planet (Kepler-47 c) has a $\sim 303$ day orbit, placing it well within the habitable zone (Haghighipour & Kaltenegger, 2013), i.e., the region surrounding a star where water could exist in a liquid state on a terrestrial planet (Kopparapu et al., 2013). The existence of the Kepler-47 system, particularly with the newly discovered middle planet described herein, shows that multi-planetary systems can form and survive around close binary stars. This despite theory suggesting that strong perturbations from the binary on the protoplanetary disk most likely inhibit in situ formation in regions close to the binary (Meschiari, 2012; Paardekooper et al., 2012; Marzari et al., 2013; Dunhill et al., 2013; Pierens & Nelson, 2013; Rafikov, 2013; Lines et al., 2015a, b, 2016).

In the system discovery paper (Orosz et al., 2012b), eighteen transits of the inner planet and three of the outer planet were detected. An additional 0.2% deep “orphan” transit event was noted, a feature that could not be attributed to the two planets known at that time. Further Kepler observations revealed two more transits not attributable to the two known planets, suggesting the existence of a third planetary body and allowing for a robust determination of its orbital period. This, in turn, allowed us to identify three much weaker transits in the earlier data. In total, six transits were observed that can be attributed to an additional planet, Kepler-47 d, orbiting between planets b and c. In this work, we provide an updated analysis of the five-body dynamics. The paper is organized as follows. We discuss the available observational material in §2. The transits of the new planet and transit times for all of the planets are discussed in §3. Measurements of times when the primary and secondary stars eclipse one another are given in §4. Our five-body photodynamical model is described in §5. The mass estimates for the planets, the long-term dynamical stability of the system, the evolutionary status of the two stars, and the locations of the middle and outer planets to the habitable zone of the binary are discussed in §6. We summarize our results in §7.

2 Observational Data

2.1 Optical Light Curves

Kepler-47 (KIC 10020423, KOI 3154, 2MASS J19411149+4655136) was observed in long-cadence mode ($\approx 30$ minute samples) from Kepler Quarter 1 (BJD 2,454,964.51 or 2009 May 13) through the end of the nominal spacecraft operation in the second month of Quarter 17 (BJD 2,456,424.00 or 2013 May 11). There were a total of 65,428 long-cadence observations during this span (including flagged data), for a duty cycle of 91.6%. Starting in Quarter 14 (BJD 2,456,107.13 or 2012 June 28), Kepler-47 was on the list of short-cadence targets ($\approx 1$ minute samples), where it remained until the cessation of normal spacecraft operations. There were a total of 409,046 short-cadence observations, including flagged data. We downloaded the FITS tables from the Mikulski Archive for Space Telescopes (MAST). These files were processed using the SOC release 9.0.3, and they include a $\approx 1$ minute correction to the time stamps that was not available at the time of original work reported in Orosz et al. (2012b).

The Kepler light curves of Kepler-47 required detrending, owing to instrumental trends and modulations due to star spots. We arrived at our final detrended light curves using a two-step process. First, an interactive technique described in Orosz et al. (2012b) and Bass et al. (2012) was used. In this step, the light curve was broken up into segments using data gaps and sudden jumps in the flux as end points. For each segment, data in the eclipses and transits were masked out, and cubic splines were fit to the remaining data. The segments were normalized by the spline fits, and pieced together to produce a detrended light curve. While effective, this technique is not entirely reproducible; it is also difficult to use with the short-cadence data, owing to limitations in our implementation of the method. Thus, an automated and reproducible detrending algorithm was devised as follows: The light curve that was detrended manually was modeled as described in §5. Once a good fit was found, the model was used to precisely determine the times of the eclipses and transits and the durations of these events. The times and durations were then used to determine which data points occurred during an eclipse or a transit event, which allowed us to assign those points zero weight during the normalization.

In our original work, only long-cadence Kepler data were available to us. As noted above, Kepler-47 was on the list of short-cadence targets for Kepler Quarters Q14 through Q17. In our view, using short-cadence data (when available) instead of long-cadence data is preferable because the ingress and egress phases of eclipses and transits are better-resolved. Thus, a “raw” light curve using long-cadence data from Quarter 1 through Quarter 13 and short-cadence data from Quarter 14 through Quarter 17 was made, and data with quality flags with values greater than 16 were eliminated (in our experience, obviously poor cadences almost always have data quality flags greater than 16). The overall duty cycle of the retained data was 91.6%. Segments centered at the mid-eclipse and transit times with widths of three times the duration of the event in question were extracted from the raw light curve. In a few cases, these segments included both an eclipse event and a transit event; also in a few cases, the segments were truncated owing to gaps in the data. A fifth-order Legendre polynomial was fit to each segment, where data that occurred during an eclipse or a transit event were given zero weight. The normalized segments were assembled to give the final normalized and trimmed light curve. We also produced a normalized light curve consisting entirely of long-cadence observations for the purposes of plotting the relatively weak transits.

2.2 Radial Velocities

Owing to the relative faintness of the secondary star (it is about 0.5% as bright as the primary in the optical), Kepler-47 is a single-lined binary, so only radial velocity measurements for the primary star are available. We used 11 radial velocity measurements of the primary in the original work (Orosz et al., 2012b). Four of these came from the HRS spectrograph on the Hobby–Eberly Telescope (HET), six of these came from the Tull Coudé spectrograph on the Harlan J. Smith 2.7m telescope (HJST) at McDonald Observatory, and one measurement came from the HIRES spectrograph on the Keck I telescope. Since that time, Kostov et al. (2013) have provided an additional four measurements with the SOPHIE spectrograph on the 1.93m telescope at Haute-Provence Observatory. These measurements are summarized in Table 1 and displayed in Figure 1. Small differences were found in the velocity zero-points between the different observatories, so these were fitted for and removed. To arrive at our final radial velocity curve, we subtracted systemic offsets of 4.5680 km s*-1* for the HJST radial velocities, 4.5999 km s*-1* for the HET radial velocities, and 4.2685 km s*-1* for the Kostov et al. (2013) velocities. Given that we have only one measurement from Keck, we cannot determine the systematic offset and we do not use that measurement in our analysis. The uncertainties in the individual measurements in each set were scaled to give $\chi^{2}=N$, where $N=14$ is the number of radial velocity measurements. The offset velocities with the scaled uncertainties are also given in Table 1.

3 New Transits and Transit Times

Figure 2 shows the transits of the middle planet across the primary star. During the nominal Kepler mission, this planet transited the primary eight times. However, we were able to observe only six of these transits. The first transit near day 45.15 was blended with an eclipse of the primary, and the second transit near day 231.58 was lost due to a gap in the data. As discussed further below, the primary transits of the middle planet were getting deeper with time, owing to nodal precession of its orbit. The three weak primary transits near days 418.0, 604.5, and 790.9 were not identified until a preliminary orbital solution based on the last three observed primary transits was available. According to our best-fitting model discussed below, the middle planet did not transit the secondary star during the nominal Kepler mission. Even if transits of the secondary star did occur, they would be nearly impossible to see in the data, owing to the faintness of the secondary star.

There were a total of 25 transits of the inner planet Kepler-47 b across the primary observed in the complete Kepler data set, seven of which are shown here for the first time (see Figure 3). Six primary transits of the inner planet were missed due to gaps in the data. A total of four transits of the outer planet Kepler-47c across the primary are available in the complete Kepler data set (see Figure 4), and fortunately, none were missed. According to our best-fitting model discussed below, the inner planet had seven transits across the secondary star, and the outer planet had one transit across the secondary star during the nominal Kepler mission. All of these secondary transits had impact parameters larger than 0.75 and are not seen above the noise level in the data, owing to the extreme faintness of the secondary star. Hereafter, when we use the term “transit,” it is understood to be a transit of a planet across the primary star.

We measured the transit times for each of the planets using the ELC code (Orosz & Hauschildt, 2000). Separate data files were made from the normalized light curve for each event, where between 0.5 and 1.0 days of out-of-eclipse data were kept on either side of the transit event, when possible. In cases where the transit events occurred near eclipses, the eclipses were trimmed from the data.

To model the profile of a planet transit across a single star (assuming a circular orbit), one needs to specify the time of mid-transit, the orbital period, the inclination (or impact parameter), the radius of the star, the ratio of the stellar and planetary radii, and the limb darkening parameters. Given these parameters, the Mandel & Agol (2002) algorithm can be used to produce a model light curve. Obviously, in the case of Kepler-47, the primary is part of a binary system and not a single star. Nevertheless, the shapes of the transits can still be matched with a suitable change in the orbital period of the simple model. Because the transit shapes are well-matched, our model of a single planet transiting a single-star should produce reliable times of mid-transit.

The model fits to each individual transit were optimized using both a genetic algorithm and a Differential Evolution Monte Carlo Markov Chain (DE-MCMC) routine (Ter Braak, 2006). Once a best-fitting model was found, the uncertainties on the observations were scaled to give a reduced $\chi^{2}$ of 1 and the model was re-optimized. The lower ($\sigma_{\rm low}$) upper ($\sigma_{\rm high}$) $1\sigma$ uncertainties were taken to be the interval in the time where $\chi^{2}=\chi^{2}_{\rm min}+1$ on the low and high sides of the best-fitting time, respectively. Table 2 gives the measured times and durations along with their uncertainties, where the adopted uncertainty on each time was the larger of $\sigma_{{\rm low}}$ or $\sigma_{{\rm high}}$. For completeness, we also give the corresponding times and durations found from the best-fitting photodynamical model described below.

There is a point that was affected by a cosmic ray (data quality flag 8192) close to the middle of the transit of the middle planet, near day 418.0 (see Figure 2). That point was not used in the analysis. When that transit is plotted with all of the points connected by lines and without the bad point, it appears the transit is $\approx 30-45$ minutes late relative to the model. The formal uncertainty on the measured transit time with the bad point excluded is about none minutes. There is a corrected flux value for the time in question, and the black point with the error bar in Figure 2 shows where that corrected flux would have been had that observation been in the detrended light curve. The transit looks less convincing when that point is shown, because its flux value is a bit higher than the flux values on either side. We did some Monte Carlo simulations where the bad point was put back into the light curve, but with a flux value drawn from a Gaussian distribution centered on the corrected value and with a standard deviation equal to the reported uncertainty. The measured transit time was shifted earlier by about three minutes from where it was without the questionable point included.

As discussed in Orosz et al. (2012b), the planet transit times are not expected to follow a simple linear ephemeris–owing mostly to the motion of the primary star around the system barycenter (because the location of the primary star changes as it orbits, transits can occur early or late relative to a stationary star), and to a lesser extent, perturbations from the binary (Agol et al., 2005). We fitted each set of transit times to linear ephemerides, and formed “Observed” minus “Computed” (“$O-C$”) diagrams that are shown in Figure 5. The middle planet has $O-C$ variations of a few hours, which is strong evidence that the transits are due to a circumbinary body and are not from a background blend (a common false-positive scenario for planets around single stars). The inner planet has maximum $O-C$ variations of about 25 hr, and the outer planet has maximum $O-C$ variations of nearly 50 hr. Moreover, the $O-C$ variations for all three planets vary cyclically on the binary star’s orbit. We define phase 0.0 as the time of the primary eclipse. Given this definition, the $O-C$ variation is zero at phases 0.0 and 0.5 because the projected displacement of the star from the center of mass is at a minimum. Near phases 0.25 and 0.75, the star is at the projected ends of its orbit, and the timing variations are largest. If the binary system and the planet both have circular orbits, then the $O-C$ variation with orbital phase will be sinusoidal; see Figure 5. Note that the transits of the middle planet span a limited range in binary phase due to the roughly 25:1 ratio of the planet’s orbital period to the binary orbital period, so only a small portion of the expected sinusoidal variation is seen. This limited range in binary phase has no adverse effect on the mass determination from the photodynamical model discussed below.

Finally, the planet transit durations depend on the relative projected velocities of the primary star and the planet. If a transit occurs near primary eclipse, the bodies are moving in opposite directions, resulting in a narrow (short-duration) transit. For transits near secondary eclipse, the bodies are moving in the same direction, yielding a longer duration transit; see Figure 5. Note that, if the planet transits near the stellar limb rather than the center, then the duration can be short at any phase–this is the case for the middle planet. Kostov et al. (2014) give an analytic expression to compute the transit durations as a function of binary phase, assuming a circular orbit for the planet and a constant impact parameter. Figure 5 shows the model duration curves for each planet, where we used $b=0.45$, 0.63, and 0.92 for the inner, middle, and outer planets, respectively. The model curves do a reasonable job of fitting the observed durations, although the scatter is somewhat large given that individual transits can have an impact parameter that is quite different than the mean; see Table 2.

4 Stellar Eclipse Times and

Corrections

We measured the times of the primary and secondary eclipses using the technique outlined in Welsh et al. (2012) and Orosz et al. (2012b), and the results are given in Table 3. Note that the cycle numbers given for the secondary eclipses are not exactly half integers because the orbit is eccentric.

The times of the primary eclipses were fitted to a linear ephemeris. The $O-C$ residual times were computed and are shown in Figure 6. As was the case in our earlier work (Orosz et al., 2012b), there are coherent modulations of up to about two minutes, with a quasiperiod near 178 days. These timing anomalies are caused by star spots on the primary that are partially or fully covered during the primary eclipse. When this occurs, the eclipse profile is not symmetric in time. As discussed in Orosz et al. (2012b) and in Welsh et al. (2015), the shift in the measured eclipse time will depend on where the spot is seen on the face of the primary star and when the secondary passes over it. Also, the slope of the light curve near primary eclipse depends on where the spot is on the primary as it rotates into and out of view. Therefore, a correlation is expected between the local light curve slope and the $O-C$ residual (Holczer et al., 2015; Mazeh, Holczer & Shporer, 2015), which is shown in Figure 7. A linear function was fitted to the data shown in Figure 7 and used to statistically correct the times of the primary eclipses. A new linear ephemeris was fitted to the corrected eclipse times, and the $O-C$ residuals are shown in the middle panel of Figure 6. Apart from a feature near day 1000, the scatter in the residuals has been greatly reduced compared to the the residuals shown in the top panel of the figure. No spot-induced variations were seen in the times of the secondary eclipse, so no correction to those times was applied.

The best-fitting ephemerides for the corrected primary eclipse times and the secondary eclipse times are:

[TABLE]

The difference between the primary and secondary periods is $0.02\pm 0.17$ s. The fitted period from the secondary eclipses is formally longer than the fitted period from the primary eclipses, although the difference is not significant. This lack of a significant difference between the two periods means the effect of the planets on the eclipse timings is unmeasurable, and the masses of the planets are constrained by planet-planet interactions, not planet–binary interactions (see §6.1).

5 Photodynamical Modeling of the Light and Velocity

Curves

The complete Kepler-47 light and velocity curves were modeled using the ELC code. Given five bodies with some initial positions and velocities, the Newtonian equations of motion can be integrated, thereby giving the positions and velocities of each body as a function of time. Then, given the position of each body on the plane of the sky at a specific time (corrected for light travel time) and information about their radiative properties, the light curve can be computed. In the discussion below, we give details of the photodynamical model and its application to the Kepler-47 data.

5.1 ODE Integrator

The numerical integrator that ELC uses is a symplectic, 12th-order Gaussian Runge-Kutta (GRK) routine based on methods and codes devised by Hairer, Lubich & Wanner (2006). For situations where only the mutual gravitational forces between the five bodies are considered, the ordinary differential equations (ODEs) are second-order and can be written such that only the positions of the bodies appear explicitly. GRK is a collocation method in which a system of equations based on the Gaussian quadrature nodes must be solved at each time step. In particular, for the 12th-order method, there are six nodes that determine the coefficients of the collocation polynomial. These coefficients are solved for iteratively. With a suitable initial guess, the convergence of this iteration is quadratic–and hence, very fast. If needed, the equations of motion can be modified to account for extra effects such as precession due to General Relativity (GR) and/or precession due to tidal bulges on the stars (Eggleton, Kiseleva & Hut, 1998; Mardling & Lin, 2002). In the case of the GR correction, the velocities of the bodies appear explicitly in the equations of motion, and the first-order ODEs must be solved using a different iteration scheme. This iteration scheme is only to first order, so the convergence is somewhat slower than the second-order scheme.

Many of the symplectic integrators used in the field of solar system dynamics (e.g., Wisdom & Holman, 1991) rely on a factorization (or splitting) of the Hamiltonian in a Lie algebra sense. These integrators are very fast, but are best applied to problems where there is one dominant central mass. One advantage of GRK integrators is that they do not approach the issue as a Lie algebra factorization and hence do not have the restriction of just one large mass. The trade-off is that the step size is determined by the shortest period, which results in longer run-times (we typically use a time step $h$ that is about 400 times smaller than the smallest period in the system). However, we note that, even with the relatively low speed compared to other integrators, it takes much less than one second to solve the coupled equations of motion for the five bodies in Kepler-47 over a 1500 day time span using the GRK method.

The GRK scheme does have another advantage, in that it is possible to find positions and velocities at intermediate times to high accuracy with relatively few computations. In our implementation, the solutions at the six internal nodal points used for the collocation for each time step are saved. Using the collocation points with divided differences, it is easy to determine the positions and velocities of each body at intermediate times to the same order of accuracy (namely, 12th) as at the integration times (e.g., $t_{\rm start}$, $t_{\rm start}+h$, $t_{\rm start}+2h$, etc.). The number of operations needed to achieve the 12th-order accuracy is comparable to what one would need for a spline interpolation scheme, which would give the intermediate values with a much lower-order accuracy.

The integrator works in Cartesian coordinates relative to the system’s barycenter. The unit of mass is the solar mass, the unit of distance is the astronomical unit, the unit of time is the day, and the unit of velocity is astronomical units per day. The adopted value of the Gaussian gravitational constant is $k=0.01720209895$ (Clemence, 1965). Newton’s gravitational constant, expressed in $({\rm AU})^{2}{M_{\odot}}^{-1}({\rm day})^{-2}$, is then $G=k^{2}$. For convenience, instantaneous Keplerian parameters valid at some reference epoch for each body are used to specify the initial conditions (our adopted reference time is $T_{\rm ref}=$ BJD 2,454,965.000). These parameters are the orbital period $P$, the time of barycentric transit (which is the same as the time of inferior conjunction) $T_{\rm conj}$, the inclination $i$, the eccentricity $e$, the argument of periastron $\omega$, and the nodal angle $\Omega$. The coordinate system is Jacobian, so that the orbital parameters of the secondary star are given relative to the primary, the orbital parameters of the first planet are given relative to the binary’s center of mass, and so on. Because of this, the specified conjunction time for the binary occurs very near a primary eclipse, whereas nothing observable has to happen at the conjunction times of the planets (the actual transits can occur either earlier or later than the barycentric conjunction times). The Keplerian parameters are converted to Cartesian coordinates using the algorithms given in Murray & Dermott (1999).

After the ODE is solved, the plane-of-sky positions of each body are corrected for light travel time following the method outlined in Carter et al. (2011). The times of the primary eclipses and secondary eclipses, the times of transits of each planet across the primary secondary are then computed by finding the times when the projected separation between two given bodies in the plane of the sky is minimized. For computational convenience, we minimize the square of the plane-of-sky distance vector, which is just the dot product of that distance vector with itself. That dot product is minimized when its derivative is zero. That derivative is two times the dot product of the distance vector with its velocity vector (which itself is readily found using the appropriate components of the velocities that are output from the ODE solver), and is zero when the distance vector itself is zero or when the distance vector is perpendicular to its velocity vector. Several iterations of the secant method are used to find the zeros of the derivative. By comparing the results found using the first-order ODE solver with the results found using the second-order ODE solver, we estimate the times of conjunction to be accurate to $\approx 1\mu$s or better. An eclipse or a transit will occur at a conjunction if the projected minimum separation on the sky is less than the sum of the radii of the two bodies.

For the specific case of Kepler-47, the corrections due to GR result in an apsidal advance of about 0.00017 degrees per orbital period. The rate of apsidal advance due to tidal bulges on each star is between 10 and 50 times smaller, depending on the values of the tidal Love numbers used (we considered $0.0\leq k_{2}\leq 0.01$ for the primary and $0.0\leq k_{2}\leq 0.2$ for the secondary). In what follows, we have included the GR correction in the ODEs and neglected the corrections for the apsidal advance due to tidal bulges on the stars.

5.2 Light Curve Synthesis

For situations where the bodies are spherical with linear or quadratic limb darkening laws, ELC can use the algorithm of Mandel & Agol (2002) or the algorithm of Giménez (2006a) to compute the light curves during eclipses or transits. In the present work, we have used the Mandel & Agol (2002) routine with a quadratic limb darkening law, because it is much faster than the Giménez (2006a) routine.

5.3 Model Setup for Kepler-47

Our model as applied to Kepler-47 has 42 free parameters in total. Each orbit needs six Keplerian parameters for a total of 24. However, the nodal angle of the binary is fixed at zero, so there are really only 23 free parameters to describe the orbits. For each orbit, we used the combinations $e\cos\omega$ and $e\sin\omega$ rather than $e$ and $\omega$ separately, because the former pair of variables usually are less correlated with each other compared to the latter pair of variables. Because there are five bodies, five parameters describing the masses of the bodies and five parameters describing the sizes of the bodies are needed. For the binary, we used the primary mass $M_{1}$ and the binary mass ratio $Q_{2,1}\equiv M_{2}/M_{1}$. For the two stars, we used the fractional radii $R_{1}/a$ and $R_{2}/a$ to parameterize their sizes. To parameterize the planetary radii, we used the ratio $R_{1}/R_{p}$, which is the ratio of the primary star’s radius to the particular planet radius in question. There are two radiating bodies in the system, so a total of six parameters are required to specify the radiative properties: the two temperatures and two limb darkening coefficients for each star. In this case, the temperature of the primary was fixed at 5636 K (Orosz et al., 2012b), so only five free parameters were used. For convenience, the temperature ratio $T_{2}/T_{1}$ was used instead of directly using $T_{2}$. ELC has a table of model atmosphere specific intensities (for solar metallicity) derived from the NextGen models (Allard et al., 1997; Hauschildt, Allard, & Baron, 1999), and these were used to compute the baseline fluxes of the primary and secondary star (in the Kepler bandpass) given their temperatures and gravities. The standard quadratic limb darkening law given by

[TABLE]

was used (where $\mu=\cos\theta$ is the projected distance from the center of the stellar disk), but with the “triangular” sampling technique of Kipping (2013) with coefficients given by $q_{1}=(u_{1}+u_{2})^{2}$ and $q_{2}=0.5u_{1}(u_{1}+u_{2})^{-1}$. Finally, to account for light from other sources in Kepler’s aperture, four seasonal contamination parameters were used.

5.4 Model Optimization

As was shown in §4, many of the primary eclipse profiles show clear evidence of spot crossing events. Because our light curve model does not include star spots, fitting the primary eclipses where there are clear spot crossing events might lead to biased results. To avoid this possible bias in the fitting, we fit the statistically corrected primary eclipse times along with only a small subset of the light curve. We chose ten primary eclipse profiles (seven in long cadence and three in short cadence) that do not have spot crossing events to fit, as well as ten secondary eclipses that are close in time to the ten selected primary eclipses. We also fit the 25 transits of the inner planet, the six transits of the middle planet, and the four transits of the outer planet. Finally, we also included the radial velocities of the primary star in the fit.

We used four algorithms to explore parameter space: an adaptation of a simple grid-search algorithm (Bevington, 1969), a genetic algorithm (Charbonneau, 1995); a simple search based on Monte Carlo Markov Chain (MCMC) algorithm (Tegmark et al., 2004), and a Differential Evolution MCMC (DE-MCMC) algorithm (Ter Braak, 2006; Nelson et al., 2014). All four of these algorithms are bounded, meaning the fitting parameters never wander outside predefined ranges. In addition, the grid-search code, the genetic code, and the DE-MCMC code can be run in parallel on multiple CPU cores. We assumed the likelihood function follows a $\chi^{2}$ distribution, where $\chi^{2}$ is based on the sum of squares of differences between observed quantities (a photometric measurement, a radial velocity measurement, or an eclipse or transit time) and predicted quantities normalized by the measurement uncertainty for each quantity. Each of these algorithms needs to have fixed ranges for each free parameter. Some parameters, like some of the orbital parameters for the binary, are constrained reasonably well. We also had some good starting solutions for the parameters of the inner planet from our previous work (Orosz et al., 2012b), whereas some of the orbital parameters (for example, the eccentricity parameters) for the outer planet were constrained much less well. Apart from the rough value of the orbital period and the fact that the inclination of the orbit must be near $90^{\circ}$ in order for transits to occur, the orbital parameters for the middle planet were initially not known.

Because the parameter space is vast, we looked for optimal solutions in three main stages: (i) updating the model from the previous work by adding the middle planet to the model and finding some initial fits; (ii) exploring parameter space to find something close to the “best” model, get some rough idea of the uncertainties in the fitted and derived parameters, and also look for possible correlations between various parameters; and (iii) finding realistic uncertainties on the fitted and derived parameters. The three stages are described below.

Stage i: The inner planet has seven additional transits and the outer planet has one additional transit since the previous work, so the orbital parameters for those two planets were revised first using the DE-MCMC code. Next, the middle planet was added to the model, and initial orbital parameters were found “by hand.” Educated guesses for the time of barycentric conjunction, the inclination, the nodal angle, etc. were used to initialize the simple MCMC code, and the model fits were inspected by eye after running a few iterations of that optimizer. It was not unduly difficult to find a model that produced transits of the middle planet at roughly the correct times.

Stage ii: After a few decent models were found, longer runs using the genetic algorithm and the DE-MCMC were done. Both of these codes require an initial “population” of models (we typically used population sizes between 100 and 200), and we have various ways of generating the initial population: (a) totally random values between the specified lower and upper parameter bounds for each free parameter; (b) one or more “elite” models (e.g. ones that are already a good match to the data) along with randomly generated models; and (c) one or more “elite” models along with “mutated” copies of the elite models where a few randomly chosen parameters are “tweaked” by adding or subtracting small offsets. In practice, options (b) and (c) work the best for cases where one has a few dozen or more free parameters.

For this intermediate stage of fitting, we adopted ranges for the free parameters given in Table 4. We ran the genetic algorithm or the DE-MCMC code for a few pilot runs of several thousand generations to confirm that the prior ranges included support for the entire range with nontrivial likelihood. After each pilot run, plots of parameter values versus the generation number and plots of the $\chi^{2}$ versus each parameter were inspected to find instances where the allowed range of a parameter was either too large or too small. We verified that the likelihood falls to extremely small values by the time model parameters reach any of these boundaries (with the possible exception of hard physical boundaries).

After this intermediate state of optimization, several million models have been computed. The genetic algorithm and the DE-MCMC code explore parameter space in very different ways; and after several long runs of each, we can be reasonably sure the optimal region in parameter space has been found.

Stage iii: We used a brute-force “stepper” to generate initial “seed” models for the final run of the DE-MCMC code. The stepper works as follows. Take the best overall model and choose a key fitting parameter, such as the primary mass. Offset that parameter by a small amount (say, 0.5% of the parameter value) and hold it fixed while optimizing the other parameters using the simple parallel grid-search code, where the parameter values from the optimal model are used as the initial guess. After that optimization is done, offset the same key parameter again and repeat the optimization, using the previous solution as the initial guess. After the stepper is done, you have a set of reasonably optimal models with a range of different values of the chosen key parameter. We ran steppers on the mass parameters (primary mass, binary mass ratio, and the planet masses), the radius parameters (primary and secondary fractional radii and planet radius ratios), and the planet orbital periods.

For the last run of the DE-MCMC code, we used a population size of 1600 models per generation with 900 “seed” models found from the steppers. The other 700 models were mutated copies of randomly selected seed models with five parameters in each model randomly tweaked. The code was ran for 17,000 generations, giving about 27 million models in total. Judging from plots of the parameters vs. the generation number, convergence was achieved after about 4000 generations. However, to be conservative, the first 5000 generations were discarded, leaving about 19 million models. Figure 1 shows the best-fitting model radial velocity curve. The best fits to the transits of the middle planet, inner planet, and the outer planet are shown in Figures 2, 3, and 4, respectively. Figure 8 shows the best-fitting model with the “clean” primary eclipses and Figure 9 shows the best-fitting model with the corresponding secondary eclipses. The best-fitting input parameters are given in Table 5, and Table 6 summarizes several derived parameters of interest. In order to allow others to reproduce our best-fitting model using other codes, we give the initial Cartesian barycentric coordinates for the best-fitting model in Table 7, and the initial orbital elements (which are traditionally used in dynamical studies) for the best-fitting model in Table 8. These coordinates and orbital elements are valid for the reference time BJD 2,454,965.000.

For each parameter’s posterior distribution, we computed the mode and the lower and upper boundaries of the region that contains 68.3% of the area under the curve. We show plots of posterior distributions for several fitting and derived parameters in Figures 10-17, and plots of the two-parameter joint posterior distributions for the 28 major orbital parameters in Figures 18-20 (Table 9 gives the displayed parameters and their ranges).

When one examines the posterior distribution of a given parameter, the mode of the distribution might be taken as the “most likely” value of that parameter. We note, however, that a model constructed using the modes of each distribution will not necessarily provide a good fit to the data. Therefore, for the final adopted parameters, we take the parameters from the best-fitting model.

6 Discussion

6.1 Constraints on the Planetary Masses

For several CBPs (Doyle et al., 2011; Welsh et al., 2012; Kostov et al., 2016), the gravitational pull of the planet is sufficiently large to measurably alter the positions and velocities of their host stars. These perturbations can be detected by measuring slight deviations from strict periodicity in the times of the stellar eclipses. From these eclipse-timing variations, the mass of the planet can be estimated. For Kepler-47, only upper limits on the masses could be made at the time of their discovery: $<$ 2.0 Jupiter masses ($<635\ M_{\oplus}$) for the inner planet and $<$ 28 Jupiter masses ($<8900\ M_{\oplus}$) for the outer planet (Orosz et al., 2012b).

The presence of the third planet dramatically changes this situation. The two outer planets gravitationally perturb one another, resulting in dynamical transit-timing variations. Due to their proximity, the planets perturb each other more than they affect the stars. In the discovery paper (Orosz et al., 2012b), the planets were assumed to be massless in the modeling process, as the interactions between the planets and the binary–and between each other–was determined to be small. In the present model, all five bodies have mass and the mutual interactions between each body are fully accounted for. We find that we can place meaningful constraints on the masses of the middle and outer planets. The $1\sigma$ range of the mass of the inner planet is between zero and $25.8\,M_{\oplus}$, whereas the mass of the middle planet is different from zero at the $\approx 5\sigma$ limit and the mass of the outer planet is different from zero at the $\approx 2\sigma$ limit.

We computed the planet transit times of the best model with a zero-mass middle planet and the best model with a zero-mass outer planet, and compared those times to the transit times from the overall best-fitting model. The difference in the respective transit times are shown in Figure 21. In either case, the transit times of the inner planet do not change by more than $\approx 1$ minute. On the other hand, a massless middle planet results in changes in the outer planet’s transit times by up to $\approx 20$ minutes. Likewise, a massless outer planet results in changes in the middle planet’s transit times by up to $\approx 10$ minutes. These changes in the transit times can be detected at the significance of a few $\sigma$, leading to the mass constraints that we have on the middle and outer planets.

Uncertainties in the planet masses have improved by over an order of magnitude compared with the discovery paper (Orosz et al., 2012b). The mass estimates for the inner and outer planets are now $<26\,M_{\oplus}$ and $\sim 2$-$5\,M_{\oplus}$, at the $1\sigma$ level. The middle planet has a mass $\sim 7$-$43\ M_{\oplus}$ ($1\sigma$ range). As usual, the radii are determined much better than the masses: 3.05, 7.0, and 4.7 $R_{\oplus}$ for the inner, middle, and outer planets, respectively. The middle and outer planets have low bulk densities, $<0.68$ and $<0.26$ g cm*-3* at the $1\sigma$ level (Saturn’s density is 0.69 g cm*-3*). The densities of all currently known low-mass planets are shown in Figure 22; models indicate that such low-density planets must have substantial hydrogen and helium atmospheres (Lopez & Fortney, 2014; Jontof-Hutter et al., 2016). There is also an apparent trend that highly irradiated planets tend to have high densities while planets with low incident fluxes have a range of densities. Kepler-47c and Kepler-47 d fit this latter trend.

Several studies have shown that the formation of CBPs at close distances to the binary is not expected to proceed efficiently, and that CBPs are likely formed at large distances and migrated to their current orbits (Pierens & Nelson, 2008, 2013; Meschiari, 2012; Paardekooper et al., 2012; Pelupessy & Portegies Zwart, 2013; Martin, Armitage & Alexander, 2013; Rafikov, 2013; Kley & Haghighipour, 2014, 2015; Lines et al., 2015a). Planets that survive the migration phase are expected to have orbits with small eccentricities and small mutual inclinations relative to the orbital plane of the binary stars. The planets in Kepler-47 have low masses compared with Jupiter or Saturn, and the orbits have low eccentricities: $e<$ 0.030, 0.024, 0.05, 0.07 (at the +1$\sigma$ level) for the binary and the inner, middle, and outer planets, respectively. All four orbits have mutual inclinations aligned to within 1.6 degrees of one another (+1$\sigma$ level). This nearly circular, co-planar, packed configuration is unlikely to have arisen as an outcome of strong gravitational scattering of the planets into their current orbits. Rather, the observations suggest that this planetary configuration is the result of relatively gentle migration in a circumbinary protoplanetary disk.

6.2 Stability of the Orbits

Numerical integrations have shown that the region around the binary where the orbit of a planet will become unstable extends to $\approx 0.18$ au (Hinse et al., 2015). Thus the semimajor axis of the innermost planet (0.2877 au) is much farther beyond the outer edge of the unstable region. The integrations of the observed system (using the nominal masses given in Table 6) also show evidence of stability, as the variations in the semimajor axis, eccentricity, ascending node, and inclination relative to the binary plane, $i_{\rm rel}$, do not have appreciable secular changes and are remarkably flat in their maximum and minimum values over the 100 Myr integration time scale (see Figure 23). The angular momentum deficit (AMD) of the system is used to estimate the changes of secular perturbations that could induce instability on a long time scale. We scale this quantity relative to the present-day solar system value for the terrestrial planets. The Kepler-47 planets are more massive than our terrestrial system, thereby allowing for higher values in AMD. Finally, we show a periodogram for the cyclical variations in inclination within the system and find the precession time scales of the inner, middle, and outer planets to be $\sim 10$ yr, 245 yr, and 738 yr, respectively (see Figure 23). Figure 24 shows the evolution of the impact parameters (for transits across the primary star) for the three planets. The horizontal dashed lines in the figure denote impact parameters where transits of the primary star could occur. As one might expect, the precessionary motion affects the duration of time for which a given planet can transit relative to our line of sight (Schneider, 1994).

We examined the possibility of additional planets within the Kepler-47 system. Before exploring the possibility of the existence of an actual planet, we began by integrating the orbits of a large battery of test particles in between the system’s planets. The numerical scheme for these simulations used a modified version of the orbital integration package, mercury6, that is specifically designed to evaluate the orbits of CBPs efficiently (Chambers et al., 2002). We chose to investigate the test particles initially on coplanar (relative to the binary plane), circular orbits and varied the initial phase of the test particles between $0^{\circ}$ and $360^{\circ}$, in increments of $2^{\circ}$. Additionally we varied the starting semimajor axis of the test bodies from 0.25 to 1.05 au, such that we could evaluate whether non-transiting planets between either the inner-middle pair or the middle-outer pair could be stable for 10 Myr. Overall, we integrated test particles corresponding to a grid of $801\times 181$ initial conditions. From this study, we found that virtually all of the test particles between the middle and outer planets became unstable (experiencing ejection from the system or a collision with another body), indicating that no stable (prograde) planetary orbits exist between these two planets. In other words, the region between the middle and outer planets is dynamically full (i.e., these planets are in a packed orbital configuration). However, there remained a broad region of parameter space where the test particles were stable between the inner and middle planets. Despite this stability, there is insufficient evidence to suggest that a massive body exists, because of the many dynamical interactions that would influence the observations of the transiting planets and thereby alert us to its presence. There are many regions where the test particles’ eccentricities were significantly increased due to mean motion resonances with the middle planet and we expect these locations to be eroded over time in a similar manner to how the Kirkwood gaps in the solar system are.

Based on analytical studies, we used the estimate of Chambers, Wetherill & Boss (1996), assuming that the binary can be approximated as a single star with a large oblateness, and found that the separation in mutual Hill radii of the inner-middle pair and middle-outer pair to be $32.90\pm 17.44$ $R_{H,m}$ and $12.34\pm 4.43$ $R_{H,m}$, respectively. These values are well above the critical value of $2\sqrt{3}\approx 3.46$ for single-star systems (Gladman, 1993) or the critical range of values (5–7) suggested by Kratter & Shannon (2014) for circumbinary systems. The mutual Hill radii of the middle-outer pair are quite close, even when considering the uncertainties, and this proximity does not allow for any additional planets to orbit between 0.7 and 0.96 au. Based on our numerical tests, assuming nominal mass estimates, the middle and outer planets are “packed.” This conclusion is supported by other numerical studies (Smith & Lissauer, 2009; Kratter & Shannon, 2014; Quarles & Lissauer, 2018).

We have calculated a dynamical Mean Exponential Growth factor of Nearby Orbits (MEGNO) map (Cincotta, Giordano & Simó, 2003; Goździewski et al., 2001) of the middle planet using the best-fit osculating (Jacobian) initial conditions presented in this work. Initial values for the mean anomalies of the binary and each planet were set to zero. The MEGNO factor $\langle Y\rangle$ is a numerical technique to detect and differentiate between aperiodic (chaotic) and quasiperiodic (regular/stable) motion. The computation of the variational vector necessary to calculate MEGNO was done by solving the variational equations of motion in parallel to the five-body equations of motion. These equations of motion were solved using the accurate adaptive time-step ODEX integrator222http://www.unige.ch/$\sim$hairer/prog/nonstiff/odex.f. A given integration was terminated once $\langle Y\rangle$ was greater than five. The result of computing MEGNO over a grid in semimajor axis versus eccentricity space is shown in Figure 25. The color coding follows the numerical value of MEGNO attained at the end of the integration. Yellow colors indicate chaotic and a blue ($|\langle Y\rangle-2.0|<0.001$) color indicate quasiperiodic dynamics of the third planet. The binary and remaining two planets were started at their best-fit osculating Jacobian elements.

Recently, the MEGNO technique was used to study possible positions of the middle planet in the framework of the five-body problem (Hinse et al., 2015). Compared to that study, the quantitative picture of the semimajor axis vs. eccentricity phase-space structure has changed significantly in the present analysis. The main difference is that the eccentricity of the outer planet has decreased from 0.41 to almost zero in the newly found best-fit model. This demonstrates the overall destabilizing effect an eccentricity of 0.41 for the outer planet orbit has, as it would greatly limit the number of possible stable orbits of the middle planet. This is a fine example of the possibility of planet packing as demonstrated in Kratter & Shannon (2014). Finally, the location of mean-motion resonances (vertical structures in Figure 25) might change for different initial phases of the planets (when varied within their parameter uncertainties). However, the overall global topology structure of chaotic/quasiperiodic orbits of the third planet would not change dramatically.

While there is no dynamical evidence for an additional planet in the Kepler-47 system, another (low-mass) planet between the inner and middle planets could exist. Also, more planets could have stable orbits sufficiently far outside the orbit of the outer planet. Thus, it would be worthwhile to obtain additional high-precision photometry of Kepler-47 whenever possible, to look for additional transits. From Figure 24, we see that the precession timescales for the three known planets are relatively short, and that the middle and outer planets can transit the primary star for only a small fraction of their precession cycles ($\sim 23\%$ and $\sim 12\%$ of the time for the middle and outer planets, respectively). If there are additional planets in the Kepler-47 system with similar precession cycles, they might precess into view at a later time.

6.3 Comparison with Stellar Evolution Models

Our new values for the masses of the stars in Kepler-47 are about 6% smaller than in the discovery paper, and the new radii are about 3% smaller than previously reported. We therefore update our analysis and discussion of the evolutionary state of the Kepler-47 binary. Figure 26 compares the new measurements against Dartmouth models (Dotter et al., 2008) for the nominal spectroscopic metallicity of $[{\rm Fe/H}]=-0.25$ (Orosz et al., 2012b), in the same way as we did previously. There is still an age discrepancy for the primary star: the mass/radius combination yields an age of 3.5 Gyr (solid line in top panel), whereas the mass/temperature combination gives 11.5 Gyr. If the radii are given more weight than the temperatures, then the secondary star is now slightly larger than the theoretical predictions, as is seen in many other M dwarfs. As we noted in the discovery paper, the primary temperature derived from color indices points to a higher value than the spectroscopy does (which would actually be consistent with the 3.5 Gyr age). However, given the uncertainties in the reddening, the spectroscopy is perhaps more trustworthy.

By way of comparison, Figure 27 shows the Dartmouth models for solar metallicity. The primary age from mass/radius is now 4.5 Gyr (solid line), and mass/temperature gives an age consistent with this within the error bar. The secondary radius would be perfectly consistent with the prediction for this metallicity, as would the secondary temperature. However, solar metallicity is formally ruled out by the spectroscopic analysis (Orosz et al., 2012b) at the 3$\sigma$ level, so a scenario where the stars have solar metallicity is unlikely.

We note also that, in its latest data release (DR2), the Gaia mission (Gaia Collaboration et al., 2016, 2018) has provided an accurate parallax for Kepler-47 (source ID 2080506523540902912) of $\pi=0.948\pm 0.029$ mas, corresponding to a distance of $1025^{+31.5}_{-29.7}$ pc (Bailer-Jones et al., 2018), that enables a check on our dynamical radius measurement for the primary star. With a reddening estimate to the system of $E(B-V)=0.056\pm 0.020$ from Green et al. (2015), the assumption that $A_{V}=3.1E(B-V)$, a bolometric correction from Flower (1996) based on the temperature, and the visual magnitude of the star ($V=15.395\pm 0.010$, Henden et al., 2015), we obtain $R=0.928\pm 0.072\,R_{\odot}$. This is less precise than our dynamical value of $R=0.936\pm 0.005\,R_{\odot}$, but is in excellent agreement.

6.4 The Habitable Zone

The habitable zone (HZ) is usually defined as the region around a star such that the stellar incident flux on the planet (i.e. the insolation) would allow an Earth-like planet with a CO2/H2O/N2 atmosphere to permanently maintain liquid water on its solid surface. The situation is more complicated for a CBP system, in that (i) there are two stars, typically of different temperatures and luminosities; and (ii) the orbital motion of the stars causes the HZ to “orbit” with the binary. Even for a planet on a purely circular orbit, the planet–star distances are always rapidly changing. The HZ is no longer static or spherically symmetric.

Therefore, for a CBP system, the definition of the HZ has to be slightly generalized. Here, we define a circumbinary HZ as the region in space where the time-averaged weighted insolation would allow conditions conducive to liquid water on the surface of an Earth-like planet. The weighted insolation is not simply the total flux incident at the top of the planet’s atmosphere (though it reduces to this if one of the stars’ flux contribution is negligible compared to the other, as is the case for Kepler-47). Adding the flux from each star is not appropriate because it is not simply a matter of the total energy; the spectral energy distribution is highly relevant. The planet’s atmosphere acts to filter the energy received from each star before it reaches the planet’s surface, and this filtering is strongly wavelength-dependent. The spectral energy distribution of the light emitted by the star can be reduced to a simple function of the star’s effective temperature ($T_{\rm{eff}}$) to a sufficient level of approximation. Thus, the $T_{\rm{eff}}$-weighted sum of the fluxes from each star is used as the total flux received by the planet (Haghighipour & Kaltenegger, 2013). To determine the inner and outer boundaries of the circumbinary HZ, we equate the weighted incident flux with the insolation corresponding to those boundaries for a single star as computed by Kopparapu et al. (2013, 2014).

Figure 28 shows a face-on view of the Kepler-47 system. The binary stars are shown in black, the planets in blue, the HZ in green, and the red circle marks the critical instability radius (Holman & Wiegert, 1999). The darker green shaded area corresponds to the conservative HZ as defined by the runaway greenhouse and maximum greenhouse conditions. The lighter green regions show the optimistic HZ boundaries, as defined by the recent Venus and early Mars conditions (Kopparapu et al., 2013). The arrow shows the direction of the line of sight from the observer to Kepler-47. The sizes of the stars and planets are not to scale, and their locations are for the epoch of the orbital elements presented in Table 5. The outer planet is within the limits of the conservative HZ, and Planet D (the middle planet) skirts the inner boundary of the optimistic HZ333An animated version of this figure can be found at http://astro.twam.info/hz-ptype/.

In Figure 29, we show the total (unweighted) insolation incident at the top of the atmosphere of the middle planet (in blue) and the outer planet (in black) as a function of time. The dotted lines show the limits of the optimistic HZ and the solid green lines mark the limits of the HZ defined by the moist greenhouse condition. The dashed line shows the boundary corresponding to the runaway greenhouse condition, which is also often used as the boundary for the inner edge of the conservative HZ. The three panels show three different time scales. The left panel spans two binary star periods, and the variations in insolation are dominated by the orbital motion of the primary star. Note the primary and secondary eclipses, as seen from the location of the planets. The middle panel shows one orbital period of the outer planet, and the insolation variations now show the effect of the eccentricities of the planet’s orbit. The right-hand panel shows five orbits of the outer planet. The time-averaged total insolation for the outer planet is $S=0.865$ in Sun–Earth units, with an rms variation of 7.2% about the mean. For the middle planet (d), the time-averaged insolation is $S=1.67$, with an rms variation of 5.5%. For comparative purposes (only), if one assumes a Bond albedo of 0.34 (similar to that of Jupiter or Saturn; see de Pater & Lissauer, 2015), and complete redistribution of the incident energy, the equilibrium temperature is $T_{\rm eq}\approx 241$ K for the outer planet. In the course of this discussion, we must keep in mind that these are low-density planets, with thick H, He atmospheres, and therefore not likely to be habitable.

6.5 Times, Durations, and Impact Parameters of Future Eclipses

and Planet Transits

Ground-based observations of eclipse and transit events would be useful in constraining the dynamical model of the Kepler-47 system. The times of primary and secondary eclipses are well-fit by the linear ephemerides given in Equation (1). Because the times of the planet transits do not follow a linear ephemeris, we computed (using ELC) the times, impact parameters, and depths of future transits (for the years 2019 through 2026) for each planet using the best-fitting model given in Table 5. The uncertainties on these quantities were estimated using sets of initial conditions drawn from the posterior samples of solutions found from the DE-MCMC run, and are given in Table 10. The formal uncertainties on the predicted times generally vary between about one to two hours for the inner planet, and about 0.75 to two hours for the middle planet and outer planet. With a depth of about 15%, the primary eclipses should be observable by modest-sized ground-based telescopes. On the other hand, owing to the relative faintness of the system ($V=15.395\pm 0.010$, Henden et al., 2015), the 0.5% deep secondary eclipses will be more difficult to observe. With depths between about 0.1% and 0.4%, the transits of three planets will be more challenging still. According to our best-fitting model, the inner planet always transits during its $\sim 10$ yr precession cycle, the middle planet transits about 23% of the time during its $\sim 245$ yr precession cycle, and the outer planet transits about 12% of the time during its $\sim 738$ yr precession cycle. We used the REBOUND software (Rein & Liu, 2012) to integrate a set of 368 initial conditions drawn from the posterior sample to estimate the uncertainties on the predicted transit times near the times when the transits stop. We find that the uncertainties of the predicted transit times for the middle planet are about 1.4 days by the time this planet stops transiting in the year $2037.7\pm 3.4$ (because the transits are grazing before they stop, the predicted stopping date has a much larger spread). Similarly, the uncertainties of the predicted transit times for the outer planet are about 2.1 days by the time this planet stops transiting in the year $2049.1\pm 6.8$.

7 Summary

We have identified a third planet in the Kepler-47 system. Although Kepler-47 d is the largest of the three planets, it was not previously detected because its transits were much weaker in the early Kepler data, as the transit impact parameter was near unity. The precession of its orbital plane is such that the transits are currently growing deeper with time. We modeled the Kepler light curve and radial velocity curve simultaneously using a five-body photodynamical model with 42 free parameters. The revised system parameters are generally in agreement with the previously published values (Orosz et al., 2012b). The tight constraints placed by the planets on the orbits of the stars enable highly precise measurements of their masses and radii (see Table 6). The stars have masses of $0.957\pm 0.014$ and $0.342\pm 0.003$ $M_{\odot}$ and radii of $0.936\pm 0.005$ and $0.338\pm 0.002$ $R_{\odot}$, consistent with an age of 3.5-11.5 Gyrs and slightly favoring an age in the younger portion of its range. Uncertainties in the planet masses have improved by over an order of magnitude compared with the 2012 discovery paper (Orosz et al., 2012b). The $1\sigma$ mass estimates for the planets are now $<26\ M_{\oplus}$, $\sim 7$-$43\ M_{\oplus}$, and $\sim 2$-$5\ M_{\oplus}$ for the inner, middle, and outer planets respectively. The radii are much better determined than the masses, and we find radii of 3.05, 7.0, and 4.7 $R_{\oplus}$ for the inner, middle, and outer planets, respectively. The middle and outer planets have low bulk densities, $<$ 0.68 and $<$ 0.26 g cm*-3* at the $1\sigma$ level. With such low densities, these two planets must have substantial hydrogen and helium atmospheres (Lopez & Fortney, 2014).

All of the orbits in the Kepler-47 system have eccentricities less than 0.07 at the $+1\sigma$ level, and all four orbits have mutual inclinations aligned to within 1.6 degrees of one another at the $+1\sigma$ level. This nearly circular, co-planar, packed configuration is unlikely to have arisen as an outcome of strong gravitational scattering of the planets into their current orbits. Rather, the observations suggest that the planetary configuration is the result of relatively gentle migration in a circumbinary protoplanetary disk.

We found that the configuration is dynamically stable for at least 100 Myr. The new planet conforms with expectations based on an earlier analysis of a two-planet system (Hinse et al., 2015). As clearly demonstrated by the appearance and growth of the middle planet’s transits, the planet’s orbit precesses. The precession time scales for the inner, middle, and outer planets are calculated to be $\sim$ 10, 245, and 738 yr. The middle and outer planets are separated by $\sim$ 12 mutual Hill radii, and we verify the earlier estimation (Kratter & Shannon, 2014) that they are dynamically packed. This is the first detection of a dynamically packed region in a circumbinary system, and it further confirms suspicions that planet formation and subsequent migration can proceed much like that around a single star, at least when far from the binary (Pierens & Nelson, 2008, 2013; Kley & Haghighipour, 2014, 2015). We also find that, although they are close to having integer commensurate periods, the middle and outer planets are not in a mean-motion resonance–and yet they are gravitationally interacting and exchanging angular momentum, as indicated by their antiphased oscillations in inclination and eccentricity.

The average insolation (incident radiative flux) the outer planet receives from its two suns is 86.5% of the Sun–Earth value, placing it within the boundaries of the circumbinary habitable zone (Haghighipour & Kaltenegger, 2013). The middle planet’s orbit straddles the “recent Venus” hot edge of the system’s habitable zone. However, given their densities, neither planet is likely to host life.

Kepler was selected as the 10th mission of the Discovery Program. Funding for this mission is provided by the NASA Science Mission Directorate. The Kepler data were obtained from the Mikulski Archive for Space Telescopes (MAST). The Space Telescope Science Institute (STScI) is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NXX09AF08G and by other grants and contracts. This work is based in part on observations obtained with the Hobby–Eberly Telescope, which is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximilians-Universität München, and Georg-August-Universität Göttingen. This research made use of data from the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund and NSF AST-1412587. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, particularly the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. Some of the simulations in this paper made use of the REBOUND code, which can be downloaded freely at http://github.com/hannorein/rebound. J.A.O. and W.F.W. acknowledge support from the Kepler Participating Scientist Program via NASA grant NNX14AB91G and support from NSF grant AST-1617004. J.A.O., W.F.W., G.W., and B.Q. also gratefully acknowledge support from NASA via grant NNX13AI76G. D.C.F. acknowledges support from the Alfred P. Sloan Research Fellowship. D.C.F. and S.M.M. were supported by NASA through the Kepler Participating Scientist Program award NNX14AB87G. N.H. acknowledges support from the NASA ADAP program under grant NNX13AF20G, NASA PAST program grant NNX14AJ38G, and NASA XRP program under grant 80NSSC18K0519. T.C.H. acknowledges KASI research grant #2015-1-850-04. B.Q. gratefully acknowledges support by an appointment to the NASA Postdoctoral Program at the Ames Research Center, administered by Oak Ridge Associated Universities through a contract with NASA. G.T. acknowledges partial support for this work from NSF grant AST-1509375. E.B.F. acknowledges support from NASA Participating Scientists Program awards NNX12AF73G and NNX14AN76G, NASA ADAP program under grant NNX13AF20G, as well as the University of Florida and the Pennsylvania State University’s Center for Exoplanets and Habitable Worlds. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. J.J.L.’s work on this project was supported by NASA’s Astrophysics Data Analysis Program under grant no. 16-ADAP16-0034. Finally, we are deeply grateful to John Hood Jr. for his generous support of exoplanet research at San Diego State University.

\startlongtable

\startlongtable

\startlongtable

\startlongtable

Appendix A Planetary mass and radius data

Table 11 lists the planetary masses, radii, and insolation values used in Figure 22. Note that the Marcy et al. (2014) sample includes many nondetections and upper limits that we do not plot.

We estimate the insolation each planet received from its host star using the cited $T_{\rm eff}$, $R_{*}$ and $M_{*}$ of each host such that:

[TABLE]

where the semimajor axis $a$ is calculated from the measured mass and orbital period (eccentricities are neglected). Uncertainties are added in quadrature.