The Kepler-19 System: A Thick-envelope Super-Earth with Two Neptune-mass Companions Characterized Using Radial Velocities and Transit Timing Variations

After the discovery of thousands of planets with radii smaller than 2.7 R_⊕ from the NASA Kepler mission (Borucki et al. 2011; Coughlin et al. 2016), a consistent effort has been devoted to understanding the formation scenario and chemical composition of such planets (e.g., Weiss & Marcy 2014; Dressing et al. 2015 and Wolfgang & Lopez 2015). To distinguish between a rocky composition and the presence of a thick envelope of water or volatile elements, the radius derived from the transit depth must be coupled with a precise mass determination (better than 20%), either from radial velocity (RV) measurements or transit timing variations (TTVs). Planets that have been characterized with such a level of precision appear to fall into two populations:  one following an Earth-like composition and a second one with planets larger than 2 Earth radii, requiring a significant fraction of volatiles (e.g., Rogers 2015; Gettel et al. 2016 and López-Morales et al. 2016). Recently, improved mass and radius determinations of known planets have uncovered the existence of super-Earths that fall between these two populations, such as 55 Cancri e (Demory et al. 2016) and Kepler-20b (Buchhave et al. 2016).

We need more planets in that mass regime with precise mass and radius measurements to understand the undergoing physics. For this reason we carried out an RV follow-up of Kepler-19b (hereafter K19b), a planet with a period of 9.287 days and a radius of 2.209 ± 0.048 R_⊕, orbiting a relatively bright (V = 12.1, K = 11.9) solar-type star (T_eff = 5541 ± 60 K, $\mathrm{log}g=4.59\pm 0.10$, [Fe/H] = −0.13 ± 0.06). The planet was detected by Borucki et al. 2011 and subsequently validated by Ballard et al. (2011, hereafter B11) using adaptive optics and speckle imaging to exclude a secondary source in the Kepler light curve, Spitzer observations to verify the achromaticity of the transit and Keck-HIRES high-resolution spectroscopy to rule out the presence of massive, non-planetary perturbers. From BLENDER analysis (Torres et al. 2011) the probability of a false-positive scenario was constrained to less than 1.5 × 10⁻⁴. RVs measured on high-resolution spectra were consistent with a mass of 1.5 M_⊕ (0.5 m s⁻¹) and an activity-induced RV jitter of 4 m s⁻¹, but ultimately they lacked the required precision for a robust determination of the planetary mass, and only an upper limit of 20.3 M_⊕ was set. The existence of an additional planet, Kepler-19c (hereafter K19c), with a period ≲160 days and mass ≲6 M_jup, and further confirmation of the planetary nature of K19b, were inferred by B11 from the presence of TTVs on eight quarters of the K19b light curve.

In this paper we couple high-precision RV measurements obtained with HARPS-N with updated measurements of transit times (T₀) encompassing all 17 quarters of Kepler data to determine the orbital parameters of K19b, K19c, and a third, previously unknown planet in the system, Kepler-19d (hereafter K19d). TTV and RV data sets are analyzed independently, to understand which constraints they can provide to the characterization of the system. Subsequently, a simultaneous TTV and RV fit is performed using dynamical simulations to take into account gravitational interactions between planets. We perform this analysis under the assumption of coplanarity between planets and then investigate the effect of different mutual inclinations on the goodness of the fit. We confirm that only the inner planet is seen transiting the host star. A comparison between the RV obtained from dynamical simulations and that when assuming non-interacting planets is performed. We conclude by describing the role of K19b in understanding the bulk densities of small planets.

We collected 101 spectra using HARPS-N at the Telescopio Nazionale Galileo, in La Palma. The observations spanned over two years, from 2012 June to 2014 November, overlapping the Kepler observations during the first year. Every observation consisted of a 30 minute exposure, with a median signal-to-noise ratio of 37 at 550 nm, corresponding to a RV nominal error of 2.8 m s⁻¹. Given the faintness of the target, observations were gathered with the objAB setting, i.e., the second fiber (fiber B) was observing the sky instead of acquiring a simultaneous thorium–argon (ThAr) lamp spectrum. Several observations demonstrated that the stability of the instrument over 24 hr is within 1 m s⁻¹ (e.g., Cosentino et al. 2014), thus the precision of the measurements was dictated largely by photon noise. Data were reduced using the standard data reduction software (DRS) using a G8 flux template (the closest available one to the spectral type of the target) to correct for variations in the flux distribution as a function of the wavelength, and a G2 binary mask was used to perform the cross-correlation (Baranne et al. 1996; Pepe et al. 2002). The resulting RV data, with their formal 1-σ uncertainties, the FWHM of the cross-correlation function (CCF) and its contrast (i.e., the depth normalized to the continuum), the bisector inverse span (BIS), and the $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$ activity index, are listed in Table 1.

Table 1.  HARPS-N Radial Velocities and Ancillary Measurements of Kepler-19

BJDUTC	RV	σRV	FWHM	Contrast	BIS	$\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$	${\sigma }_{{\mathrm{logR}}_{\mathrm{HK}}^{\prime }}$	Moon flag
(days)	(m s−1)	(m s−1)	(km s−1)		(m s−1)	(dex)	(dex)
2456100.608	−10608.31	2.31	6.745	48.83	−41.90	−4.975	0.039	0
2456100.629	−10610.99	2.09	6.735	48.91	−40.73	−5.039	0.031	0
2456101.606	−10613.07	3.56	6.732	48.70	−55.78	−4.961	0.060	0

Note. This table presents the epoch of the observations, the RV with associated noise, the FWHM and contrast of the CCF, the inverse bisector span, the $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$ activity index with the associated error, and a flag indicating if the data have been contaminated by the Moon (1) or not (0); see Section 2.1.

Download table as:  ASCIITypeset image

2.1. Effect of Moon Illumination

Kepler-19 was initially observed in long-cadence (LC) mode during quarters 0–2, and then in short-cadence (SC) mode from quarter 3 until the end of the mission in 2013 (at quarter 17). At the time of publication, B11 had at their disposal only the first 8 quarters. In an effort to make use of any additional information coming from Kepler photometry, we redetermined the transit times for all the quarters now available, in both LC and SC light curves. Quarters already analyzed by B11 were examined as well, to validate our T₀ determination and to provide a homogeneous set of measurements.

Transit identification was performed by propagating the linear ephemeris of B11, with the inclusion of 3 hr of pre-ingress and post-egress around the expected transit time.

For each transit time, we first detrended the transit light curve with a polynomial between the 1st and 10th degrees, with the best-fit degree chosen according to the Bayesian information criterion (BIC). Then, we determined a new T₀ guess with an automatic selection among different search algorithms,¹⁹ fitting the Mandel & Agol (2002) transit model, implemented in PyTransit ²⁰ (Parviainen 2015), fixing all other parameters to the literature value. Finally, the T₀s were refined using the JKTEBOP program (Southworth et al. 2004) and the associated errors were determined with a classical bootstrap approach.

Transit times from LC and SC light curves were matched together, keeping the SC measurements when available. Transit times are reported in Table 2. A comparison with B11 measurements of the observed minus predicted time of transit (O − C), using their linear ephemeris for both data sets, is shown in Figure 1. The scatter of the residuals, well within the error bars, shows that the methodologies are perfectly consistent, i.e., that we are limited by photon noise, data sampling, and/or unknown systematics rather than the exact procedure followed to measure the transit times.

Zoom In Zoom Out Reset image size

Table 2.  Transit Times of Kepler-19 from Q0-Q17

Transit Number	T0 (BJDUTC)	${\sigma }_{{T}_{0}}$ (days)
0	2454959.7074	0.0014
1	2454968.9935	0.0023
2	2454978.2801	0.0020

Download table as:  ASCIITypeset image

Due to an error in the Kepler archiving system, at the time of B11's publication the time stamps of all the Kepler light curves were reported in the Coordinated universal Time system (UTC) instead of the Barycentric Dynamical Time system (TDB).²¹ While this error was did not affect the internal consistency of the B11 analysis, it must be taken into account when comparing time series with timing accuracy better than a few minutes. We corrected for this error before comparing B11 data with our new T₀ measurements.

Atmospheric stellar parameters of Kepler-19 were determined in B11 using Spectroscopy Made Easy (SME, Valenti & Fischer 2005). Since this method may suffer from correlation between derived parameters (Torres et al. 2012), and since we have several high-resolution spectra from HARPS-N at our disposal, we decided to carry out an independent determination with an alternative approach, i.e., equivalent width measurements of individual spectral lines instead of fitting of the whole spectrum. We used all the spectra free from sky contamination to obtain a coadded spectrum with an average S/N of 350.

Stellar atmospheric parameters were determined using the classical line-of-growth approach. For this purpose we used the 2014 version of the line analysis and synthetic code MOOG ²² (Sneden 1973), which works under the assumption of local thermodynamic equilibrium, and we used the ATLAS9 grid of stellar atmosphere models from Castelli & Kurucz (2004), with the new opacity distribution functions and no convective overshooting. Equivalent width measurements were carried out with the code ARESv2 ²³ (Sousa et al. 2015) coupled with the updated linelist of Malavolta et al. (2016), where the oscillator strength of the atomic lines has been modified to correctly take into account the chemical abundances from Asplund et al. (2009).

Temperature and microturbulent velocity were determined by minimizing the trend of iron abundances from individual lines with respect to excitation potential and reduced equivalent width, respectively, while the gravity $\mathrm{log}g$ was adjusted by imposing the same average abundance from neutral and ionized iron lines. For a detailed description of the procedure for the atmospheric parameters and associated errors we refer the reader to Dumusque et al. (2014). The derived atmospheric parameters are summarized in Table 3.

Table 3.  Astrophysical Parameters of the Star

Parameter	B11	This Work
Teff (K)	5541 ± 60	5544 ± 20
$\mathrm{log}g$	4.59 ± 0.10	4.51 ± 0.03
ξt (km s−1)	⋯	0.88 ± 0.05
[Fe/H]	−0.13 ± 0.06	−0.08 ± 0.02
M⋆ (M⊙)	0.936 ± 0.040	⋯
R⋆ (R⊙)	0.859 ± 0.018	⋯
Age (Gyr)	1.9 ± 1.7	⋯
$\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$	−4.95 ± 0.05	−5.00 ± 0.04

Download table as:  ASCIITypeset image

Our stellar atmospheric parameters agree within the uncertainties with the ones determined by B11, including the surface gravity that is usually the parameter most difficult to derive, and there is only a difference of 3 K in T_eff despite the use of two complementary approaches and independent data sets. For this reason we adopted their determination for the mass and radius of the star and the physical radius of K19b based on light curve analysis.

Considerable effort has recently been devoted to analyzing the effect of stellar activity on RV measurements and T₀ determination (Mazeh et al. 2015; Ioannidis et al. 2016). The HARPS-N DRS automatically delivers several diagnostics for activity such as the FWHM of the CCF, the bisector inverse span, and the $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$ index, while several other chromospheric indexes such as the Hα index (Gomes da Silva et al. 2011; Robertson et al. 2013) can be determined from the spectra themselves.²⁴ In Figure 2 analyses of BIS and $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$ are reported as representative of all the indexes. For each index we checked the presence of any correlation with time, either by visual inspection (upper panels of Figure 2) or with the Generalized Lomb–Scargle (GLS) periodogram (Zechmeister & Kürster 2009), where the 1% and 0.1% false alarm probabilities (FAPs) have been computed with a bootstrap approach (middle panels of Figure 2). Finally, the presence of any correlation with RVs was verified by calculating the Spearman rank correlation coefficient ρ, the slope of the linear fit m with its error, and the p-value using the weighted least-square regression²⁵ (lower panels of Figure 2). We omitted the FWHM from the analysis since a few changes to the spectrograph focus during the first year of observations modified the instrumental profile,  and hence the measured FWHM (however, the measured RVs were unaffected.

Zoom In Zoom Out Reset image size

The absence of significant peaks in the periodogram of the indexes under analysis around the expected rotational periods of 32 days (from B11 following Noyes et al. 1984), 34 ± 6 days, and 36 ± 3 days (from HARPS-N $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$, following Noyes et al. 1984 and Mamajek & Hillenbrand 2008, respectively), and the lack of statistically significant correlations between the indexes and RVs, confirmed the low activity level of the star that was already deduced by B11, and is consistent with $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$ = −5.00 ± 0.04.²⁶

We also searched for stellar variability on the most recent Kepler Pre-search Data Conditioning (PDC) light curve, which presents several improvements that correct instrumental trends (and thus is better suited to search for activity modulation) with respect to the light curve available at the time of B11's publication. Most of the time the star is photometrically quiet, while in some parts of the light curve a clear signal, likely due to stellar activity, is detected. We applied the autocorrelation function technique over these portions of light curve and we estimated a rotational period of 30 days. This signal is characterized by a short timescale of decay (a few rotational periods) and a rapid loss in coherence. We comment on the impact on the RV in Section 7.

Dynamical analysis of the system was performed with TRADES ²⁷ (TRAnsits and Dynamics of Exoplanetary Systems, Borsato et al. 2014), an N-body integrator with the capability of fitting RVs and T₀s simultaneously to determine the orbital parameters of the system through χ² minimization.

Since only the innermost planet is transiting, it is extremely difficult to constrain the orbital parameters of all the planets in the system by TTV alone. In fact, attempts to fit the T₀s with a two-planet model resulted in a strong degeneracy between the mass and the period of the second planet, i.e., O–C diagrams with similar shapes could be produced by jointly increasing the mass and the period of K19c. The amplitude and shape of the T₀s, however, can still give us upper limits on the mass and period of the non-transiting planet if we compare the outcome of the dynamical simulations with the maximum mass compatible with the observed semi-amplitude of the RVs.

We proceeded as follows. We used TRADES to perform a fit of the T₀s by assuming a two-planet model and fixing the mass of K19b to a grid of values between 2.5 and 20 M_⊕ and a spacing of 2.5 M_⊕, with its period already measured from the Kepler observations. To each point on this grid we assigned several values for the mass of K19c randomly selected between 5 and 250 M_⊕. For each combination of K19b and K19c masses, the other orbital parameters of the system were left free to vary and their best-fitting values were determined by χ² minimization through the Levenberg–Marquardt algorithm. We discarded those solutions with at least one planet having an eccentricity greater than 0.3, assuming that interacting planets meeting this condition are likely to be unstable.

In Figure 3 we show the results obtained for the mass of K19c as a function of its period. The expected RV semi-amplitudes of K19c as a function of mass and period are superimposed. We note that many of the orbital configurations reported in the plot could be unstable, since dynamical stability was not yet checked  at this stage (stability analysis is introduced in Section 8). We can attempt to estimate a lower limit to the periods and masses of the non-transiting planets according to the observed semi-amplitudes of the RVs by taking advantage of the correlation between the mass and period of K19c. Our RVs have a peak-to-peak variation of 23 m s⁻¹, so if we take into account the additional signal of K19b (a few m s⁻¹ in the case of a Neptune-like density), K19c amplitude should necessarily lie below the K = 10 m s⁻¹ line. This fact suggests that a short period (≲50 days) for K19c should be expected, while nothing can be said regarding K19b from TTV alone. We remind the reader that this analysis cannot be considered conclusive due to the reduced number of points and their scatter around the best linear fit.

Zoom In Zoom Out Reset image size

While only one planet is transiting Kepler-19, the presence of at least one additional planet was inferred by the presence of TTVs. To reveal such planets, we first performed a frequentist analysis on the RV data set by computing the GLS periodogram and iterating over the residuals until no significant periodicity was present. This analysis revealed two signals at 28.6 and 62.3 days with FAPs lower than 1%, while the signal corresponding to K19b at 9.3 days was barely detected (Figure 4).

Zoom In Zoom Out Reset image size

The signal at 28 days is very close to the second-order 3:1 mean motion resonance (MMR), which is among the possible configurations listed by B11 as the cause of the TTVs of K19b. While the FAP of this signal is below the traditional threshold for RV planet detection claims (≃10⁻³), we can compare it with the probability of observing K19c near another MMR resonance. In order to do so, we scrambled once again the RV observations and calculated the fraction of periodograms that had a stronger peak with respect to the untouched data set, in a frequency range of 5% around the interior and exterior 1:2, 2:3, and 3:4 first-order resonances; the 1:3 and 3:5 second-order resonances; and the 1:4 and 1:5 higher-order resonances. The FAP of the signal at 28 days computed in this way decreases to 8 × 10⁻⁴, i.e., it is unlikely that the signal at 28 days is a spurious signal caused by a planet in another MMR configuration.

The signal at ≃28.6 days is consistent with the rotational period of the star obtained from the active regions on the Kepler light curve, which, however, have the characteristic of rapidly decaying and reappearing later at different phases and intensities. We checked if the RV signal had the same properties by performing a jacknife analysis by splitting the RV data set into two parts (at BJD = 2456700 days) and determining the phase and amplitude of the signal in each data set.²⁸ We obtained K = 2.6 ± 1.1 m s⁻¹, ϕ = 3.0 ± 0.3 rad for the first half of the data set, and K = 3.0 ± 0.6 m s⁻¹, ϕ = 2.8 ± 0.2 rad for the second half, i.e., the signal is stable over several years, and therefore unlikely to be due to stellar activity. As an additional check, we computed the stacked Bayesian Generalized Lomb–Scargle (BGLS) periodogram (Mortier & Collier Cameron 2017). The S/N is increasing with the square of the number of observations, as expected for a coherent (planetary) signal, although the plot is heavily affected by the poor sampling and the signal of K19b near the 3:1 MMR, as we verified by performing the same analysis on synthetic data sets with the same temporal sampling. The absence of RV modulation due to stellar activity is supported by the analysis performed in Section 5, where no periodicity or correlation with RVs is seen for any of the activity indexes.

The semi-amplitudes of these signals, however, cannot be inferred by the frequentist analysis alone, since an offset in RV between the data taken before and after 2012 September may exist due to the failure of the first CCD of HARPS-N (see Bonomo et al. 2014). The value of this offset cannot be determined a priori, even when observations of RV standard stars are available, since it may depend primarily on the spectral type of the star (i.e., the two CCDs may have a different efficiency as a function of wavelength). Introducing an RV offset as a free parameter is a common procedure for non-overlapping data sets, e.g., see Benatti et al. (2016).

We then performed a tentative fit using the MCMC code PyORBIT ²⁹ (Malavolta et al. 2016), allowing for two different systemic velocity γ of the star and using a three-planet model to fit the data. We attempted two different fits with both circular and Keplerian orbits. For planets c and d we set uniform priors in the ranges of [10, 50] days and [50, 90] days, respectively, while the other parameters were left to vary within their physically meaningful range (e.g., positive-definite RV semi-amplitude). The results are shown in Table 4.

Table 4.  Orbital Parameters for a Three-planet Model for the Kepler-19 System Obtained from RV-only, except for the Period of K19c, which is Constrained by the Kepler Light Curve

	K19b	K19c	K19d
	Circular Orbits
Period (days)	9.28699 ± 10−5	28.61 ± 0.24	63.0 ± 0.3
K (m s−1)	2.3 ± 0.5	1.7 ± 0.8	3.8 ± 0.6
ϕ (deg)	194.4 ± 0.3	185 ± 43	174 ± 6
Mass (M⊕)	7.4 ± 1.7	7.8 ± 3.6	22.5 ± 3.8
	Keplerian Orbits
Period (days)	9.28699 ± 10−5	28.54 ± 0.27	62.9 ± 0.3
K (m s−1)	2.6 ± 0.6	1.95 ± 0.8	4.0 ± 0.7
ϕ (deg)	198 ± 8	192 ± 29	172 ± 8
$\sqrt{e}\cos \omega$	0.0 ± 0.2	0.0 ± 0.4	${0.31}_{-0.31}^{+0.20}$
$\sqrt{e}\sin \omega$	0.0 ± 0.3	0.0 ± 0.3	$-{0.25}_{-0.22}^{0.35}$
Mass (M⊕)	8.0 ± 1.8	8.6 ± 3.5	23.1 ± 3.8
e	≤0.13	≤0.30	≤0.32
ω (deg)	unconstrained	unconstrained	$-{36}_{-35}^{+65}$

Note. The reference time for the orbital elements is T_ref = 2456624.82263024 days.

Download table as:  ASCIITypeset image

The eccentricities and arguments of pericenter for each of the three planets are poorly constrained by the RVs alone, thus affecting the precision on mass of the presumed planets. A combined analysis of RVs and TTVs is then required to unambiguously detect and characterize the planets in the system.

The amplitude of dynamical perturbations between planets is very sensitive to the eccentricity and angular parameters (i.e., the argument of pericenter ω and mean anomaly at reference time ${{ \mathcal M }}_{0}$) of the planetary orbits, which, however, can be only poorly constrained by the RVs, especially for planets in nearly circular orbits. For this reason we expect an overall improvement on the precision of the orbital parameters by simultaneously fitting RVs and TTVs.

Dynamical simulations are extremely time-consuming, and we have to use all the information at our disposal to reduce the extension of the parameter space. We used the results from the RV fit in Section 7 to put a constraint to the ranges of period, mass, and orbital phase of each planet.

We followed an iterative approach to avoid being trapped in a local minimum of the χ². We started 10 separate runs of TRADES with the Particle Swarm Optimization algorithm and loose priors on the periods (±10 days around the expected period of each planet from the RV analysis), masses (between 0 and 40 M_⊕), and eccentricities (e ≤ 0.5), taking into account the limits imposed by photometry and RVs. We checked the stability of the outcome of each run and then we ran TRADES again on a range of parameters that was half the size of that from the previous run and that was centered of the outcome of the previous run with the lowest χ² among those that satisfied the stability requirement. Convergence was considered achieved when all the runs resulted in similar parameters (within 5% from the mean) and similar χ² (10% from the mean).

Following Gladman (1993), during the numerical integration we checked the stability criterion for each pair of planets: ${\rm{\Delta }}=2\sqrt{3}{R}_{{\rm{H}}}(i,j)$, where ${\rm{\Delta }}={a}_{j}-{a}_{i}$ is the semimajor axis difference between the jth and ith planet, and R_H(i, j) is the mutual Hill radius between planets i and j. At the end of each TRADES fit we performed an N-body integration with SymBA (Duncan et al. 1998) and checked the stability of the result with the Frequency Map Analysis tool (Laskar et al. 1992; Laskar 1993a, 1993b) with the prescriptions of Marzari et al. (2002). A system is considered stable if the coefficient of diffusion is lower than 10⁻⁵, in an unstable or chaotic state otherwise.

We used the value measured by B11 for the inclination i of K19b. Several attempts to fit the mutual inclinations of the non-transiting planets along with the other parameters resulted in unstable solutions with high mutual inclinations. We took advantage of the low mutual inclinations of planetary orbital planes inferred from Kepler multi-planet systems (Fabrycky et al. 2014) and the additional information coming from systems characterized with high-precision RVs (Figueira et al. 2012) to impose coplanarity with K19b for the other planets, while the longitude of the ascending node Ω was fixed to zero for all the planets. This assumption allowed us to drastically reduce computational time. The orbital period and mass of the planet, the eccentricity e, the argument of periastron ω, and the mean anomaly at the reference epoch ${ \mathcal M }$ were left as free parameters.

Differing from Borsato et al. (2014), we fitted $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$ instead of $e\cos \omega$ and $e\sin \omega$, and the mean longitude at the reference epoch $\lambda =\omega +{ \mathcal M }+{\rm{\Omega }}$ (where Ω is the longitude of the ascending node, fixed to zero to be unconstrained by the data) instead of ${ \mathcal M }$. We used scaled stellar and planetary masses, as commonly done in TTV analysis (e.g., Nesvorný et al. 2012). An RV offset between the data taken before and after 2012 September was included as a free parameter to take into account the change of the CCD (see Section 7).

We used the solution obtained with the global exploration of the parameter space as a starting point for the Bayesian analysis. For this purpose we expanded TRADES functionalities with the emcee package (Foreman-Mackey et al. 2013), an affine invariant MCMC ensemble sampler, and made it available to the community in the TRADES repository. Following Feigelson & Babu (2012), we calculated the log-likelihood $\mathrm{ln}{ \mathcal L }$ from the χ² using Equation (1), where dof is the degree of the freedom of the problem.

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\mathrm{ln}(2\pi )\displaystyle \frac{\mathrm{dof}}{2}-\displaystyle \frac{\sum \,\mathrm{ln}{\sigma }^{2}}{2}-\displaystyle \frac{{\chi }^{2}}{2}.\end{eqnarray} \tag{ 1 }$$

We tested three different scenarios. The first model assumed that the signal at ≃28 days (the first one being detected in the RV periodogram) is the only planet in the system other than K19b. The second model is still a two-planet model, but here we assumed that the signal at ≃28 days is due to activity (without any effects on the T₀s) and that the system consists of two planets at ≃9.2 and ≃62 days. The third model assumed that all the RV signals have planetary origins. The results are presented in the following subsections. We note that a planet in a strong MMR with K19b could still produce the TTVs while having a RV semi-amplitude below our detection sensitivity. Following B11, the TTVs of K19b could be explained by a planet in [2:1, 3:2, 4:3] MMRs with masses of [4, 2, 1] M_⊕, respectively. The signal of the perturbing planet would have an RV semi-amplitude on the order of [1, 0.6, 0.3] m s⁻¹, i.e., beyond the reach of modern velocimeters given the magnitude of our target. However, this scenario requires an activity origin for the 28-day signals, which is not supported by the analysis of Section 5 and the coherent nature of the RV signal over years of observation, and is in opposition to the short timescale decay of the spots observed in Kepler photometry, as shown in Section 7. For these reasons we can regard this scenario as unlikely.

8.1. Two-planet Model

We tested the two-planet model following the steps described in Section 8, with the only additional constraint being that K19c should have a period lower than 40 days. We ran the MCMC sampler for 50,000 steps, with a number of chains being twice the dimensionality of the problem. We checked the convergence of the chains using the Gelman–Rubin statistic ($\hat{R}\lt 1.03$, Gelman & Rubin 1992; Ford 2006), and we built the posterior distributions with the last 20,000 steps and a thinning factor of 200.

Our results are listed in Table 5. Rather than using the median or the mode, we summarize the outcome of the analysis by selecting from the chains the sample with the nearest χ² to the median of its distribution and selecting parameters within their confidence intervals. Confidence intervals are computed by taking the 15.87th and 84.14th percentiles of the distributions, and are reported in the table with respect to the selected sample.

Table 5.  Orbital Parameters of the Planets in the Kepler-19 System Obtained from TTV+RV MCMC Analysis Using Different Assumptions for the Number of Planets and the Stellar Activity The Reference Time for the Orbital Elements is T_ref = 2456624.82263024 days

Planet	Period (days)	Mass [M⊕]	$\sqrt{e}\cos \omega$	$\sqrt{e}\sin \omega$	λ(deg)	e	ω(deg)	${ \mathcal M }(\deg )$
Two-planet Model, Pc < 40 days
K19b	${9.287108}_{-0.00003}^{+0.00006}$	${8.2}_{-1.4}^{+1.6}$	${0.23}_{-0.02}^{+0.04}$	${0.21}_{-0.05}^{+0.03}$	${188.9}_{-1.5}^{+0.9}$	${0.10}_{-0.01}^{+0.01}$	${42.5}_{-11.6}^{+5.7}$	${146.4}_{-5.1}^{+10.5}$
K19c	${28.723}_{-0.007}^{+0.0005}$	${17.0}_{-2.4}^{+2.4}$	${0.48}_{-0.04}^{+0.01}$	${0.15}_{-0.11}^{+0.02}$	${175.5}_{-7.0}^{+3.0}$	${0.25}_{-0.04}^{+0.01}$	${17.3}_{-12.7}^{+2.9}$	${158.1}_{-1.0}^{+6.8}$
Two-planet Model, Pc > 40 days
K19b	${9.286970}_{-0.00005}^{+0.00009}$	${8.6}_{-1.7}^{+1.9}$	$-{0.09}_{-0.06}^{+0.11}$	${0.64}_{-0.05}^{+0.04}$	${201.6}_{-6.3}^{+3.5}$	${0.42}_{-0.05}^{+0.06}$	${97.9}_{-9.9}^{+5.5}$	${103.7}_{-2.2}^{+3.6}$
K19c	${63.128}_{-0.007}^{+0.010}$	${21.7}_{-3.7}^{+1.8}$	$-{0.11}_{-0.02}^{+0.25}$	$-{0.59}_{-0.02}^{+0.02}$	${167.7}_{-0.1}^{+14.4}$	${0.36}_{-0.02}^{+0.011}$	$-{100.4}_{-1.2}^{+23.9}$	${268.1}_{-10.1}^{+0.2}$
Two-planet Model, GLS-cleaned RV Data Set
K19b	${9.28696}_{-0.00006}^{+0.00006}$	${8.4}_{-1.6}^{+1.3}$	$-{0.08}_{-0.17}^{+0.02}$	${0.74}_{-0.11}^{+0.02}$	${201.1}_{-1.5}^{+10.2}$	${0.551}_{-0.11}^{+0.03}$	${95.9}_{-1.7}^{+14.8}$	${105.2}_{-4.8}^{+0.3}$
K19c	${63.128}_{-0.008}^{+0.011}$	${17.4}_{-2.9}^{+2.9}$	$-{0.07}_{-0.08}^{+0.11}$	$-{0.57}_{-0.03}^{+0.01}$	${169.3}_{-3.3}^{+10.2}$	${0.332}_{-0.006}^{+0.035}$	$-{97.4}_{-7.4}^{+10.9}$	${266.7}_{-3.4}^{+6.26}$
Two-planet and Activity Model
K19b	${9.28696}_{-0.00004}^{+0.00008}$	${7.8}_{-1.9}^{+2.1}$	$-{0.02}_{-0.11}^{+0.19}$	${0.71}_{-0.06}^{+0.09}$	${197.9}_{-11.3}^{+6.4}$	${0.50}_{-0.07}^{+0.17}$	${91.7}_{-14.8}^{+9.2}$	${106.2}_{-1.8}^{+2.6}$
Act	${28.57}_{-0.11}^{+0.08}$	$({2.3}_{-0.3}^{+1.1})$ a	${0.22}_{-0.68}^{+0.31}$	${0.34}_{-0.23}^{+0.31}$	${178.4}_{-23.6}^{+2.0}$	${0.16}_{0.01}^{+0.45}$	${56.4}_{-36.1}^{+84.4}$	${122.0}_{-86.7}^{+31.4}$
K19c	${63.139}_{-0.019}^{+0.003}$	${16.2}_{-2.7}^{+2.4}$	${0.11}_{-0.24}^{+0.02}$	$-{0.58}_{-0.02}^{+0.04}$	${178.3}_{-15.6}^{+1.16}$	${0.35}_{-0.04}^{+0.02}$	$-{79.0}_{-23.6}^{+1.3}$	${257.4}_{-1.7}^{+9.0}$
Three-planet Model
K19b	${9.28716}_{-0.00006}^{+0.00004}$	${8.4}_{-1.5}^{+1.6}$	${0.17}_{-0.03}^{+0.05}$	${0.29}_{-0.06}^{+0.04}$	${190.3}_{-1.9}^{+1.0}$	${0.12}_{-0.02}^{+0.02}$	${59.1}_{-11.7}^{+68.9}$	${131.2}_{-6.1}^{+10.4}$
K19c	${28.731}_{-0.005}^{+0.012}$	${13.1}_{-2.7}^{+2.7}$	${0.42}_{-0.03}^{+0.04}$	${0.19}_{-0.05}^{+0.10}$	${181.5}_{-4.5}^{+6.0}$	${0.21}_{-0.07}^{+0.05}$	${23.8}_{-6.8}^{+11.1}$	${157.6}_{-6.6}^{+3.1}$
K19d	${62.95}_{-0.30}^{+0.04}$	${22.5}_{-5.6}^{+1.2}$	${0.13}_{-0.38}^{+0.17}$	${0.18}_{-0.19}^{+0.21}$	${170.1}_{-0.3}^{+12.3}$	${0.05}_{-0.01}^{+0.16}$	${55.2}_{-57.4}^{+76.8}$	${114.9}_{-70.5}^{+63.7}$

Notes.

^aSemi-amplitude of the signal in m s⁻¹.

Download table as:  ASCIITypeset image

8.2. Two-planet and Stellar Activity Model

The presence of the activity signal can potentially affect the analysis, but TRADES is not equipped to deal with non-planetary signals in RV data sets. For this reason we decided to remove such signals from the RV time series before starting the global exploration of the parameter space, using as a data set the residuals of the first iteration of the GLS analysis in Section 7. The global solution was used as starting point for an MCMC analysis with PyORBIT, this time using the original RV data set and the activity signal modeled with a Keplerian curve (as similarly done in Pepe et al. 2013). RV and T₀ calculations regarding the two interacting planets were performed by calling the dynamical integrator of TRADES through a FORTRAN90 wrapper. As additional test cases, we ran TRADES, using the two-planet model and the no-activity model on both the original RV data set (imposing a lower limit of 40 days for the period of K19c, to avoid repeating the case analyzed in Section 8.1) and the GLS-cleaned data set. In all cases we repeated the analysis without imposing any constraint on eccentricity, since the eccentricities of both planets were moving toward the upper boundary imposed in Section 8. We followed the same methodology described in Section 8.1 to run the MCMC and extract the results reported in table.

The eccentricity of K19b is extremely high in all the cases we considered. The most likely explanation is that high eccentricities for both planets are required to produce the same T₀s while keeping their masses within the boundaries set by the RVs. Following Burke (2008) we computed the ratio τ of the transit duration of an eccentric orbit with respect to a circular orbit, and we obtained τ = 0.64 ± 0.04 for the original RV data set, τ = 0.58 ± 0.06 for the GLS-cleaned case, and 0.54 ± 0.10 for the planets+activity case. All these values are well below the value of τ = 0.7 that B11 considered as the minimum reasonable value for the transit duration ratio due to the eccentricity of the planet from the analysis of the Kepler light curve. We note that the period-mass combinations we obtained for the outer planet fall in an empty region of Figure 3, which was obtained by selecting those solutions with e < 0.3 for at least one of the planets.

Zoom In Zoom Out Reset image size

8.3. Three-planet Model

Bayesian analyses for the three-planet model were performed using TRADES combined with emcee (see Section 8). From a preliminary analysis we noted that the chains were affected by poor mixing, with a Gelman–Rubin statistic $\hat{R}\simeq 1.3$ (Gelman & Rubin 1992; Ford 2006). We decided to run 100 chains for an extensive number of steps (250,000) to overcome the poor mixing and perform a proper exploration of the parameter space, knowing that we already reside near the global minimum of the χ². After the first 150,000 iterations we did not see any variation of the posterior distributions while increasing the length of the chains, which made us confident of the robustness of our result. We finally built the posterior probability by drawing 40,000 independent samples from the chains, after removing the burn-in part and applying a thinning factor equal to their auto-correlation time (≃100 steps).

Our results are listed in Table 5. The posterior distributions of the fitted parameters are shown in Figure 5. The confidence intervals of the posteriors are computed by taking the 15.87th and 84.14th percentiles of the distributions, and they are reported as error bars around a sample solution selected as in 8.2. From now on, we will use the reported values as a representative solution of the orbital parameters of the planets. In Figure 6 we show the solutions overplotted on T₀s and RVs, with their respective residuals.

Zoom In Zoom Out Reset image size

In Table 6 we compare the outcomes of the different models under examination. The three-planet solution is the favorite one according to the BIC (having a ΔBIC > 10 with respect to all the other models under analysis, Kass & Raftery 1995), and the considerations at the end of Section 8.2 further strengthen this result. Our three-planet solution has a total reduced ${\chi }_{r}^{2}=1.25$ (χ² = 272.7, BIC = 365.5), where the contributions from RV and TTV are respectively 0.37 and 0.88.³⁰ The standard deviation of the residual RVs is 2.9 m s⁻¹, consistent with the average error of the measurements, and while several peaks can be found in the periodogram, none of them reach the 1% false alarm probability threshold (Figure 7). Similarly, the standard deviation of the TTV of 2.2 minutes is consistent with the average error of 2.1 minutes.

Zoom In Zoom Out Reset image size

Table 6.  Statistical Indexes for the Different Models Under Examination

Model	N Parameters	dof	χ2	${\chi }_{\mathrm{red}}^{2}$	$\mathrm{ln}{ \mathcal L }$	BIC
2 planets Pc < 40 days	12	223	${323.8}_{-7.4}^{+5.5}$	${1.45}_{-0.02}^{+0.01}$	${498.6}_{-3.7}^{+1.0}$	${389.4}_{-1.9}^{+5.5}$
2 planets, Pc > 40 days	12	223	${324.0}_{-3.4}^{+4.4}$	${1.45}_{-0.01}^{+0.02}$	${499.5}_{-2.2}^{+1.6}$	${389.5}_{-3.4}^{+4.4}$
2 planets, GLS-cleaned RVs	12	223	${318.8}_{-3.5}^{+4.4}$	${1.43}_{-0.02}^{+0.02}$	${502.1}_{-2.2}^{+1.8}$	${384.3}_{-3.5}^{+4.4}$
2 planets and activity	17	218	${319.5}_{-6.1}^{+8.8}$	${1.47}_{-0.03}^{+0.04}$	${506.3}_{-4.4}^{+3.0}$	${412.3}_{-6.1}^{+8.8}$
3 planets	17	218	${276.7}_{-5.27}^{+6.9}$	${1.27}_{-0.02}^{+0.03}$	${527.7}_{-3.4}^{+2.6}$	${369.5}_{-5.3}^{+6.9}$

Note. The number of parameters in the fit, the degree of freedom (dof), the χ² and its reduced value, the log-likelihood $\mathrm{ln}{ \mathcal L }$, and the Bayesian information criterion (BIC), are reported.

Download table as:  ASCIITypeset image

We determined the masses of the three planets with precisions of ${\sigma }_{{M}_{b}}=1.6$ M_⊕ (error of 19% on planetary mass) for K19b, ${\sigma }_{{M}_{c}}=2.7$ M_⊕ (error of 21%) for K19c, and ${\sigma }_{{M}_{d}}=3.4$ M_⊕ (error of 17%) for K19d. We included in the computation the uncertainty on the stellar mass and the effect of orbital inclinations, assuming as representative distributions for the latter a normal distribution with i_c = i_d = 8994 and dispersions of σ_i = 5° for K19c and σ_i = 150 for K19d, following Section 9. Knowing that K19c and K19d should have a transit depth greater than 0.5 mmag (Winn 2010), a visual inspection of the Kepler light curve confirmed that the planets are not transiting.

Our solution agrees with the solution in Table 4, except for the mass of K19c, with a 2σ difference between the RV-only and RV+TTV determinations. The origin of this discrepancy is likely due to the uncertainty in the orbital phase and eccentricity from the RV-only analysis, with one planet absorbing the signal of the other, possibly coupled with the effect described in Section 10.

The RV offset between the old and new CCD is ≃2 m s⁻¹, i.e., within the RV noise. We do not expect any influence of the RV offset on the outcome of the analysis, since the three planets have periods fully contained within the time span of both old and new CCD data sets.

In Section 8 we assumed coplanarity with K19b (the only planet with known inclination) to derive the orbital parameters of the additional planets in the system. To check if this assumption was still valid after determining the orbital period of the additional planets, we determined the upper limit on inclination for which a transit is visible for a given planet by inverting Equation (7) of Winn (2010) with the assumption for the impact parameter b_tra = 1. To properly take into account variation with a time of e and a ω induced by dynamical interactions, we selected 1000 (e, ω) pairs for each planet by randomly sampling in time the integration of our best solution over the Kepler observational time span. Error in the stellar radius was included by generating random samples from a normal distribution with mean R_⋆ and standard deviation ${\sigma }_{{{\rm{R}}}_{\star }}$ as reported in Table 3. We obtained i_min = 8872 ± 003 for K19c and i_min = 8924 ± 002 for K19d. Since we assumed i_b = 89.94 (from B11),  the coplanarity assumption cannot hold for K19c. We assessed the influence of orbital inclinations on the validity of our solution by running dynamical simulations for a grid of i_c and i_d (from 60 to 120 degrees, with steps of 025 for i_c and 05 for i_d), with the remaining orbital parameters fixed to our best solution; and we determined the reduced χ² as in Section 8. As shown in Figure 8, the reduced χ² increases rapidly, with K19c going farther from coplanarity with K19b, reaching a value of ≈1.6 for $| {i}_{b}-{i}_{c}| \simeq 5^\circ$, while K19d can span a larger interval in inclinations without affecting the outcome (although increasing $| {i}_{b}-{i}_{d}|$ would negatively affect the stability of the system). Assuming i_c = 88.72 (grazing scenario) the outcome of the fit is nearly the same, ${\chi }_{\mathrm{red}}^{2}=1.26$, which is very close to the value obtained when assuming coplanarity (χ² = 1.25; see Section 8). It is likely that the system is very close to orbital alignment, as observed for the majority of transiting multi-planet systems, e.g., Figueira et al. 2012; Fabrycky et al. 2014; Ballard & Johnson 2016; Becker & Adams 2016.

Zoom In Zoom Out Reset image size

Dynamical simulations include by definition the effects of gravitational interactions between planets. This is different from the usual approach followed in the exoplanet literature, where a series of non-interacting Keplerian orbits is used to derive the planet parameters in multiple systems. While the assumption of negligible interactions between planets may hold in most of the cases, in the presence of TTVs we know that such interactions are happening. It is then worthwhile to analyze the differences in the RVs between the two approaches, i.e., dynamical versus non-interacting Keplerian orbits. In order to do so, we have simulated the expected RVs of the Kepler-19 system at the observational epochs of our data set, using the planetary parameters in Table 5 and assuming non-interacting orbits, and we have subtracted the outcome from the dynamically derived RVs from the same orbital parameters. The results are shown in Figure 9. The ${\mathrm{RV}}_{\mathrm{TTV}}-{\mathrm{RV}}_{\mathrm{Kep}}$ residuals show a peak-to-peak variation of 0.30 m s⁻¹, with a prominent periodicity at 29.3 days, i.e., very close to the orbital period of K19c, although these values depend strongly on the assumed orbital parameters. Note that the difference between dynamical and Keplerian RVs diverges the further you move from the reference time of the orbital parameters T_ref. This occurs because of small variations of the orbital parameters with time caused by planet interactions, which are implicitly taken into account in dynamical integration but are ignored in the Keplerian approach. For the Kepler-19 system, after two years of observations from T_ref, the peak-to-peak difference has already reached ≃0.6 m s⁻¹, i.e., a value that can negatively affect the mass determination of the planets. Assuming that non-interacting planets when modeling the RVs can thus mislead the determination of the orbital parameters, e.g., in the case of a data set spanning several years.

Zoom In Zoom Out Reset image size

The transiting planet Kepler-19b was previously validated by Ballard et al. (2011). A period of 9.23 days and an upper limit of 20 M_⊕ were determined from the analysis of the light curve and high-resolution spectroscopy. From the presence of TTV in 8 Kepler quarters they deduced the presence of a second, non-transiting planet with a period of ≲160 days and mass of ≲6 M_jup. In this paper we presented the first precise mass measurement for K19b (8.4 ± 1.6 M_⊕) and the characterization of two non-transiting Neptune-mass companions, K19c (P = 28.73 ± 0.01 days, M = 13.1 ± 2.7 M_⊕), and K19d (P = 62.9 ± 0.2 days, M = 20.3 ± 3.4 M_⊕), obtained by simultaneously modeling TTVs and RVs through dynamical integration. We excluded stellar activity as a possible origin for the RV signal at P ≃ 28 days, using all the activity diagnostics at our disposal, including the latest reduction of the Kepler light curve. Nevertheless, we performed standard model selection between the three-planet and two-planet models, either with the hypothesis that the outer planet causes the TTVs and different assumptions for the activity signal, or  the hypothesis that the perturber has period of ≃28 days and no outer planets. In all cases the three-planet model was favored, with a high degree of confidence. A planet in a strong MMR with K19b could still produce the TTVs while having an RV semi-amplitude below our detection sensitivity, but only under the condition that the 28-day signal is due to stellar activity, which is not supported by our data.

With a period ratio P_c/P_b = 3.09, the system is very close to a 3:1 MMR. The sinusoidal shape of the TTV is then induced by the 3:1 MMR of the inner planets, with a modulation caused by the outer planet. We performed a comparison between dynamical RVs and those calculated assuming non-interacting planets, using the orbital parameters of the Kepler-19 system, and showed that, for this specific system, the difference for a data set spanning several years is at the limit of detection with the state-of-the-art instruments used for planet search and characterization.

Our new determination of K19b mass disagrees with the most likely value of 1.6 M_⊕ (semi-amplitude of 0.5 m s⁻¹), which was obtained by B11 while attempting a fit with only 8 Keck-HIRES RVs and assuming only 1 planet in the system. As a consequence of this assumption, they derived a likely RV jitter of ≃4 m s⁻¹, which is not confirmed by our analysis. In general, RV analysis performed on a small number of measurements should always be handled with extreme care, since the results could be affected by additional non-transiting planets that have little influence on the TTV of the transiting planets but have significant RV semi-amplitude, such as K19d.

K19b falls in the region of super-Earths with rocky cores and a significant fraction of volatiles or H/He gas, in opposition to the low-density planets characterized by TTVs only (Weiss & Marcy 2014; Jontof-Hutter et al. 2016), and well separated from the group of rocky planets with radii smaller than 2 R_⊕, as can be seen in Figure 10. If we assume that the planet composition is a mixture of an H/He envelope with a solar composition atop a rocky core with Earth-like rock/iron abundances, we can estimate the internal structure of K19b. By employing theoretical models from Lopez & Fortney (2014), and assuming an age of 2 Gyr, an envelope of 0.4 ± 0.3% of the total mass is required to explain the observed mass and radius of K19b. Despite the fairly high level of irradiation, its atmosphere is only moderately vulnerable to photoevaporation, due to the relatively large mass of the planet (Lopez et al. 2012). Employing the same models, we found that the mass of the primordial envelope was approximately twice (≃1%) that of the current envelope mass.

Zoom In Zoom Out Reset image size

Although the radii of the Neptune-mass planets are unknown, we can still speculate on their possible internal composition. By using the gas accretion scaling relations from Lee & Chiang (2015), we expect for K19c and K19d to have accreted significantly larger envelopes than K19b, likely unaffected by photoevaporation due to the larger masses of these two planets. Using Equation (22) from Lee & Chiang (2015), we can get a rough estimate of the sizes of these envelopes and find that K19c should have an H/He envelope of ≃4.5% of its total mass and a radius of ≃3.2 R_⊕, while K19d should have a fraction of volatiles around 10%–20%, and a radius of 4–5 R_⊕.

Our results confirm that TTV and RV techniques can converge to planetary densities similar to the ones obtained only with the RV data set, if enough data are available from both sides. This result supports the analysis of Steffen (2016), where the discrepancy between the planetary density obtained from the TTV and RV noted by Weiss & Marcy (2014) is explained as a selection effect rather than as an intrinsic problem of one of the two techniques. Previous notable examples of agreement between RV and TTV masses are represented by the WASP-47 (Becker et al. 2015; Weiss et al. 2016) and Kepler-18 systems (Cochran et al. 2011), and the independent RV confirmation by Dai et al. (2016) of the TTV-derived masses of the planets in the K2-19 system (Barros et al. 2015). It should be noted, however, that by combining TTV and RV measurements, Nespral et al. (2016) derived a higher mass for K2-19b, consistent within 1σ with the value calculated by Barros et al. (2015) only. In the case of the KOI-94 system, the RV (Weiss et al. 2013) and TTV (Masuda et al. 2013) planetary masses agree except for planet d, where the TTV mass measurement is half that of the mass obtained by high-precision RVs. Until now, most of the targets characterized with TTVs were too faint for a precise, independent mass characterization with RVs. Future space-borne missions such as TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2014) and PLATO (PLAnetary Transits an stellar Oscillations, Rauer et al. 2014) will finally shed light on this problem by providing a large number of targets bright enough for mass measurement with both techniques independently.

We thank the anonymous referee for the prompt and careful review and for providing detailed and useful suggestions. The HARPS-N project was funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen's University Belfast, and University of Edinburgh.

The research leading to these results received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement number 313014 (ETAEARTH).

This work was performed in part under contract with the California Institute of Technology/Jet Propulsion Laboratory funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute.

A.V. is supported by the NSF Graduate Research Fellowship, grant No. DGE 1144152.

This publication was made possible by a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. This material is based upon work supported by NASA under grant No. NNX15AC90G issued through the Exoplanets Research Program.

P.F. acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contract of reference IF/01037/2013 and POPH/FSE (EC) by FEDER funding through the program "Programa Operacional de Factores de Competitividade—COMPETE", and further support in the form of an exploratory project of reference IF/01037/2013CP1191/CT0001.

X.D. is grateful to the Society in Science–Branco Weiss Fellowship for its financial support.