THE MASS OF Kepler-93b AND THE COMPOSITION OF TERRESTRIAL PLANETS^*

Small planets are abundant in the galaxy, but the compositional diversity of small planets is not well understood. Theoretical models of planet formation predict that planets intermediate in size between Earth and Neptune could be gaseous "mini-Neptunes," water worlds, or rocky "Super-Earths" (Kuchner 2003; Léger et al. 2004; Valencia et al. 2006; Seager et al. 2007; Fortney et al. 2007; Rogers et al. 2011; Lopez et al. 2012; Zeng & Sasselov 2013). Recent studies have explored the compositional diversity of small planets using hierarchical Bayesian modeling of the observed planet radii and measured planet masses (Rogers 2014) or theoretical models (Wolfgang & Lopez 2014), but a thorough investigation of planet densities is hindered by the small number of small planets with well-measured masses and radii. There are currently only nine planets smaller than 2.7 R_⊕ with masses measured to 20% precision: 55 Cnc e (Gillon et al. 2012; Nelson et al. 2014), CoRoT-7b (Barros et al. 2014; Haywood et al. 2014), GJ1214b (Charbonneau et al. 2009), HD 97658b (Dragomir et al. 2013), HIP116454b (Vanderburg et al. 2014), Kepler-36b (Carter et al. 2012), Kepler-78b (Pepe et al. 2013; Howard et al. 2013), and Kepler-10b and 10c (Dumusque et al. 2014).

The host star Kepler-93 (KIC 3544595, KOI 69) is one of the brightest stars observed by Kepler (V = 10.2, Kp = 9.93), enabling very high precision photometry of 17 ppm on six-hour timescales (Christiansen et al. 2012). Kepler observed Kepler-93 throughout the baseline mission (Quarters 0–17) and conducted observations at short cadence (exposure time of 58.5 s) beginning in Quarter 2 and extending until the end of the mission. Due to the high photometric precision of the Kepler-93 observations, the planet was detected in the first four months of Kepler data (Borucki et al. 2011). Marcy et al. (2014) acquired 32 Keck HIRES radial velocity (RV) observations of Kepler-93 from 2009 July–2012 September and provided an estimate of 2.6 ± 2.0 M_⊕ for the mass of Kepler-93b. Marcy et al. (2014) also noted a large linear RV trend of 11.2 ± 1.5 m s⁻¹ yr⁻¹ and calculated lower limits on the mass and period of the perturbing companion of M > 3 M_Jup and P > 5 yr. Incorporating an additional 14 spectra from the 2013 observing season, the HIRES mass estimate for Kepler-93b increased to 3.8 ± 1.5 M_⊕ (Ballard et al. 2014). Nonetheless, the 40% error on the mass measurement allows a wide range of planetary compositions including a rocky body, an ice world, and even a substantial primordial envelope of hydrogen and helium (Ballard et al. 2014).

In contrast, the properties of the host star Kepler-93 are well-constrained. Using 37 months of Kepler short cadence data, Ballard et al. (2014) conducted an asteroseismic investigation to characterize Kepler-93 in exquisite detail. They estimated an average stellar density of 1.652 ± 0.006 g cm⁻³, a stellar mass of 0.911 ± 0.033 M_☉, and a stellar radius of 0.919 ± 0.011 R_☉. Adopting priors from their asteroseismic investigation, they fit the Kepler photometry to obtain a precise radius estimate of 1.478 ± 0.019 R_⊕ for Kepler-93b.

In addition to characterizing the host star, Ballard et al. (2014) present a variety of evidence that Kepler-93b is a bona fide planet rather than an astrophysical false positive. First, they report that the steep shape of the ingress and egress portions of the Kepler-93b light curve cannot be reproduced by a non-planetary companion. Second, they note that the infrared transit depth they measured with the Spitzer Space Telescope is consistent with the planetary interpretation of Kepler-93b. Third, they place stringent limits on the presence of nearby stars based on Keck AO images (Marcy et al. 2014). Fourth, they state that the stellar density derived from the transit duration (Seager & Mallén-Ornelas 2003; Nutzman et al. 2011) is consistent with the asteroseismic stellar density constraint, indicating that the planet likely orbits the target star rather than the companion causing the large RV trend.

In this paper, we refine the mass measurement of Kepler-93b from 2.5σ to 6σ by analyzing two seasons of HARPS-N radial velocities in addition to the publicly available HIRES data. We discuss these observations and our data reduction methods in Section 2. In Section 3, we develop a model to fit the observed radial velocities. Finally, we discuss the implications of the resulting planet mass and present our conclusions in Section 4.

We obtained 86 spectra of Kepler-93 using the HARPS-N spectrograph on the 3.57 m Telescopio Nazionale Galileo (TNG) at the Observatorio del Roque de los Muchachos. HARPS-N is a high-precision, vacuum-stabilized, high-resolution (R ≃ 115, 000) echelle spectrograph. The design is very similar to the design of the original HARPS instrument at the ESO 3.6 m (Mayor et al. 2003). The main differences are that HARPS-N is fed by octagonal fibers rather than circular fibers to improve the scrambling of the light and features a monolithic 4096 × 4096 CCD instead of the dual CCD configuration used for the HARPS focal plane (Cosentino et al. 2012).

We acquired 38 and 49 HARPS-N observations of Kepler-93 during the 2013 and 2014 observing seasons, respectively. In most cases, we used an exposure time of 30 minutes and achieved a mean S/N per extracted pixel of 103 at 550 nm. (Four of the spectra had an exposure time of 15 minutes and one had 27 minutes; these were all gathered in 2013 July.) One of the observations collected in 2013 was contaminated by light from a mercury lamp and was therefore removed from the analysis. The final HARPS-N data set analyzed in this paper consists of 86 spectra. In most cases (75 of 86 spectra), we observed Kepler-93 using simultaneous thorium argon (observing mode HARPN_ech_obs_thosimult). The remaining 11 observations were obtained without simultaneous thorium argon in observing mode HARPN_ech_obs_objAB.

We reduced the data with the standard HARPS-N pipeline by cross-correlating the observed spectra with a numerical mask based on the spectrum of a G2V star (Baranne et al. 1996; Pepe et al. 2002). We provide the resulting RVs and their 1σ errors in Table 1 along with the observation BJDs, exposure times, bisector spans, and stellar activity levels as measured by the Ca ii $\log (R^{\prime }_{\rm HK}$) activity indicator (Noyes et al. 1984). The BJDs in Table 1 are provided in UTC, but we converted the times to TDB (the units used by the Kepler mission) using the IDL routine utc2bjd.pro¹⁷ prior to fitting the RVs. We did not find evidence for a correlation between RV and bisector span or $\log (R^{\prime }_{\rm HK})$.

Table 1. HARPS-N Radial Velocity Observations of Kepler-93

BJDUTC	RV (m s−1)		Bisector	log($R^{\prime }_{{\rm HK}}$) (dex)		texp
−2,450,000	Value	Error	(m s−1)	Value	Error	(s)
2456462.686262	27335.24	1.02	−28.39	−5.01	0.01	1800
2456463.584483	27337.75	0.94	−31.63	−5.00	0.01	1800
2456464.609617	27331.57	1.75	−27.24	−5.04	0.02	1800
2456465.606438	27342.34	1.00	−32.74	−5.02	0.01	1800
2456466.608850	27337.58	0.86	−29.26	−5.00	0.01	1800
⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Our full data set included RVs from four seasons of HIRES observations (2009 July–2012 September) and two seasons of HARPS-N observations (2013 June–2014 October). We fit the combined HARPS-N and HIRES data set by incorporating a single offset RV_off between the HIRES and HARPS-N data. We used the following general model:

$$\begin{eqnarray} \mathcal {M}(t_i)& =& \gamma + {\rm RV}_{\rm off} + \beta (t_i) \nonumber\\ && +K \left[\cos (\theta (t_i, T_C, P, e) + \omega) + e \cos {\omega }\right], \end{eqnarray} \tag{ 1 }$$

where γ is the systemic velocity of Kepler-93, RV_off = RV_HARPS-N − RV_HIRES is the offset between the HIRES and HARPS-N RVs, β(t_i) is a long-term RV trend due to a third component in the system, K is the semi-amplitude due to Kepler-93b, and ω is the argument of periastron. The function θ is the true anomaly of Kepler-93b at time t_i and depends on the period P, epoch of transit T_C, and eccentricity e. When fitting eccentric orbits, we used IDL routine keplereq.pro written by Brian Jackson¹⁸ to solve Kepler's equation for the eccentric anomaly. The routine uses the method suggested by Mikkola (1987) as an initial guess.

We considered linear and quadratic parameterizations of the long-term trend β(t_i) and circular and eccentric orbits for Kepler-93b. For all models, we determined an initial solution using the Levenberg–Marquardt minimization algorithm as implemented by lmfit in IDL. We then explored the region of parameter-space near the best-fit solution using a Bayesian Markov Chain Monte Carlo (MCMC) analysis with a Metropolis–Hastings acceptance criterion (Metropolis et al. 1953). We initialized N chains, where N was twice the number of free parameters in the chosen model. We selected different initial positions for each chain by perturbing each free parameter of the best-fit solution by a random number drawn from a distribution with a width of five times the step size. We tuned the step sizes such that the acceptance fractions for each parameter were 10%–30%. For the MCMC analysis, we set uniform priors for all parameters except the orbital period and epoch of transit. We allowed only non-negative values for the RV semi-amplitude K and the separate stellar jitter terms σ_sj for the HIRES and HARPS-N observations (see below).

The Kepler photometry places tight constraints on the period and epoch of transit (Ballard et al. 2014). We incorporated this knowledge into our MCMC analysis by including Gaussian priors on the period and transit epoch in the likelihood calculation. As shown in Dumusque et al. (2014), using the tight prior from Kepler photometry when fitting a circular model to RV observations of a planet in a circular orbit yields a result very similar to that from a combined photometric and spectroscopic fit. We also tested fitting the data while allowing the epoch of transit to float and find the epoch of transit at BJD = 2,454,944.29514. This epoch differs from the value determined by Ballard et al. (2014) by 4 minutes (0.3σ). The possible shift in the transit center is therefore insignificant. Accordingly, we adopt the photometric ephemeris determined by Ballard et al. (2014).

In our calculations, we shifted the epoch of transit close to the start of the HARPS-N RVs to reduce error propagation. We increased the efficiency of our model fits by parameterizing eccentric models using $\sqrt{e} \cos (\omega)$ and $\sqrt{e} \sin (\omega)$ rather than varying e and ω directly (Ford 2006; Eastman et al. 2013). As in Dumusque et al. (2014), we accounted for stellar activity by incorporating a stellar jitter term σ_sj in our adopted likelihood $\mathcal {L}$:

$$\begin{equation} \mathcal {L}= \prod _{i=1}^{N}\left({\frac{1}{\sqrt{2 \pi \left(\sigma _i^2 + \sigma _{sj}^2\right)}} \exp \left[ - \frac{({\rm RV}(t_i) - \mathcal {M}(t_i))^2}{2 \left(\sigma _i^2+\sigma _{sj}^2\right)}\right] }\right), \end{equation} \tag{ 2 }$$

where RV(t_i) is the measured RV at each time t_i in the set of N observations, $\mathcal {M}$ is the model, σ_i is the instrumental noise listed in Table 1, and the stellar jitter noise σ_sj is allowed to adopt a different constant value for the HARPS-N and HIRES data.

We ran each chain for a minimum of 10⁴ steps and checked for convergence by computing the Gelman–Rubin potential scale reduction factor $\hat{R}$ for each parameter (Gelman et al. 2004). We stopped the MCMC analysis when $\hat{R}< 1.03$ for all parameters. Next, we accounted for "burn-in" by identifying the point in each chain at which the likelihood first became higher than the median likelihood of the chain and removing all earlier steps. After combining all of the chains, we selected the median values of each parameter as the best-fit value and assigned symmetric errors encompassing 68% of values closest to the adopted best-fit value.

We then used Bayesian statistics to determine which of the models considered best describes the data. We followed the method of Chib & Jeliazkov (2001) as described in the Appendix of Haywood et al. (2014) to calculate the Bayes factor between pairs of models using the posterior distributions and acceptance probabilities from our MCMC analyses. This method was previously used by Dumusque et al. (2014) to compare RV models of the Kepler-10 system. We found that penalties incurred by the additional complexity of fitting the orbit of Kepler-93b with an eccentric model or fitting the long-term trend with a quadratic model outweighed the improvement in the likelihood. We also compared the models by holding the stellar jitter terms fixed to σ_{sj, HARPS-N} = 1.56 m s⁻¹ and σ_{sj, HIRES} = 2.03 m s⁻¹ and computing the Bayesian Information Criterion (BIC; Schwarz 1978) and finite sample Akaike Information Criterion (AIC_C; Hurvich & Tsai 1989). When considering only the HARPS-N data, we found that the model with a quadratic trend and a circular orbit for Kepler-93b was preferred over the models with a linear trend and circular orbit (ΔBIC = 6.1, ΔAIC_C = 8.3), linear trend and eccentric orbit (ΔBIC = 12.4, ΔAIC_C = 10.4), or quadratic trend and eccentric orbit (ΔBIC = 6.6, ΔAIC_C = 2.5).

Nonetheless, when we included the HIRES data we found that the simplest model (linear trend and circular orbit) was preferred over the models with a linear trend and eccentric orbit (ΔBIC = 5.6, ΔAIC_C = 0.7), quadratic trend and circular orbit (ΔBIC = 4.7, ΔAIC_C = 2.1), or quadratic trend and eccentric orbit (ΔBIC = 10.5, ΔAIC_C = 3.0). We therefore treat the perturbation from Kepler-93c as a linear trend and model the orbit of Kepler-93b as circular. (For the eccentric fits, we found a median eccentricity of 0.15 and an upper limit of e < 0.31 with 95% confidence.) We present the resulting system properties including a mass estimate for Kepler-93b of 4.02 ± 0.68 M_⊕ in Table 2 and display the measured RVs and the best-fit model in Figure 1. As highlighted in Figure 2, the HARPS-N residuals are Gaussian with a distribution centered on zero and containing 68% of the data within a half width 1.6 m s⁻¹. For the HIRES residuals the region encompassing 68% of the data has a half width of 3.4 m s⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Parameters for the Kepler-93 System

Parameter	Value and	Ref.
	1σ Errors
Kepler-93 (star) = KIC 3544595 = KOI 69
Right ascension	$19^{\rm h}25^{\rm m}40\buildrel{\mathrm{s}}\over{.}39$	1, 2
Declination	$+38^{\rm d}40^{\rm m}20\buildrel{\mathrm{s}}\over{.}45$	1, 2
Kepler magnitude	9.931	3
2MASS K	8.370	3
Teff (K)	5669 ± 75	4
R* (solar radii)	0.919 ± 0.011	4
M* (solar masses)	0.911 ± 0.033	4
[Fe/H]	−0.18 ± 0.10	4
log g	4.470 ± 0.004	4
Age (Gyr)	6.6 ± 0.9	4
Systemic Velocitya (m s−1)	27337.89 ± 0.51	5
HIRES Offset (m s−1)	27304.1 ± 1.5	5
RV Jitter (HARPS-N)	1.58 ± 0.19	5
RV Jitter (HIRES)	2.09 ± 0.71	5
Kepler-93b (planet) = KOI 69.01
Transit and orbital parameters
Orbital period P (days)	4.72673978 ± 9.7 × 10−7	4
Transit epoch TC (BJD)	2454944.29227 ± 0.00013	4
Rp/R*	0.014751 ± 0.000059	4
a/R*	12.496 ± 0.015	4
Inc (deg)	89.183 ± 0.044	4
Impact parameter	0.1765 ± 0.0095	4
Orbital eccentricity e	0 (fixed)	5
RV semi-amplitude K (m s−1)	1.63 ± 0.27	5
Planetary parameters
Rp( R⊕)	1.478 ± 0.019	4
Mp( M⊕)	4.02 ± 0.68	5
ρp(g cm−3)	6.88 ± 1.18	5
log gp (cgs)	3.26 ± 0.07	5
a (AU)	0.053 ± 0.002	5
Teq (K)b	1037 ± 13	4
Kepler-93c (companion)
Fit parameters
Acceleration (m s−1 yr−1)	12.0 ± 0.4	5
Companion limits
Mass (MJ)	>8.5	5
Orbital period P (yr)	>10	5

Notes. ^aSystemic velocity at BJD 2,456,461.57573945. ^bAssuming a Bond albeo of 0.3. References. (1) Høg et al. 1998; (2) Høg et al. 2000; (3) Brown et al. 2011; (4) Ballard et al. 2014; (5) this paper.

Download table as:  ASCIITypeset image

The expected circularization timescale for Kepler-93b is significantly shorter than the 6.6 ± 0.9 Gyr age of the star (Ballard et al. 2014). Following Goldreich & Soter (1966), we calculated a tidal circularization timescale of 75 Myr for a 4.02 M_⊕, 1.48 R_⊕ planet in an orbit with a = 0.053 AU around a 0.91 M_☉ star. We assumed Q = 100 based on the tidal quality factors estimated for terrestrial planets in the Solar System (Yoder 1995; Henning et al. 2009). Obtaining a tidal circularization timescale similar to the age of the system would require Q = 9000, comparable to the estimate for Neptune (Zhang & Hamilton 2008). Although the tidal circularization argument is consistent with the preference for a circular orbit, we caution that the tidal quality factors for exoplanets are largely unknown.

3.1. Limits on the Properties of Kepler-93c

The baseline of our RV data is too short to measure the period and minimum mass of the perturber responsible for the long-term trend, but we can place lower limits on the companion properties. Wang et al. (2014) conducted a similar analysis of the properties of Kepler-93c based on AO observations and Keck/HIRES RVs. They found a linear RV trend of 12.2 ± 0.2 m s⁻¹ yr⁻¹ and argued that Kepler-93c is most likely to have a mass below 101 M_J and a semi major axis a = 15.5–33 AU if it is a stellar companion. For the substellar case, they found limits of a = 5.5–27.6 AU and M = 10–80 M_J.

Our additional two years of HARPS-N observations have allowed us to further restrict the allowed parameter space for Kepler-93c. We measured a linear trend of 12.0  ±  0.4 m s⁻¹ yr⁻¹ for 5 years, implying that Kepler-93c has P  >  10 yr and M  >  8.5 M_J. Assuming the 100 pc distance to Kepler-93 estimated by Ballard et al. (2014), the resulting semimajor axis a > 4.5 AU corresponds to an angular separation of 0045. At this separation, the detection limit from Keck AO imaging is 1.7 Ks magnitudes fainter than Kepler-93. We can therefore place an upper limit of Ks > 10.1 on Kepler-93c unless Kepler-93c happened to have an orbital geometry precluding detection at the epoch of the Keck observations. Converting the Ks upper limit into a mass limit via the Delfosse et al. (2000) relation¹⁹ and the distance, we found a mass upper limit of 0.64 M_☉ for angular separations beyond 0045. We display the combined limits from the AO and RV data in Figure 3. In the future, astrometric measurements from Gaia (Perryman et al. 2001) will likely provide additional constraints on the properties of the Kepler-93 system. We will then be able to investigate the dynamical history of the system and test whether Kepler-93c might be responsible for scattering Kepler-93b inward onto a short-period orbit.

Zoom In Zoom Out Reset image size

Combining our estimate of 4.02 ± 0.68 M_⊕ for the mass of Kepler-93 with the radius estimate of Rp = 1.478 ± 0.019 R_⊕ from Ballard et al. (2014), we find a density of 6.88 ± 1.18 g cm⁻³. In Figure 4, we show Kepler-93b on the mass-radius diagram. In this diagram we plot only those planets smaller than 2.7 R_⊕ and with masses determined to a precision better than 20%. In addition to Venus and the Earth, there are 10 such planets. We observe that Kepler-93b falls in a cluster of planets with radii 50% larger than that of the Earth, all of which have extremely similar densities: Kepler-10b (ρ = 5.8 ± 0.8 g cm⁻³; Dumusque et al. 2014), Kepler-36b ($\rho = 7.46^{+0.74}_{-0.59}$ g cm⁻³; Carter et al. 2012), and CoRoT-7b (ρ = 6.56 ± 1.40 g cm⁻³; Barros et al. 2014; Haywood et al. 2014). This cluster falls upon a relation that includes Earth, Venus, and Kepler-78b (Howard et al. 2013; Pepe et al. 2013), which is itself only 20% larger than the Earth. To investigate this further, we used the two-component iron-magnesium silicate models of Zeng & Sasselov (2013) to see if we could find a single composition that explained these seven worlds.

Zoom In Zoom Out Reset image size

For the solar system planets, we artificially include mass and radius errors equal to the mean fractional errors for the exoplanets considered so that they do not have undue influence on the resulting fit. We find the lowest χ² for a model composition of 83% MgSiO₃ and 17% Fe. We arrive at the same best-fit relation when we exclude Earth and Venus. We caution that the two-component models used in this analysis make two simplifying approximations about the interior structure of planets that cause the core mass fraction to be underestimated: (1) the core contains only iron and the mantle contains only magnesium silicate and (2) the planet is completely dry with no water content. Accordingly, we expect the actual core mass fraction to be slightly higher by 5%–8% to account for incorporation of lighter elements like oxygen, sulfur, and silicon in the core and the inclusion of water in the mantle. In addition, there could be a change of roughly 2% toward higher or lower core fractions due to uncertainties in the equations of state used in the model calculations. Our purpose in this exercise is to test whether we can find one composition that successfully explains all seven planets, not to place stringent constraints on the abundance of magnesium silicate or iron.

Intriguingly, all of these planets, which are smaller than 1.6 R_⊕, have a tight dispersion around this best-fit compositional curve, suggesting that the distribution of small planet compositions has low intrinsic scatter. In the solar system, the strong agreement between abundance ratios of elements in meteorites and those of the solar photosphere (Lodders 2003) is a key constraint by which we deduce the composition of the interior of the Earth. Therefore, we might look to the bulk abundances of exoplanet host stars for similar constraints on the interior compositions of their terrestrial planets. Grasset et al. (2009) use a set of planetary models to investigate the dependence of planet radii on elemental abundances. Varying the ratios of iron to silicate and magnesium to silicate within the range observed for the photospheric abundances of nearby exoplanet host stars (Beirão et al. 2005; Gilli et al. 2006), Grasset et al. (2009) predicted that the radii of terrestrial planets would vary by roughly 2% at a given mass. Our findings are in agreement with this picture: We measure a mean absolute deviation of 1.9% between the estimated planet radii and the values predicted by a 83% MgSiO₃/17% Fe model for planets less massive than 6 M_⊕. Indeed, rocky planets very close to their host stars seem to obey a well-defined relationship between radius and mass, although with only five such examples outside the solar system, the immediate task is to characterize other terrestrial exoplanets with similar precision. Increasing the sample of small planets with well-constrained masses and radii will allow us to learn whether additional rocky planets could also be explained by a single mass-radius relation and investigate whether the relation found for close-in planets extends to planets in more distant orbits.

Our mass-radius diagram also includes five planets more massive than 6 M_⊕: 55 Cnc e, GJ1214b, HD 97658b, HIP 116454b, and Kepler-10c. In contrast, none of these more massive planets have a high density consistent with the best-fit magnesium silicate/iron composition described above. In agreement with Rogers (2014), we find that planets larger than approximately 1.6 R_⊕ (e.g., more massive than approximately 6 M_⊕) contain significant fractions of volatiles or H/He gas. These planets appear to have a diversity of compositions that is not well-explained by a single mass-radius relation (Wolfgang & Lopez 2014).

The discussion above focused exclusively on planets smaller than 2.7 R_⊕ with masses measured to better than 20%. Some low-mass worlds with very low densities are known, notably the Kepler-11 system (Lissauer et al. 2013) and KOI-314c (Kipping et al. 2014). Thus we are not proposing that all planets less massive than 6 M_⊕ obey a single mass-radius relation; rather, we suggest that the rocky analogs of the Earth might do so.

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), and the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen's University Belfast, and University of Edinburgh. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 313014 (ETAEARTH). C.D. is supported by a National Science Foundation Graduate Research Fellowship. X.D. would like to thank the Swiss National Science Foundation (SNSF) for its support through an Early Postdoc Mobility fellowship. P.F. acknowledges support by Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contracts of reference IF/01037/2013 and POPH/FSE (EC) by FEDER funding through the program "Programa Operacional de Factores de Competitividade - COMPETE." This publication was made possible through the support of 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.