Two Super-Earths in the 3:2 MMR around KOI-1599

TL;DR

This study confirms the presence of two super-Earths in a 3:2 mean motion resonance around KOI-1599, analyzing transit timing variations to determine their masses, orbital configurations, and resonance stability.

Contribution

It provides the first detailed dynamical analysis of KOI-1599's planetary system, revealing two super-Earths in a stable 3:2 MMR with dual mass solutions.

Findings

Planets are in 3:2 mean motion resonance with anti-aligned apsides.

Mass estimates depend on eccentricity priors, with a preferred smaller-mass solution.

The system's resonance is dynamically robust across different configurations.

Abstract

We validate the planetary origin of the KOI-1599 transit time variations (TTVs) with statistical and dynamical tests. We re-analysed KEPLER Q1-Q17 light-curves of the star, and we independently derived the TTVs. They appear as strongly anti-correlated, suggestive of two mutually interacting planets. We found similar radii of the candidates, $1.9 \pm 0.2\, \mbox{R}_{\oplus}$ for the inner KOI-1599.02, and $1.9 \pm 0.3 \, \mbox{R}_{\oplus}$ for the outer KOI-1599.01. The standard MCMC TTV analysis constrains the planet masses safely below the dynamical instability limit of $\simeq 3\, \mbox{M}_{\mbox{Jup}} $. The best-fitting MCMC model yields $(9.0\pm 0.3)\, \mbox{R}_{\oplus}$, and $(4.6\pm0.3),\ \mbox{R}_{\oplus}$, for the inner and the outer planet, respectively. The planets are trapped in 3:2 mean motion resonance (MMR) with anti-aligned apsides ($\Delta \varpi = 180^{\circ}$) at low-eccentric ($e\simeq 0.01)$ orbits. However, we found that the TTV mass determination depends on eccentricity priors with the dispersion in the (0.01,0.05) range. They permit a second family of TTV models with smaller masses of $\simeq 7\,\mbox{R}_{\oplus}$, and $\simeq 3.6\,\mbox{R}_{\oplus}$, respectively, exhibiting two modes of $\Delta \varpi = 0^{\circ},180^{\circ}$ librations. The 3:2 MMR is dynamically robust and persists for both modes. In order to resolve the mass duality, we re-analysed the TTV data with a quasi-analytic model of resonant TTV signals. This model favours the smaller masses. We also reproduced this model with simulating the migration capture of the system into the 3:2 MMR.

Tables

Table 1. Table 1: Parameters inferred from the analysis of the folded light-curve and their uncertainties. The mid-transit time T mid . subscript 𝑇 mid T_{\rm mid.} and mean period P mean subscript 𝑃 mean P_{\rm mean} from linear ephemeris, planet-to-star radii ratio p 𝑝 p , semi-major axis in stellar units a / R ⋆ 𝑎 subscript 𝑅 ⋆ a/R_{\star} , and we report the limb darkening coefficients q 1 , q 2 subscript 𝑞 1 subscript 𝑞 2 q_{\rm 1},q_{\rm 2} . For all parameters, we estimate the uncertainties as the 16-th, and 84-th percentile of the MCMC samples.

Planet	KOI-1599.02	KOI-1599.01
$P_{\mathrm{mean}}$[d]	13.6164 $\pm$0.0001	20.408 $\pm$ 0.0004
$T_{\mathrm{mid}.}$[BJD-2454900]	60.32 $\pm$ 0.01	52.73 $\pm$ 0.02
$p$	0.0181 $\pm$ 0.0003	0.0180 $\pm$ 0.0004
$a/R_{\star }$	24.89 $\pm$ 0.11	32.63 $\pm$ 0.42
$I$ [deg]	88.60$\pm$0.06	89.78$\pm$0.16
$q_1$	0.2602 (fixed)
$q_2$	0.4096 (fixed)

Table 2. Table 2: Keplerian osculating orbital elements and masses of four representative solutions for KOI-1599 derived by the GEA optimisation algorithm. We also report the best-fitting period and initial mid-transit time for each solution. The reference epoch is Barycentric Kepler Julian Day (BKJD) − - 139 days. The mass of the KOI-1599 star is 1.02 ​ M ⊙ 1.02 subscript M direct-product 1.02\,\mbox{M}_{\odot} (Rowe et al., 2015 ) .

Model/	GEA I (low $e$, $\Delta \omega \simeq \pi$)		GEA II (low $e$, $\Delta \omega \simeq 0$)		GEA III (moderate $e$, $\Delta \omega \simeq \pi$)		GEA IV (high $e$, $\Delta \omega \simeq \pi$)
Planet	KOI-1599.02	KOI-1599.01	KOI-1599.02	KOI-1599.01	KOI-1599.02	KOI-1599.01	KOI-1599.02	KOI-1599.01
$P$[d]	13.6094	20.4386	13.6099	20.4376	13.6099	20.4361	13.6100	20.4355
$T$[d]	74.011	72.961	74.009	72.951	74.012	72.952	74.012	72.953
$m_p\left[{\text{M}}_{\oplus }\right]$	8.3	4.3	7.6	3.7	7.0	3.6	6.8	3.5
$a\left[\text{au}\right]$	0.11230	0.14727	0.11230	0.14726	0.11230	0.14725	0.11230	0.14725
$e$	0.0067	0.0177	0.0062	0.0444	0.0250	0.0500	0.050	0.050
$\omega$[deg]	40.6	-139.9	-177.9	176.7	-24.3	155.6	-31.2	148.8
$ℳ$[deg]	176.8	31.4	34.0	71.4	-116.3	92.3	-107.2	99.4

Table 3. Table 3: The TTV and light-curve model parameters and their uncertainties from the MCMC sampling. For the TTV model, the eccentricity priors σ x i , y i = 0.0033 subscript 𝜎 subscript 𝑥 𝑖 subscript 𝑦 𝑖 0.0033 \sigma_{x_{i},y_{i}}=0.0033 , resulting in Δ ​ ϖ = 180 ∘ Δ italic-ϖ superscript 180 \Delta\varpi=180^{\circ} and χ ν 2 ≃ 1.1 similar-to-or-equals subscript superscript 𝜒 2 𝜈 1.1 \chi^{2}_{\nu}\simeq 1.1 . For all parameters, we estimate the uncertainties as the 16–th, and 84–th percentile of the samples. The T 0 subscript 𝑇 0 T_{0} epoch is BKJD − 139 139 -139  days. The star mass is m ⋆ = 1.02 ​ M ⊙ subscript 𝑚 ⋆ 1.02 subscript M direct-product m_{\star}=1.02\,\mbox{M}_{\odot} , its radius R ⋆ = 0.972 ​ R ⊙ subscript 𝑅 ⋆ 0.972 subscript R direct-product R_{\star}=0.972~{}\mbox{R}_{\odot} (Rowe et al., 2015 ) . We inferred the orbital elements a 𝑎 a , e 𝑒 e , ϖ italic-ϖ \varpi and ℳ ℳ {\cal M} (the mean anomaly at the epoch) from the primary parameters.

Planet	KOI-1599.02	KOI-1599.01
$P$[d]	13.6088 $\pm$0.0006	20.4415 $\pm$ 0.0013
$x\equiv e\cos \varpi$	0.007 $\pm$ 0.003	-0.009 $\pm$ 0.003
$y\equiv e\sin \varpi$	0.009 $\pm$ 0.003	-0.011 $\pm$ 0.003
$T$[d]	74.012 $\pm$ 0.006	72.946 $\pm$ 0.008
mass $m\left[{\text{M}}_{\oplus }\right]$	9.0 $\pm$ 0.3	4.6 $\pm$ 0.3
$a\left[\text{au}\right]$	0.112293	0.147280
$e$	0.0114	0.0140
$\varpi$ [deg]	49.951	230.175
$ℳ$[deg]	167.660	22.141
$R\left[{\text{R}}_{\oplus }\right]$	1.9$\pm$0.2	1.9$\pm$0.3

Table 4. Table SM2: Transit midpoint times, TTVs and relative errors obtained from the MCMC light-curve analysis.

Planet	Transit	Transit time	TTV	$+1\sigma$	$-1\sigma$
[KOI]	number	[BJD-2454900]	[min]	[min]	[min]
1599.02	57	850.07208	-19.886	26.32	27.67
	58	863.68848	-25.747	11.09	11.74
	59	877.30487	-80.165	54.96	52.08
	60	890.92133	-92.448	55.38	27.55
	61	904.53772	16.157	32.64	63.01
	62	918.15411	-35.309	28.63	18.12
	63	931.77051	31.694	30.90	30.87
	65	959.00330	-10.267	19.63	23.08
	66	972.61970	2.851	36.83	33.88
	67	986.23608	-0.346	43.42	38.72
	68	999.85248	41.242	51.36	25.76
	69	1013.46887	40.003	8.12	15.26
	73	1067.93445	-30.370	29.15	27.43
	74	1081.55090	55.022	25.70	17.77
	76	1108.78369	19.613	42.65	67.22
	77	1122.40015	90.288	82.45	61.21
	78	1136.01648	19.886	23.07	26.80
	79	1149.63293	93.802	47.79	34.76
	80	1163.24927	35.021	16.96	15.75
	81	1176.86572	47.333	18.26	16.81
	82	1190.48206	62.467	20.33	91.28
	83	1204.09851	4.6800	19.02	12.90
	84	1217.71484	35.381	25.11	31.38
	85	1231.33130	36.446	32.13	24.77
	86	1244.94775	-9.461	52.92	34.07
	87	1258.56409	92.909	19.60	41.59
	88	1272.18054	57.370	47.48	58.45
	89	1285.79688	100.901	61.20	43.57
	90	1299.41333	26.770	19.12	60.91
	91	1313.02966	61.272	73.15	50.44
	92	1326.64612	47.160	17.01	12.73
	93	1340.26245	29.491	63.12	99.36
	94	1353.87891	-80.093	23.10	31.26
	95	1367.49524	62.712	23.46	36.91
	96	1381.11169	80.611	32.72	37.64
	97	1394.72815	30.442	15.51	21.10
	99	1421.96094	-18.360	31.67	494.97
	100	1435.57727	4.018	23.95	18.48
	101	1449.19373	-5.155	19.63	26.14
	102	1462.81006	2.203	23.88	64.70
	103	1476.42651	12.658	14.01	13.78
	104	1490.04285	-21.082	31.70	35.78
	105	1503.65930	-11.635	12.33	19.02
1599.01	0	73.14220	-263.002	29.48	18.81
	1	93.55050	-242.870	21.77	33.16
	2	113.95880	-205.661	37.12	20.82
	3	134.36710	-181.958	17.97	18.54
	4	154.77541	-165.600	27.56	42.97
	5	175.18370	-77.818	20.39	29.22
	6	195.59200	-91.858	20.13	22.23
	7	216.00030	-40.910	23.54	28.34
	8	236.40860	-48.312	33.67	17.61
	9	256.81690	19.584	80.48	50.67
	10	277.22520	9.677	24.72	21.86
	11	297.63351	9.620	33.22	46.18
	12	318.04181	59.414	14.00	41.34
	13	338.45010	10.757	46.18	17.38
	14	358.85840	69.509	21.72	17.70
	15	379.26670	108.950	29.38	28.18
					Continue

Table 5. Table SM3: Transit midpoint times, TTVs and relative errors obtained from the MCMC light-curve analysis.

Planet	Transit	Transit time	TTV	$+1\sigma$	$-1\sigma$
[KOI]	number	[BJD-2454900]	[min]	[min]	[min]
1599.01	16	399.67500	100.051	15.51	39.77
	17	420.08331	89.410	21.90	24.70
	18	440.49161	122.342	43.59	24.51
	19	460.89990	129.082	16.20	23.24
	20	481.30820	115.474	19.70	27.14
	21	501.71649	94.291	24.71	47.33
	22	522.12482	134.525	22.61	42.60
	23	542.53308	142.891	20.02	16.80
	25	583.34967	153.274	30.27	43.06
	26	603.75800	98.986	15.48	20.53
	27	624.16632	152.554	26.31	8.09
	28	644.57459	126.850	25.76	37.27
	30	685.39117	120.715	22.84	15.25
	31	705.79950	116.525	77.17	13.69
	32	726.20783	105.106	24.68	30.43
	33	746.61609	65.866	37.24	28.93
	34	767.02441	71.755	24.83	51.75
	35	787.43268	75.802	62.48	22.02
	36	807.84100	94.061	72.29	73.74
	37	828.24933	63.245	19.31	11.03
	38	848.65759	47.808	29.88	41.18
	39	869.06592	65.088	36.92	41.64
	40	889.47418	18.677	22.92	31.45
	41	909.88250	-0.562	32.80	23.72
	42	930.29077	5.112	15.11	16.57
	43	950.69910	5.904	18.07	26.68
	44	971.10742	-42.509	31.90	33.24
	45	991.51569	-42.206	28.71	33.93
	46	1011.92401	-51.624	39.14	41.36
	48	1052.74060	-19.598	77.99	21.70
	49	1073.14893	-15.062	21.54	18.81
	50	1093.55725	-53.136	81.56	64.84
	51	1113.96545	-46.843	21.87	14.08
	52	1134.37378	-73.930	21.89	21.18
	53	1154.78210	-82.642	19.21	18.82
	54	1175.19043	-56.938	15.13	22.42
	55	1195.59875	-128.750	27.35	30.90
	56	1216.00696	-68.242	17.19	30.82
	57	1236.41528	-98.654	14.96	30.11
	58	1256.82361	-77.242	11.17	27.06
	59	1277.23193	-107.136	26.60	47.51
	60	1297.64026	-71.280	22.09	16.46
	61	1318.04846	-113.717	33.83	45.85
	62	1338.45679	-28.843	17.27	43.55
	63	1358.86511	-64.094	25.73	13.69
	64	1379.27344	1.454	26.96	8.04
	65	1399.68164	-7.978	23.59	21.33
	67	1440.49829	21.542	10.83	12.99
	68	1460.90662	87.394	35.29	47.29
	69	1481.31494	42.768	64.24	19.24
	70	1501.72314	48.989	21.43	13.38
	71	1522.13147	51.250	38.33	28.32
					Continue

Equations

$$\mathcal{H}=\sum _{i=1,2}\left(\frac{{\mathbf{p}}^2}{2{\beta }_i}-\frac{{\mu }_i{\beta }_i}{|{\mathbf{r}}_i|}\right)-k^2\frac{m_1{m}_2}{|{\mathbf{r}}_1-{\mathbf{r}}_2|}+\frac{{\mathbf{p}}_1\cdot {\mathbf{p}}_2}{m_{\star }},$$

$${\mathbf{f}}_i=-\frac{{\mathbf{v}}_i}{2{\tau }_i}-\frac{{\mathbf{v}}_i-{\mathbf{v}}_{c,i}}{{\tau }_i{\kappa }_i^{-1}},$$

$$t_n=t_n^{\text{mbox}(lin)}+\Delta t_n^{\text{mbox}(s-p)}+\Delta t_n^{\text{mbox}(res)},$$

$$t_n^{\text{mbox}(lin)}=n⟨P{⟩}_{\text{mbox}obs}+T$$

$$\Delta t_n^{\text{mbox}(s-p)}=A_{\text{mbox}s-p}\cos (N_{\text{mbox}s-p}{t}_n^{\text{mbox}(lin)}+{\Phi }_{\text{mbox}s-p}),$$

$$\Delta t_n^{\text{mbox}(res)}=A_{\text{mbox}res}\cos (N_{\text{mbox}res}{t}_n^{\text{mbox}(lin)}+{\Phi }_{\text{mbox}res}).$$

$$\frac{A_{\text{mbox}s-p}}{⟨P{⟩}_{\text{mbox}obs}}=\frac{⟨e⟩}{\pi },$$

$${\Phi }_{\text{mbox}s-p,1}=\frac{3}{2}\pi -{\varphi }_1^{\text{mbox}(eq)}+q{n}_1{T}_1-(q+1)n_2{T}_2,$$

Full text

††pagerange: Two Super-Earths in the 3:2 MMR around KOI-1599–References††pubyear: 2019

1 Introduction

The Kepler mission (Borucki et al., 2011) discovered hundreds of extra-solar multiple planetary systems ((Akeson et al., 2013), https://exoplanetarchive.ipac.caltech.edu/). The distribution of the period-ratio in the Kepler sample (Lissauer et al., 2011; Fabrycky et al., 2014; Delisle & Laskar, 2014) shows a paucity of systems near to first-order mean motion resonances (MMRs), with a significant peak close to the 3:2 MMR. Such features of the period-ratio distribution may be related to the formation history and dynamic evolution of multiple planet systems (e.g., Lithwick & Wu, 2012; Batygin & Morbidelli, 2013; Papaloizou, 2015). It is therefore critical to determine whether a multiple planetary system is dynamically resonant or only close to a strictly resonant configuration (e.g., Petrovich et al., 2013; Goldreich & Schlichting, 2014).

A crucial data source regarding the multiple Kepler planet configurations are the TTV measurements published in recent catalogues by Rowe et al. (2015) and Holczer et al. (2016) (furthermore, H16). The data span the Q1-Q16 quarters of the Kepler light-curves (LCs). Through inspecting these measurements, we selected KOI-1599 with two putative planetary companions, with the inner candidate marked as a possible planet and the outer one not yet examined in the NASA Exoplanet Archive. The TTVs of KOI-1599 exhibit an anti-correlated sinusoidal trend, indicative of two gravitationally interacting objects (e.g., Steffen et al., 2012; Steffen & Hwang, 2015). Our primary motivation for investigating this putative 2-planet configuration is the proximity of their orbital periods to the 3:2 MMR. We did not find any studies aiming to characterise this interesting and likely resonant system.

Since the star is dim ($V\simeq 15$ mag), it would be a difficult target for a spectroscopic follow-up, and we aim to constrain masses of the planetary companions with the TTV orbital model (e.g. Agol et al., 2005; Holman & et al., 2010). Recently, Baranec et al. (2016) imaged $\sim 1000$ dim Kepler stars unsuitable for the spectroscopic follow-up. They detected two nearby dim field stars, yet with a substantial angular separation of $\sim 3$ arcsec from KOI-1599, which may dismiss the blend effect. We aim to verify this furthermore on the dynamical grounds, re-compute the planet-to-star radius-ratio, and the density estimates for the candidate planets.

In this work, we follow Holman & et al. (2010); Nesvorný et al. (2013); MacDonald et al. (2016), as well as Panichi et al. (2018), regarding the dynamical photometry method. In Sect. 2, we re-analysed the whole Q1-Q17 DR-25 Kepler LCs of KOI-1599, and we update the TTVs measurements. In Sect. 3, we validated the two transiting objects as planets. In Sect. 4, we derive the orbital model and masses of the 2-planet configuration. In Sect. 5, we characterise the 3:2 MMR resonant architecture. In Sect. 6 we show that the best-fitting configurations may be interpreted as the natural outcome of planetary migration. We discuss the internal compositions of the planets in Sect. 7. We present our conclusions in Sect. 8. Supplementary Material (SM) with source TTV data is presented on-line.

2 The light-curve analysis and TTVs

Aiming to validate the KOI-1599 planets, we used the Kepler photometric data only. We re-analysed the corrected, de-trended LC-INIT light-curves (LCs) from the DR-25 Kepler release, spanning the whole Q1-Q17 quarters, also in order to verify and, possibly, refine the previous TTVs measurements. To avoid confusion, the inner planet KOI-1599.02 has index “$1$” and the outer planet KOI-1599.01 is labeled with “$2$”, respectively. With the box-least-squares (BLS) algorithm (Kovács et al., 2002), we searched for periodic signals in the LCs. The two transiting objects are apparently close to the 3:2 MMR, thus their mid-transit times may be significantly shifted from the linear ephemeris.

We followed Panichi et al. (2018), to extract the TTVs from the LC-INIT LCs. We split the LCs in fragments and we selected a narrow window ($\pm 0.5$ days) at multiples of the two periods obtained with the BLS search. Each of them should contain only one transit-like signature. Once we re-normalised the out-of-transit parts of each fragments, we superimposed them in order to obtain a folded LC for both planets. We used the exonailer package of Espinoza et al. (2016) for preliminary estimations of the planet-to-star radius-ratio ($p$) and the orbital inclination ($I$). We interpolated the quadratic limb-darkening coefficients based on data in the NASA archive. We kept the photometric noise $\sigma_{w}$ fixed and we estimated it from the off-transit fragments of the LCs. After fixing these parameters, we fitted mid-transit moments ($T$) for each of the LC fragments. We applied this preliminary list of mid-transit times for re-folding each of the fragments, and we iterated the same procedure until no significant differences in the best-fitting parameters are present. The inferred parameters and related uncertainties are listed in Tab. 1. In this way, we controlled the derived TTVs, radii $R_{i}$, $(1.9\pm 0.2)$ and $(1.9\pm 0.3$) Earth radii, as well as orbital inclinations $I_{i}$. The system appears as almost co-planar, since the inner planet has $I_{\rm 1}=88.60\pm 0.06$ [deg], while for the outer planet $I_{\rm 2}=89.78\pm 0.1$ [deg]. We checked that star mass and radius from (Rowe et al., 2015) are in agreement, within the $1\sigma$ uncertainties, when compared with recent estimates in Johnson et al. (2017) and Sanchis-Ojeda et al. (2012).

We derived the median values of the photometric parameters and their uncertainties with the Markov Chain Monte Carlo (MCMC) affine sampler, developed in the emcee package by Foreman-Mackey et al. (2013). We reported the mid-transit times and TTV measurements in the SM on-line.

3 Validation of the planetary TTV origin

We applied tests originally proposed by the Kepler team to validate the majority of candidate planets in multiple systems (e.g. Lissauer et al., 2011; Steffen et al., 2012; Rowe et al., 2014).

The anti-correlation (cross-correlation) of the TTVs can be used to confirm that two objects orbit the same star (Steffen et al., 2012) and their observed anti-correlated TTV signals cannot be a random noise. This method relies on the Monte Carlo bootstrap analysis. We fitted a sinusoidal signal with a fixed modulation period to the TTV series of each planet with the Levenberg-Marquardt algorithm (Press et al., 1986). For each of the sampled periods, the cross-correlation statistics $\Xi$ (Steffen et al., 2012) measures the fit quality — the larger the $\Xi$, the better the quality. Figure 1 shows $\Xi$ for the KOI-1599 TTV signals, as a function of the modulation period. Our new TTVs dataset exhibits a peak at a position similar to the one in the H16 data, but it seems to have even larger significance.

The next step consists of the Monte Carlo bootstrap analysis (e.g., Press et al., 1986). We tested $5\times 10^{4}$ synthetic TTV datasets obtained by random shuffling of the original TTVs, with their uncertainties, and the mid-transit times. As for the real TTV data, we fitted each of the randomly generated datasets with the sine function, and the maximum value of its $\Xi$ is recorded. We illustrate the results in Fig. 1 (middle panel). The False Alarm Probability (FAP) is defined as the ratio of synthetic systems with $\Xi$ larger than that of the observed system to the total number of samples (here, $5\times 10^{4}$). The histogram indicates that unlikely the sinusoidal TTVs are artefacts. The FAP is smaller than $10^{-3}$ adopted by Steffen et al. (2012) for validating other Kepler Objects of Interest (KOIs). We conclude that the KOI-1599 data passes the TTV anti-correlation test.

As a second validation test, we used the dynamical stability constraints. We note that this experiment was not based on the TTVs measurements, rather than on canonical information regarding the orbital periods of transiting objects inferred from the LCs.

In order to conduct the test, we changed the masses, $m_{\rm 1}$ and $m_{\rm 2}$, of the two putative transiting objects in a two–dimensional grid of $512\times 512$ points. At each point, we constructed a number of synthetic configurations, by fixing their orbital periods (and semi-major axes) in accordance with the mean photometric periods. The mean anomalies and arguments of periastrons are random in the $[0,360^{\circ}]$ range, and the eccentricities are randomly sampled from $[0,0.1]$. In order to check the stability of these synthetic systems, we used the Lyapunov-based fast indicator MEGNO (Cincotta & Simó, 2000; Goździewski et al., 2008).

We present the results in Fig. 1 for the probability of picking up a stable system in the ($m_{\rm 1}$,$m_{\rm 2}$)-plane, for 100 sampled configurations with a fixed pair of masses. Clearly, stable systems are possible unless the masses are larger than $2$–$3$ mJup. Beyond this limit, the probability of guessing a stable configuration sharply decreases. Moreover, the masses in stable systems are well below the planetary threshold of $\sim 14$ mJup (e.g., Spiegel et al., 2011).

4 The best-fitting TTV models

We applied the same orbital model and the TTV model optimisation as in our earlier papers (e.g. Goździewski et al., 2016). We assumed a co-planar system, and the model parameters ${\boldsymbol{p}}=\{m_{i},P_{i},x_{i}\equiv e_{i}\cos\varpi_{i},y_{i}\equiv e_{i}\sin\varpi_{i},T_{i}\}$ for $i=1,2$, where $P_{i}$, $e_{i}$, $\varpi_{i}$ and $T_{i}$ stand for the orbital period, eccentricity, longitude of pericenter and the moment of the first transit, respectively, w.r.t. the initial epoch of $T_{0}=\mbox{BKJD}-139$ d.

We performed a preliminary optimisation of the likelihood function $\cal L({\boldsymbol{p}}$) with the evolutionary algorithms (GEA), but, as anticipated, constraining the eccentricities is difficult due to the mass-eccentricity degeneration (e.g., Hadden & Lithwick, 2014; Deck & Agol, 2015; Jontof-Hutter et al., 2016). Therefore, we restrict eccentricities of the two planets in the GEA search, in order to avoid the “over-fitting” of the TTVs (e.g., Migaszewski et al., 2017; Migaszewski & Goździewski, 2018; MacDonald et al., 2016). The limit $e_{i}<0.05$ ($i=1,2$) is typical for Earth-like multiple Kepler planets near to or involved in MMR (e.g., Kane et al., 2012; Kipping, 2014; Xie et al., 2016; Jontof-Hutter et al., 2016; Shabram et al., 2016).

In Fig. 2, we project the best-fitting GEA solutions yielding $\chi^{2}<1.15$ on the $(e_{\rm 1}+e_{\rm 2},\Delta\varpi)$-plane of the osculating elements. For these solutions, we also computed the amplitude of the secular angle $\Delta\varpi=\varpi_{2}-\varpi_{1}$, as well as of the critical angles of the 3:2 MMR, $\phi_{\rm 3:2,1}=2\lambda_{\rm 1}-3\lambda_{\rm 2}+\varpi_{\rm 1}$, $\phi_{\rm 3:2,2}=2\lambda_{\rm 1}-3\lambda_{\rm 2}+\varpi_{\rm 2}$, where $\lambda_{i}$ is the mean longitude of the $i-th$ planet, $\varpi_{i}$ its longitude of periastron, and the indexes $1,2$ are for the inner and outer planet, respectively.

The distribution of $2.5\times 10^{5}$ GEA models illustrated in Fig. 2 clusters around $\Delta\varpi\simeq 180^{\circ}$, and is qualitatively different from that one of Kepler-29, for which most of the eccentricity-unconstrained GEA models are characterised by $\Delta\varpi\sim 0^{\circ}$ (Migaszewski et al., 2017).

For the set of models illustrated in Fig. 2, we calculated the amplitude of the critical and secular angles for a fixed integration time equal to $1000$ yrs, equivalent to $\sim 2\times 10^{5}$ dynamical (outer) periods. We selected four representative solutions marked in Fig. 2. We listed their Keplerian osculating orbital elements in Tab. 2. The GEA-I solution in Tab. 2 represents low-eccentricity anti-aligned configurations, the GEA-II model is an example of moderate-eccentricity, aligned ($\Delta\varpi=0^{\circ}$) systems. The GEA-III model is a representative configuration for anti-aligned solutions with moderate eccentricity, while the GEA-IV model is a representative of the solutions with eccentricities close to the upper limit imposed in the GEA search. The anti-aligned configurations are the most frequent solutions inferred from the GEA search. We found that almost all solutions exhibit one or two critical angles librating, while only 172 models (marked in red) exhibit circulating critical angles. We distinguish between aligned models ($\Delta\varpi\simeq 0^{\circ})$ for which only one critical angle librate, while both angles librate in anti-aligned solutions ($\Delta\varpi\simeq 180^{\circ})$. Remarkably, such behaviour is independent of eccentricities and masses. Moreover, the circulating critical angles remain highly coherent. Furthermore, we found that such coherent circulations may be related to the proximity of the solutions to the separatrices of the 3:2 MMR, as will be demonstrated in Sect. 5. Therefore, the GEA experiment results indicate that the system can be considered as resonant in statistical sense.

Due to the short observational window, it is not possible to constrain the architecture of the system by inspecting the GEA models statistics. However, having the best-fitting GEA solutions as the initial guess, we used the MCMC sampling for characterising these solutions through imposing Bayesian priors on the TTV model parameters. In particular, we set Gaussian priors ${\cal N}(\mu,\sigma)$ for the $(x_{i},y_{i})$-variables with $\mu_{i}=0$ and the same $\sigma_{x_{i},y_{i}}$ for both planets. In order to assess proper values of these priors, we performed the MCMC sampling for $\sigma_{x_{i},y_{i}}\in[0.001,0.12]$, with small steps, following a similar strategy as in Migaszewski & Goździewski (2018). The priors were set uniform for all other parameters. We iterated 1024 emcee walkers around selected GEA models for up to 256,000 samples each, aiming to keep the acceptance rate between 0.2 and 0.5.

In Fig. 3, we show the posterior distributions in the $(m_{\rm 1},m_{\rm 2})$- and $(\varpi_{2},\Delta\varpi)$-plane, $\Delta\varpi\equiv\varpi_{2}-\varpi_{1}$, for two representative $\sigma_{x_{i},y_{i}}=$ 0.0066 and 0.12, respectively. For eccentricity priors smaller than the critical one, $\sigma^{C}_{x_{i},y_{i}}\simeq 0.03$ (tentatively), the posterior is single-modal (the left column of Fig. 3), with $\Delta\varpi\simeq 180^{\circ}$. We found the posterior more and more asymmetric in $\Delta\varpi$ for $\sigma_{x_{i},y_{i}}>\sigma^{C}_{x_{i},y_{i}}$, and with two clear local extrema, $\Delta\varpi\simeq 180^{\circ},0^{\circ}$ for $\sigma_{x_{i},y_{i}}>0.06$. Close to $\sigma^{C}_{x_{i},y_{i}}$, the $\Delta\varpi=180^{\circ}$ mode bifurcates, and a second mode $\Delta\varpi=0^{\circ}$ (aligned orbits) emerges. As the best-fitting TTV model, we report in Tab. 3 a low-eccentric ($e_{i}\simeq 0.01$), $\Delta\varpi=180^{\circ}$ solution for $\sigma_{x_{i},y_{i}}=0.0033$. It yields similarly small $\chi_{\nu}^{2}\simeq 1.1$ as the best GEA models, $\chi^{2}_{\nu}\sim 1$, close to a local minimum, when compared with models for $\sigma_{x_{i},y_{i}}=0.001$ ($\chi_{\nu}^{2}=1.22$) and $\sigma_{x_{i},y_{i}}=0.0066$ ($\chi_{\nu}^{2}=1.27$).

Plots in the right column of Fig. 3 illustrate the median values of eccentricities and masses derived from the posterior samples, for a number of runs with $\sigma_{x_{i},y_{i}}\in[0.001,0.12]$. There is a strong correlation of the median eccentricity with the priors, while mass estimates seem to be clustered, yet in two different regions, relative to $\sigma^{C}_{x_{i},y_{i}}\simeq 0.03$. Beyond that value, $(x_{i},x_{j})$ and/or $(y_{i},y_{j})$ are correlated, and $\Delta\varpi=180^{\circ},0^{\circ}$ modes appear clearly for $\sigma_{x_{i},y_{i}}>0.06$, a value likely dependent on the MCMC sampling strategy.

This experiment shows that it is not possible to distinguish between solutions exhibiting the two $\Delta\varpi$ modes and different eccentricities with only TTV observations. Also, the mass-eccentricity degeneracy (e.g., Hadden & Lithwick, 2014; Jontof-Hutter et al., 2016) cannot be fully removed with a “reasonable” selection of the eccentricity priors. The masses and eccentricities are globally weakly constrained, mostly due to the two-modal $\Delta\varpi=0^{\circ},180^{\circ}$ posterior. Additional constraints, like a particular type of periodic configurations (Migaszewski & Goździewski, 2018), or flowing from the migration history of the system, could be helpful for resolving this issue. Indeed, as described in detail in Sect. 6, for reasonable disk decay and migration time-scales $(\tau_{a},\tau_{e})$, we obtain an agreement between the numerical simulations and our best-fitting solutions, having synthetic systems with $\Delta\varpi=180^{\circ}$ and moderate eccentricities.

5 The 3:2 MMR dynamics

Although both the GEA and MCMC-sampling experiments make it not possible to constrain eccentricities without additional assumptions, a striking feature of the TTV models is their clustering in the anti-aligned libration mode $\Delta\varpi=180^{\circ}$. This feature is expected as a natural outcome of inward and convergent migration of two-planet systems (Lee & Peale, 2002; Batygin & Morbidelli, 2013). Combinations of low—moderate eccentricity and aligned—anti-aligned configurations, which fit the TTV observations are also possible. Also the best-fitting, low ${\chi^{2}_{\nu}}\simeq 1$ systems exhibit period ratios close to 3:2. In order to assess whether these systems are dynamically resonant, in spite of possible large-amplitude libration of the critical angles, or only their coherence (such as shown in Fig. 4), we performed additional numerical experiments regarding four classes of configurations, illustrated in Fig. 2 and listed in Tab. 2.

A clear libration of the critical angles may not be the decisive factor for identifying the MMR dynamics (e.g., Migaszewski et al., 2017; Petrovich et al., 2013; Delisle et al., 2012; Henrard & Lamaitre, 1983). Here, we follow the strictly dynamical understanding of the resonance (MMR), as the 3:2 commensurability region in the parameter space, in which the proper (fundamental) frequencies are closely commensurate and dynamically bordered by separatrices (boundaries between different modes of orbital evolution). The presence of separatrices in multi-dimensional planetary systems usually leads to zones of chaotic motions. In order to detect such structures in the orbital elements space, we computed dynamical maps in terms of the Maximal Lyapunov Characteristic Exponent (MLCE) expressed by the MEGNO indicator (Cincotta & Simó, 2000). We also computed the proper mean motions $n_{i}$, $i=1,2$ (fundamental frequencies associated with the orbital periods) and their ratio ($f_{2}/f_{1}\equiv n_{2}/n_{1}$).

5.1 Numerical mapping of the 3:2 MMR

A proper parametrization of the dynamical maps is required for detecting the resonance structure. Besides the common $(a_{i},e_{i})$–representation, which crosses all MMRs (Laskar & Robutel, 2001), we considered the so called representative plane of initial conditions (Michtchenko & Ferraz-Mello, 2001; Callegari et al., 2006). The equations of motion of the system are governed by the Hamiltonian expressed in Poincaré coordinates (e.g., Michtchenko & Ferraz-Mello, 2001),

[TABLE]

where $\mu_{i}=k^{2}({m_{\star}}+m_{i})$, $\beta_{i}=k^{2}m_{\star}m_{i}/\mu_{i}$. Besides the total energy integral ${\cal H}$, the total angular momentum $C=G_{1}+G_{2}$ is preserved, where $G_{i}=L_{i}\sqrt{1-e_{i}^{2}}$ and $L_{i}=\beta_{i}\sqrt{\mu_{i}a_{i}}$. Moreover, the system Eq. 1 averaged near a particular MMR $(p+q)/p$, where $p,q$ are integers (here $p=2$, $q=1$), exhibits a particular integral $K=(p+q)L_{1}+pL_{2}$ that bounds variability of semi-major axes in resonant motion.

The $C$ and $K$ integrals depend on four elements $e_{1},e_{2},a_{1},a_{2}$, and if their values are fixed, two linear equations may be solved against the two remaining variables. Since the resonant and non-resonant dynamics are governed by librations of the critical angles $\phi_{3:2,1}$ and $\phi_{3:2,2}$ around $0^{\circ}$ or $180^{\circ}$, or circulations, still involving these two critical values, the structure of the phase space restricted to the same $C$ and $K$ levels may be illustrated in a plane of eccentricities $\Sigma(e_{1},e_{2})$. This plane is composed of four quadrants with $\phi_{1},\phi_{2}$ fixed at $0^{\circ}$ or $180^{\circ}$, encoded by variables $x\equiv e_{1}\cos\phi_{1}$ and $y\equiv e_{2}\cos\phi_{2}$. Here, we fix $\varpi_{1}=0^{\circ},{\cal M}_{1}=0^{\circ}$, and then the remaining angles are $(\varpi_{2}=0^{\circ},{\cal M}_{2}=0^{\circ})$, $(\varpi_{2}=180^{\circ},{\cal M}_{2}=180^{\circ})$, $(\varpi_{2}=0^{\circ},{\cal M}_{2}=180^{\circ})$, and $(\varpi_{2}=180^{\circ},{\cal M}_{2}=0^{\circ})$, respectively. These pairs of angles define quadrants I, II, III and IV of the representative plane $\Sigma$.

We derived the proper mean motions with the frequency modified Fourier transform (Frequency Modified Fourier Transform aka Numerical Analysis of Fundamental Frequencies, NAFF, Laskar, 1990, 1993; Šidlichovský & Nesvorný, 1996) of the time series $\left\{a_{i}(t_{k})\exp[\mbox{i}\psi_{i}(t_{k})]\right\}$, where $\psi_{i}$ are appropriate angles forming the conjugate action-angle pairs with the osculating semi-major axes $a_{i}(t)$ (equivalent to rescaled actions $L_{i}$), sampled at discrete moments $t_{k}$, $k=1,2,3,\ldots,2^{K}$ (here $K=18$), inferred from the Poincaré coordinates (e.g., Morbidelli, 2002; Goździewski et al., 2008). As the conjugate angles, we may choose the mean longitudes $\psi_{i}\equiv\lambda_{i}(t)=\mathcal{M}_{i}(t)+\varpi_{i}(t)$ or the mean anomalies $\psi_{i}\equiv\mathcal{M}_{i}(t)$. The meaning of the NAFF-derived fundamental frequencies $n_{i}$ is then subtly different. In the $\psi_{i}=\lambda_{i}(t)$ settings, the proper frequencies $n_{i}$ are related to the inertial frame, while for $\psi_{i}=\mathcal{M}_{i}(t)$, they represent the orbital motion in particular rotating reference frames related to the rotations of the longitudes of pericenter of each orbit. If the orbital configuration of both planets is periodic, as explained below, then their apsides rotate with the same frequency, thus there is a common reference frame corresponding to the rotation of the planetary system as a whole.

We found the conjugate angles distinction as important for detecting regions close to the periodic orbits associated with the 3:2 MMR, or rather the 3:2 commensurability. The periodic orbits, meaning a repetitive, relative configuration of planets in the rotating reference frame with one of the planets, are associated with centers of the mean motion resonances (Hadjidemetriou, 2006). These solutions, when related to symmetric periodic orbits, may be characterised by $\Delta\varpi=0^{\circ},180^{\circ}$ and $\phi_{3:2,1}=0^{\circ},180^{\circ}$ in the averaged system (also, Voyatzis, 2008) or low-amplitude librations around these values in the original (full) system. Regarding the 3:2 MMR, these conditions correspond to a constant value of $2\mathcal{M}_{1}-3\mathcal{M}_{2}\simeq 0^{\circ},180^{\circ}$ for the periodic orbits. Therefore, by using the NAFF we may check whether $f_{1}/f_{2}\simeq 3/2$ in selected orbital parameter planes, thus detecting the regions associated with the resonant periodic orbits, besides illustrating the width and structure of the resonance.

5.2 The 3:2 MMR structure for selected models

In the top row of Fig. 4, we over-plotted the observed TTVs for both planets and synthetic TTV signals obtained from the best-fitting MCMC solution (Tab. 3). In the middle row, we show two dynamical maps in the $(a_{i},e_{i})$–plane for this model illustrating the 3:2 commensurability region. The left-hand map shows a deviation of the proper mean motions ratio relative to the exact 3:2 MMR value, computed by using the mean longitudes (from herein, $\lambda$-NAFF). In this map, the 3:2 commensurability region may be indicated by small values of $|f_{1}/f_{2}-3/2|$. We obtained the right-hand map with the FMFT analysis of time series involving the mean anomalies ($\mathcal{M}$-NAFF). In this plot, the family of symmetric periodic orbits is associated with the 3:2 MMR and close solutions appear as dark-blue and dark regions in the $(a_{1},e_{2})$–plane. Note that this plane is different from the one presented in the left panel. In the black regions in this map, the $f_{2}/f_{1}$ deviation from 3/2 may be as small as $10^{-9}$. We note that due to a complex projection of the multi-dimensional elements space onto the selected two–dimensional plane, the family of periodic orbits generally does not appear as a simple curve or an isolated region. We choose the $(a_{1},e_{2})$–plane after some experiments with the aim of crossing the phase-space of the 3:2 MMR in a representative way (yet the choice of the crossing plane is not unique).

In the bottom row of Fig. 4, we show the evolution of two resonant solutions, the nominal MCMC one (left), and the one with eccentricity for the second planet $\simeq 0.02$, i.e., the nominal value increased by $50\%$ (right). A different proximity to the periodic orbit is related to different behaviour of the critical angles. The NAFF dynamical maps reveal the MMR structure and mark the separatrices. However, while the right separatrix is clearly visible in the whole $e_{1,2}$-range, the left separatrix apparently “diffuses” and vanishes at small eccentricities. In fact, in the small-eccentricity regime the separatrix may not exist (Henrard & Lamaitre, 1983; Delisle et al., 2015), and the $\Sigma$-plane is more useful for resolving the resonance structure.

Figure 5 shows two dynamical maps in the $\Sigma$-plane computed in terms of $|f_{2}/f_{1}-3/2|$. The left-hand panel is for $\lambda$-NAFF. In this plot we also marked unstable solutions detected with the symplectic MEGNO indicator, computed for 72,000 years with a time-step of 0.5 days. They form a loop, with a filled circle inside, which marks a stable periodic orbit associated with the 3:2 MMR. The nominal model (marked with a star symbol) lies outside the loop. The right-hand panel shows this region in terms of the $\mathcal{M}$-NAFF (only the quadrant II is shown). The loop of unstable solutions overlaps with an increase of $|f_{2}/f_{1}-3/2|$, and the libration zone has a sharp border.

In the left column of Fig. 6, for a reference as well as for a better representation of the TTV models, we report two–dimensional $(a,e)$ $\mathcal{M}$-NAFF dynamical maps for the four representative GEA models in Tab. 2. For the anti-aligned, low eccentric (GEA I), moderate eccentric (GEA-III), and high eccentric (GEA-IV) models, respectively, the resonant structures are generally similar. Some differences may be observed between fine structures presented in the maps, which indicate a small dependence on eccentricities. As explained above, the left separatrix does not clearly appear at low eccentricities due to its very narrow width. Large regions of strong dynamical instability appear at moderate eccentricities beyond the resonance borders, yet the 3:2 resonance persists for eccentricities as large as 0.3. All the GEA solutions are found inside the 3:2 MMR dynamical structure bounded by two separatrices in the $(a_{1},e_{2})$– and $(a_{1},e_{1})$–planes. We recall that in the $\mathcal{M}$-NAFF maps, dark-blue zones and strips, with the lowest values of $|f_{1}/f_{2}-3/2|$ correspond to the periodic orbits in the rotating frame, and are associated with the 3:2 MMR centers for the given map coordinates.

In the middle column of Fig. 6, we calculate the amplitude of the two critical angles in one–dimensional $a_{1}$ scans across fixed, nominal eccentricities. As for the NAFF structures, these scans are similar to each other, showing critical angles librating with low amplitude close to the centre of the 3:2 commensurability region. Significant differences appear for the aligned case, in which the first critical angle ($\phi_{\rm 3:2,1}$) circulates in all the 1-dim domain. In the right column of Fig. 6, the time evolution of the critical angles is represented in the ($\phi_{\rm 3:2,1}$,$\phi_{\rm 3:2,2}$)-plane.

In Fig. 7, we show the $\lambda$-NAFF dynamical maps for the four GEA best-fitting solutions, computed at the $\Sigma$-plane. The GEA I solution is topologically similar to the MCMC one (see Fig. 5), being close to the stable periodic orbit related to the 3:2 MMR. The GEA II, GEA III, and GEA IV are topologically similar to each other. All of them are inside the dynamically resonant region determined by the very low values of $|f_{2}/f_{1}-3/2|$. For GEA III (which we chose as a representative for this group of configurations), the MEGNO unstable regions overlap with the $\lambda$-NAFF resonant borders. Interestingly, the location of the GEA III best-fitting solution coincides almost exactly with the position of the periodic orbit.

6 Migration

As we demonstrated above, the TTV fitting cannot give any unique best-fitting configuration. There are observationally permitted systems with eccentricities ranging from very small values up to $\sim 0.1$ and possibly beyond. The relative orientation of the apsidal lines cannot be constrained either, since both aligned and anti-aligned orbits are possible. We have shown also that the MCMC best-fitting configurations depend strongly on the eccentricities priors.

It is well known that convergent migration results in resonant systems with anti-aligned apsides (e.g., Batygin & Morbidelli, 2013). Moreover, if the migration is smooth and acts long enough the final systems are periodic configurations (Migaszewski, 2015). Recently, Migaszewski & Goździewski (2018) showed that anti-aligned orbits can be misinterpreted as aligned ones due to the TTV model degeneracy. This is why we limit our analysis to anti-aligned systems and we will not make any attempts to reconstruct the aligned configurations on the way of migration. However, even with this restriction, a migration-based reconstruction of the KOI-1599 system formation is a non-trivial task due to complex and non-deterministic constraints.

Inward migration of planets due to their tidal interaction with the protoplanetary disc is a widely accepted formation scenario for short-period planets. It is also well known that the convergent migration of two planets results in locking them into mean motion resonance (e.g., Snellgrove et al., 2001; Lee & Peale, 2002). In this section, we study the planetary migration as a way in which the KOI-1599 system was formed. We use the parametric model of migration (Papaloizou & Larwood, 2000; Beaugé et al., 2006; Moore et al., 2013; Voyatzis, 2016), in which the $N$-body astrocentric Newtonian equations of motion of the $i$-th planet are completed with the following acceleration which mimics the planet-disc interaction:

[TABLE]

where ${\boldsymbol{v}}_{i}$ is the $i$-th planet’s astrocentric velocity, ${\boldsymbol{v}}_{c,i}$ is its Keplerian velocity for a circular orbit of radius $r_{i}$ (the astrocentric distance of planet $i$). The time scale of migration is denoted with $\tau_{i}$, while the circularization time scale is given by $\tau_{i}$ divided by a factor $\kappa_{i}$. Within the model, we can choose the parameters freely in order to obtain the observed configuration, without considering particular disc properties.

Since, as shown above, the observational system appears as dynamically resonant, and likely has evolved into the 3:2 MMR via migration, we may a priori consistently impose the TTV signals characteristics, like their amplitudes and periodicities, as well as the period ratio. In this way, we avoid using the best-fitting Keplerian orbital elements as a target for the migration simulations, since they are dependent on additional assumptions (like eccentricity priors).

Before constrain the observables mentioned above, we ask how the TTV signals of a system which was formed on the way of migration should look like. The TTV signals of a periodic configuration should be sinusoidal with a period that equals the so-called super-period (Lithwick & Wu, 2012), which is related to the rotation of the system as a whole, i.e., $T_{{\mbox{\scriptsize s-p}}}=|q/\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}}-(q+1)/\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}}|$ for the (q+1):q resonance, where $\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ and $\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}}$ denote the mean periods of the transit times (TT) series computed over the time interval longer than the super-period. Additionally, the signals are in anti-phase.

If the system is shifted from the periodic configuration, there appears another periodicity in the TTVs, which is related to the resonant modulations of the semi-major axes. Figure 8 shows the Lomb-Scargle periodograms of example systems of masses like in the KOI-1599 representative best-fitting configuration (Tab. 3) and the $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}=1.501$ (close to the value for the best-fitting system). The thick curve in the centre of the plot corresponds to the periodic system. There is only one periodicity $\approx 10000\,$days, i.e., much longer than the observing window. The remaining six curves are for the configurations shifted from periodic, i.e., the eccentricity of the inner planet is being varied between [math] and $0.024$ (for the periodic system $e_{1}\approx 0.012$), while all other parameters are left unchanged. The super-period changes from one system to another, because when changing $e_{1}$, $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ changes as well, even though $P_{2}/P_{1}$ (the osculating period ratio) is the same. But more importantly, an additional periodicity $\sim 3000\,$days appears. The grey area indicates roughly the observed periodicity of the TTVs (see the Sect. 6.1). The amplitude of this additional signal increases with the distance from the periodic configuration. Therefore, we conclude that the KOI-1599 system is shifted from periodic, the observed TTVs are due to the resonant modulations of $a_{1},a_{2}$ and the observing window is much narrower than the super-period.

6.1 TTV signals characteristics

In the next step, we constrain the TTVs characteristics. As we mentioned, the TTV signals can be modelled with two sine (or cosine) functions. In other words, the transit times model $t_{n}$ should consist of three parts:

[TABLE]

where

[TABLE]

is the linear model of transit times, usually used in order to find the observed period $\langle P\rangle_{{\mbox{\scriptsize obs}}}$, i.e., the mean time distance between subsequent transits and the epoch of the first transit, $T$. For a system of two interacting planets those quantities, in general, depend on the width of the observing window. Therefore, the values obtained for the window of $\sim 4$ years (as it is for the Kepler candidates) usually differ from the ones computed over time longer than the longest periodicity in TTV. This means that for a wide enough observing window (in the example illustrated in Fig. 8 the window should be wider than $\sim 10000\,$days) $\langle P\rangle_{{\mbox{\scriptsize obs}}}$ and $T$ stop depending on the observing window width. Further in this section, $\langle P\rangle_{{\mbox{\scriptsize obs}}}$ and $T$ have the sense given above.

The model of the transit times series is completed with two cosine functions

[TABLE]

and

[TABLE]

The first one describes the TTV caused by the rotation of the system as a whole, with a characteristic frequency related to the super-period, $|N_{{\mbox{\scriptsize s-p}}}|\equiv 2\pi/T_{{\mbox{\scriptsize s-p}}}$, i.e., $N_{{\mbox{\scriptsize s-p}}}=(q+1)n_{2}-qn_{1}$, where $n_{1}\equiv 2\pi/\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}}$, $n_{2}\equiv 2\pi/\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}}$. The amplitude of this signal $A_{{\mbox{\scriptsize s-p}}}$ is related to the eccentricity in the following way:

[TABLE]

where $\langle e\rangle$ is the mean value of the eccentricity. The above relation stems from the fact that for the Keplerian motion (for small $e$) the true anomaly can be expressed through the mean anomaly as $\nu\approx\mathcal{M}+2e\sin\mathcal{M}$, therefore the maximal difference between anomalies with respect to the full angle is $(2e)/(2\pi)$, what gives the amplitude with respect to the period. The phase of the signal for the inner planet can be computed when $\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}}$, $\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}}$, $T_{1}$ and $T_{2}$ are known as

[TABLE]

where $\phi_{1}^{{\mbox{\scriptsize(eq)}}}$ is the equilibrium value of the resonant angle $\phi_{1}\equiv q\lambda_{1}-(q+1)\lambda_{2}+\varpi_{1}$, i.e., $\phi_{1}^{{\mbox{\scriptsize(eq)}}}=0$ for $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}>1.5$. The phase for the outer planet $\Phi_{{\mbox{\scriptsize s-p,2}}}=\Phi_{{\mbox{\scriptsize s-p,1}}}+\pi$.

The second signal $\Delta t_{n}^{{\mbox{\scriptsize(res)}}}$ stems from the resonant evolution of the system, i.e., the semi-major axes modulations. Unlike $N_{{\mbox{\scriptsize s-p}}}$ (which depends on $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ only), the frequency of the second signal $N_{{\mbox{\scriptsize res}}}$ depends on the planets’ masses and eccentricities.

The model of the transit times, Eq. 2 may be now used in order to find the best-fitting parameters $\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}}$, $\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}}$, $T_{1}$, $T_{2}$, $A_{{\mbox{\scriptsize s-p,1}}}$, $A_{{\mbox{\scriptsize s-p,2}}}$, $N_{{\mbox{\scriptsize s-p}}}$, $\Phi_{{\mbox{\scriptsize s-p,1}}}$, $\Phi_{{\mbox{\scriptsize s-p,2}}}$, $A_{{\mbox{\scriptsize res,1}}}$, $A_{{\mbox{\scriptsize res,2}}}$, $N_{{\mbox{\scriptsize res}}}$, $\Phi_{{\mbox{\scriptsize res,1}}}$, $\Phi_{{\mbox{\scriptsize res,2}}}$. In general, that would make $14$ free parameters to be found if there were no dependencies between them. However, as already mentioned $N_{{\mbox{\scriptsize s-p}}}$, $\Phi_{{\mbox{\scriptsize s-p,1}}}$ and $\Phi_{{\mbox{\scriptsize s-p,2}}}$ are functions of $\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}},\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}},T_{1},T_{2}$. The amplitudes $A_{{\mbox{\scriptsize s-p,1}}},A_{{\mbox{\scriptsize s-p,2}}}$ depend on $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ when the masses are given. Moreover, the phases of the resonant part of the model are in anti-phase, i.e., $\Phi_{{\mbox{\scriptsize res,2}}}=\Phi_{{\mbox{\scriptsize res,1}}}+\pi$. As a result there are only $8$ free parameters which have to be found. The dependence of $A_{{\mbox{\scriptsize s-p,1}}}$ and $A_{{\mbox{\scriptsize s-p,2}}}$ on $m_{1},m_{2}$ may seem a problem, since the masses are known within $\sim 50\,\%$ uncertainties. However, $N_{{\mbox{\scriptsize s-p}}}$ and the phases are independent of the masses and for $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ very close to $3/2$ the super-period is at least an order of magnitude longer than the observing window while for $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ shifted significantly from $3/2$ the amplitude $A_{{\mbox{\scriptsize s-p}}}$ is below the TTs uncertainties. As a result, the best-fitting free parameters of the model are almost independent of the masses, assuming the super-Earth mass range.

The results of fitting the model given by Eq. 2 to the transit times series of the KOI-1599 system are presented in Fig. 9. There is no minimum of the $\chi^{2}$ function. In order to illustrate the shape of $\chi^{2}$ we choose and fix $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ and $N_{{\mbox{\scriptsize res}}}$ from a grid and find a minimum of $\chi^{2}$ in the space of the remaining parameters. Next, we compute the confidence levels and plot them in Fig. 9a. The (O-C) period, $T_{{\mbox{\scriptsize res}}}\equiv 2\pi/N_{{\mbox{\scriptsize res}}}$ is correlated with $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$, although there is maximal allowed $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}\approx 1.5037$, and thus also the period can be constrained to be in a range of $(\sim 2200,\sim 4250)\,$days. For the representative configurations listed in Tab. 3 and Tab. 2 $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}\lesssim 1.5016$.

We did not check for $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}<1.5$, because the periodic configurations for $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ close to, but lower than $1.5$ are unstable. As we study the migration-induced formation of the KOI-1599 system, $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ could be lower than $1.5$ if i) the migration was divergent and stopped when $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ was just below $1.5$, ii) the convergent migration was very fast so the system passed through the resonance and stopped just below it, or iii) the convergent migration–induced resonance capture was only temporary. The two former situations are in principle possible, but the formation of this particular configuration could have happened only by pure coincidence. The latter case would require $m_{1}<m_{2}$ in order to make the capture over-stable (e.g., Delisle et al., 2015; Xu & Lai, 2016), i.e., the equilibrium should be unstable against migration, and we know from the MCMC TTV modelling that $m_{1}\approx 2m_{2}$.

Figure 9b illustrates the results of the fitting procedure for the amplitudes of the resonant modulations in TTV. They are even more strongly correlated with $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ and the allowed values range from $\sim 3$ to $\sim 18$ hours for the inner planet’s signal and from $\sim 8$ to $\sim 43$ hours for the outer planet. The amplitudes seem to be significantly overestimated when looking at the standard (O-C) diagrams presented in Fig. 4. Nevertheless, the situation becomes clear from Fig. 10. The two top panels show the (O-C) diagrams, but instead of taking $P,T$ values from the standard linear fit to the TTs series, we use $\langle P\rangle_{{\mbox{\scriptsize obs}}}$ and $T$ stemming from the model given by Eq. 2. Because, as mentioned earlier, there is no single best-fitting model, we chose arbitrarily $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ from the allowed range, i.e., $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}\approx 1.50032$ (more precisely, $\log_{10}(\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}-1.5)=-3.5$), for which the best-fitting $T_{{\mbox{\scriptsize res}}}\approx 2560\,$d, which is almost twice as long as the observing window. The super-period for this value of $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ equals $\approx 40000\,$d. The two components of the TTV model are drawn with grey and black curves, respectively. While the top panels show the model and the data in a time interval long enough to encompass the super-period, the bottom panels of Fig. 10 focus at the observing window only. The dashed lines in the bottom panels indicate the standard linear model of TTs. We see that not only the mean periods computed for the observing window differ from $\langle P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ and $\langle P_{2}\rangle_{{\mbox{\scriptsize obs}}}$, but also the TTV amplitudes one could infer from the standard (O-C) diagrams are much smaller than the actual amplitudes of the resonant modulations of the TTV.

6.2 Constraining the migration parameters

In the previous subsection, we put constraints on the observables $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}},T_{{\mbox{\scriptsize res}}},A_{{\mbox{\scriptsize res,1}}}$ and $A_{{\mbox{\scriptsize res,2}}}$. The next step is to find the migration parameters and possibly the planet’s masses which could reconstruct the observed quantities of the TTV data. As we mentioned, in order to fit the model, Eq. 2, to the data, we needed to fix the planets’ masses, because $A_{{\mbox{\scriptsize s-p,1}}},A_{{\mbox{\scriptsize s-p,2}}}$ depend on $m_{1},m_{2}$ and we could not treat the amplitudes as free parameters due to the too narrow observing window. Luckily, the resulting constraints on the TTV signal observables do not depend on the masses in the mass range expected for the KOI-1599 planets. The amplitudes $A_{{\mbox{\scriptsize res,1}}},A_{{\mbox{\scriptsize res,2}}}$, on the other hand, are free parameters of the TTs model, as is the period $T_{{\mbox{\scriptsize res}}}$. All those three quantities depend on the planets masses as well as on the eccentricities. As illustrated in Fig. 8, $A_{{\mbox{\scriptsize res,1}}}$ (and similarly $A_{{\mbox{\scriptsize res,2}}}$) increases when the system shifts away from the periodic configuration. It is known, though, that very slow convergent migration into a first order MMR, like 3:2, results in a configuration that is very close to periodic (Migaszewski, 2015). In such a case one could expect small amplitudes $A_{{\mbox{\scriptsize res,1}}}$ and $A_{{\mbox{\scriptsize res,2}}}$. Therefore, in order to obtain significant amplitudes of the TTV signal, the migration must be fast enough. Not too fast, though, since the system would pass through the resonance.

We expect from this qualitative analysis that $A_{{\mbox{\scriptsize res,1}}}$ and $A_{{\mbox{\scriptsize res,2}}}$ will help us to constrain the migration rate, while the planets’ masses could be constrained with $T_{{\mbox{\scriptsize res}}}$, since smaller masses mean slower resonant modulations of $a_{1},a_{2}$. The forth observable, $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ will be useful to constrain the $\kappa$ parameter since it governs the equilibrium period ratio, i.e., too high values of $\kappa$ would not let the system evolve towards $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ very close to $1.5$, and we know that $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}\lesssim 1.5037$.

After running a series of the migration simulations for different $\tau_{1},\tau_{2},\kappa_{1}=\kappa_{2}=\kappa$ and $m_{1},m_{2}$ we found a very good agreement with the observables for $m_{1}=7.0\,\mbox{M}_{\oplus}$, $m_{2}=3.6\,\mbox{M}_{\oplus}$, $\tau_{1}=10\,$kyr, $\tau_{2}=4\,$kyr and $\kappa=27$. The evolution of the osculating period ratio, the eccentricities as well as the critical angles is shown in subsequent panels of Fig. 11. The period ratio reaches the value of $1.5$ and librates in its vicinity with a decreasing amplitude. The eccentricities increase up to the values close to their equilibria (dashed lines), i.e., the so-called forced eccentricities (e.g., Lithwick & Wu, 2012). Note that the eccentricities as well as the masses are close to the GEA-III model (see Tab. 2). During further migration both $e_{1}$ and $e_{2}$ tend towards the dashed lines. The resonant angle $\phi_{1}$ and the difference between the longitudes of the apsidal lines $\Delta\varpi$ are slightly shifted from their equilibria, respectively [math] and $\pi$, and similarly to the eccentricities, they tend towards them during the migration.

Every $\sim 2.5\,$yr of the migration simulation we checked how would the (O-C) signal look like, by integrating the $N$-body equations of motion (without the dissipative terms) over $\sim 300\,$yr (in order to encompass the super-period), evaluating TTs and computing the observables $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$, $T_{{\mbox{\scriptsize res}}}$, $A_{{\mbox{\scriptsize res,1}}}$ and $A_{{\mbox{\scriptsize res,2}}}$. At each time the system was scaled to the size of KOI-1599. The results of the analysis are presented in Fig. 12. Despite large oscillations of $P_{2}/P_{1}$ after entering the resonance (Fig. 11a), $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ is almost constant – see Fig. 12a, note that the period ratio is presented as $\log_{10}(\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}-3/2)$ – and is kept at the value well inside the observational limits, actually very close to the example illustrated in Fig. 10 (and close to the GEA-III model). Looking at this solution globally, we may notice that it belongs to a family of $(\sim 7,\sim 3.6)$ Earth mass solutions which permit both modes of librations, as illustrated in Fig. 3.

Figure 12b illustrates the evolution of both $T_{{\mbox{\scriptsize res}}}$ and $T_{{\mbox{\scriptsize s-p}}}$. The period of the resonant modulations (the black curve) reaches a value within the observational constraints shortly after entering the resonance and stays unchanged. On the contrary, the amplitudes (shown in Fig. 12c,d) oscillate between high values (at the beginning of the resonance capture they are even higher than the observational ranges) and very low values, lower than the TTs uncertainties. Such large oscillations result from the fact that the migration is being switched-off every time we compute the synthetic TTV, and depending on the phase in which it is switched-off, we obtain higher or lower amplitudes of the resonant modulations. Therefore, one should look at the behaviour of the amplitudes in a mean sense. The ranges in which $A_{{\mbox{\scriptsize res,1}}}$ and $A_{{\mbox{\scriptsize res,2}}}$ oscillate decrease while the system reaches the periodic configuration. Therefore, the consistency between the synthetic system stemming from the migration and the observations is temporary. Nevertheless, the consistency holds for up to $\sim 4\,$kyr after entering the resonance. If the migration, which is very fast at the beginning, slowed down, the system could remain consistent with the TTVs even longer.

The evolution of the observables can be also presented in a different manner, see Fig. 13. The observational constrains, illustrated first in Fig. 9, look here as very narrow, slightly curved areas. Instead of studying the time evolution of the four observables, in Fig. 13 we show the evolution at the planes with $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}$ at the $x$-axis. The correspondence between the synthetic system and the data is clear, especially in Fig. 13a, in which the observables fits the observational constraints.

The general conclusion from this section is that the KOI-1599 system could have been formed on the way of migration. The necessary conditions are that the migration is fast enough, the eccentricity damping should not be too effective and the masses should be in the lower limits of the MCMC TTV analysis. The time-scales of the migration $\tau_{1}=10\,$kyr and $\tau_{2}=4\,$kyr, means that the time-scale of the period ratio variation, $\tau_{X}$ ($\tau_{X}=-X/\dot{X}$, where $X=P_{2}/P_{1}$), equals $10\,$kyr. The initial semi-major axis of the inner planet equals $0.4\,\mbox{au}$. Therefore, at the beginning of the simulation, i.e., shortly before entering the resonance, $\tau_{X}$ is as long as $4\times 10^{4}$ of the inner planet’s orbital period. For the migration slowed down by a factor of $2$, the amplitudes $A_{{\mbox{\scriptsize res}}}$ for the synthetic system fits the observational constrains only marginally. On the other hand, fastening the migration by a factor of $2$ makes the system pass through the resonance.

The $T_{{\mbox{\scriptsize res}}}$ dependence on the planets’ masses constrains them relatively well. In order to fit the observational limits for all the four observables, the masses should equal $7.0$ and $3.6\,\mbox{M}_{\oplus}$ within $\sim 10\,\%$ ranges. The observational constraints of $\kappa$ are also quite strong. For $\kappa\gtrsim 100$ the system only marginally reaches the observational upper limit of $\langle P_{2}/P_{1}\rangle_{{\mbox{\scriptsize obs}}}\sim 1.5037$. For lower $\kappa$ the resonant modulations are being damped slower, and the synthetic system stays consistent with the observations over a longer time interval. The $\kappa$ parameter can be therefore limited only from the top, to be $\lesssim 100$.

7 The density dichotomy of KOI-1599 planets

As we found with the photometric analysis of the LCs, the KOI-1599 Super-Earths have very similar radii of $\simeq 1.9~{}\mbox{R}_{\oplus}$ Surprisingly, our estimates differ from data in Rowe et al. (2015). They report a large, likely incorrect radius $\sim 64\mbox{R}_{\oplus}$ of the inner planet, implying its density $10^{4}$ times smaller than the density of the Earth.

The mass estimates in the TTV model in Tab. 3 have small uncertainties of $\simeq 20\%$. Their values are consistent with predictions from the mass–radius relation $M=10^{C}\times R^{E}$, where $C=0.50\pm 0.03$ and $E=0.64\pm 0.06$ (Mills & Mazeh, 2017).

Having the relatively accurate mass and radii estimates, we may infer the densities of the two planets. Given the similar radii and the inner planet twice as massive as the outer one, the internal compositions of the planets are likely different. We used two-layer theoretical models by Zeng et al. (2016) in order to localise the planets in the mass-radius diagram in Fig. 14.

We report two different internal composition models obtained for two MCMC classes of solutions. For the MCMC-model with $\Delta\varpi=180^{\circ}$ (Tab. 3), the inner planet has a bulk density of ($7.2\pm 0.3$)  g$\cdot$cm*-3*, roughly 1.5 times the density of the Earth. The density of the second planet is $(3.7\pm 0.3$)  g$\cdot$cm*-3*, roughly the density of Mars (i.e. 3.9  g$\cdot$cm*-3*). As for the GEA search, also in the case of MCMC search, we have found different local best-fitting solutions, for the dual-mode solutions ($\Delta\varpi=0^{\circ},180^{\circ}$), the masses (Fig. 3) as well as the derived densities are significantly smaller, (5.6 $\pm$ 0.5)  g$\cdot$cm*-3* and (2.9 $\pm$ 0.6)  g$\cdot$cm*-3*, respectively.

Both classes of the internal structure models predict the inner planet as a Super-Earth with a rocky envelope and a small iron-rich core. The outer planet is found in the region between rocky and icy-giant planets, populated by a new class of objects. They are called “Super-Venus” planets (e.g., Kane et al., 2013) if they are, as KOI-1599.02, close to the host star. Other multi-planet systems with planets having dissimilar densities or very different radii and in a very compact orbital configuration are known and well studied, such as Kepler-29 (Migaszewski et al., 2017), Kepler-30 (Panichi et al., 2018), and Kepler-36 (Carter et al., 2012). Among others, Kepler-36 is the most similar system to KOI-1599, models of the structure of such systems predict planets having water in the form of a thick hot atmosphere, as molecular fluid, or ionic fluid, depending on the planet’s mass, composition and UV flux (e.g., Lopez et al., 2012; Lopez & Fortney, 2013; Nettelmann et al., 2008). Recently, Bodenheimer et al. (2018) propose a new scenario for the formation of Kepler-36, likely applicable also to KOI-1599. By using this new approach, which includes the accretion and dissolution of planetesimals into the gaseous envelope of the planets, the effect of in-situ formation and migration, as well as the influence of the mass loss due to XUV radiation it may be possible to understand the difference in the internal structure of the two planets in the KOI-1599 system.

8 Conclusions

With the dynamical photometry, we aim to validate two transiting, candidate planets attributed to the KOI-1599 Kepler target.

From the Q1-Q17 DR-25 LCs, we derived the complete series of TTVs with a few hours amplitude. They exhibit a clear anti-correlation trend, indicating that the candidate planets mutually interact. Dynamical experiments make it possible to determine the upper limits of the candidate masses, roughly below 3 Jupiter masses that provide dynamically stable configurations. The photometric analysis of the LCs of KOI-1599 reveals that the planets have very similar radii, $\simeq 1.9$ times the Earth’s radius.

The orbital model of the TTVs makes it possible to constrain the masses in the Super-Earth range. The uncertainties may be as small as $\simeq 20\%$. However, we assumed that the orbits are low-eccentric ($e_{i}\simeq 0.01$) and anti-aligned, consistent with predictions of the planetary migration theory (e.g., Batygin & Morbidelli, 2013). Given the similar radii, the internal compositions of the planets may be different.

Remarkably, despite of the eccentricities and $\Delta\varpi=0^{\circ},180^{\circ}$ libration modes, the best-fitting orbital models are dynamically “easily” involved in the 3:2 MMR. Dynamical maps reveal that the best-fitting configurations are inside the separatrix region of the 3:2 MMR, although the critical angles can circulate or librate with amplitudes reaching $360^{\circ}$ (so only are coherent). The apparently natural libration–circulation criterion of the critical angles may lead to an incorrect identification of the MMR.

We determine if a particular configuration of the system is dynamically resonant by comparing its position in the phase space w.r.t. the location of the MMR separatrix. Additionally, by calculating the proper mean motions we reinforce the conclusion that a system may be dynamically resonant, having critical angles circulate.

The results of our statistical and dynamical experiments are suggestive of two real planets. By using a simple model for planet-disk interaction, we tried to reproduce the present orbital configuration of the KOI-1599 two-planet system. We find that for fast migration and slow eccentricity damping (with the migration–to–damping time scale ratio $\kappa\lesssim 100$), the synthetic systems well reproduce the best-fitting configurations characterised by librating critical angles and by low values of planetary masses equal to $7.0\mbox{M}_{\oplus}$ for the inner and $3.6\mbox{M}_{\oplus}$ for the outer planet, respectively. The fast migration ($\tau\sim 10^{4}$ orbital periods) means rapid resonance entrance, a result of which is that the system deviates from the branch of periodic configurations. The deviation is responsible for the observed TTV signal (i.e., TTVs result from the semi-major axes modulation instead of the system rotation). The requirement of $\kappa\lesssim 100$ is necessary so the system could reach an appropriate value of the period ratio, i.e., very close to the nominal value of 3:2 MMR. Additionally, the synthetic systems are dynamically resonant, have anti-aligned orbits and small eccentricities, consistent with the ones obtained by our best-fitting configurations.

9 Acknowledgements

We would like to thank the anonymous referees for their review and fruitful comments that improved the quality of our manuscript. K. G. thanks the staff of the Poznań Supercomputer and Network Centre (PCSS) for the support and CPU resources (grant No. 313). This work has been supported by Polish National Science Centre MAESTRO grant DEC-2012/06/A/ST9/00276.

On-line material

10 Transit time variations (this work)

[TABLE]