© ESO 2018

1 Introduction

Kepler-419 (KOI-1474) is a planetary system with two known giant planets. The inner one, Kepler-419b, was first discovered transiting in the Kepler photometry by Borucki et al. (2011) with a period of 69.7 days and an estimate radius of 1.0 R_J. Its size and relatively long period places it in the “Period Valley” of giant planets (Udry et al. 2003; Batygin et al. 2016). Strong transit-timing variations (TTVs) of the order of an hour were detected later by Ford et al. (2012) and Dawson et al. (2012). Dawson et al. (2012) concluded that the observed TTVs were consistent with perturbations from a massive, eccentric outer companion in the system, but they could not constrain the outer body’s orbital period or mass because of the small number of transits and poor orbit coverage Kepler data had at the time. Additionally, they validated the planetary nature of the transits, determined through the photoeccentric effect, that the orbit of the transiting planet is highly eccentric (), and they found that the host star is a rapidly rotating F7 star with a rotational period of P_rot = 4.6 ± 0.4 days. Some years later, using 11 quarters of Kepler photometry, Mazeh et al. (2015) measured a rotational period of P_rot = 4.53 ± 0.16 days.

In 2012 we started a radial velocity follow-up of Kepler-419 with the SOPHIE spectrograph in order to detect and characterise the outer companion and to refine the parameters of the inner, transiting planet. This target is part of the sample presented in Santerne et al. (2016). During our follow up campaign, Dawson et al. (2014) performed aTTV analysis of 16 quarters of Kepler data (containing all transits observed by Kepler), which allowed them to ascertain that the perturber object is a non-transiting planet, namely Kepler-419c, with a mass of 7.3 ± 0.4 M_J on an eccentric orbit (e = 0.184 ± 0.002) with a period of 675.47 ± 0.11 days. They also presented 20 radial velocities secured with the HIRES spectrograph, allowing the mass of Kepler-419b to be measured at 2.5 ±0.3 M_J and confirming the photometrically determined eccentricity. The HIRES data also showed an additional acceleration consistent with the Kepler-419c properties derived from TTVs, but were not able to detect the exterior planet independently.

The system is remarkable because it harbours an almost coplanar pair of giants for which the innermost planet has an extremely high eccentricity. Since mechanisms usually invoked to explain highly eccentric orbits require the presence of a companion on a significantly inclined orbit, Kepler-419 presents a challenge to theory to explain its origin. Moreover, the apsidal lines are close to anti-aligned, the ascending node longitudes are close to aligned, and integrations confirm that this configuration persists over secular timescales. The state of the system is highly suggestive of gentle relaxation at some point in its history, with the source of dissipation coming from either the protoplanetary disk or tides in planet b (or both). However, the evolutionary path to its current state is still not clear.

Here wepresent our SOPHIE radial velocity measurements, comprising 45 epochs over 2.2 years. Both planets are detected independently in the new data set, confirming the TTV detection of the exterior planet. In order to explore the contribution of each data set, we analyse the SOPHIE radial velocities with a simple two-Keplerian model and the Kepler photometry independently of the SOPHIE data.

On the other hand, and as described in detail by Agol et al. (2005) and Almenara et al. (2015, 2016), photometry alone makes it possible to measure the density of the bodies of the system, that is MR⁻³. However, individual masses and radii cannot be constrained, and we only have access to mass ratios and radius ratios. In other words, exactly the same light curve is obtained if lengths in the system are scaled by a factor and masses are scaled by the same factor at cubic exponent. This is called the MR⁻³ degeneracy. This degeneracy can be removed by constraining the system scale, for example by adding radial velocities, or by measuring the light-travel time, which provides access to absolute masses and radii. Then we combine all available data and use a photodynamical model (Carter et al. 2011) to derive absolute physical parameters without theoretical stellar models (Agol et al. 2005). Finally, the results from the photodynamical modelling are used to study the evolution of the system’s orbital parameters through numerical integrations over 10 kyr.

2 Data

2.1 Kepler light curve

Kepler observed the 13-magnitude star Kepler-419 from Q0 to Q17. We used Data Release 25 obtained from the Mikulski Archive for Space Telescopes (MAST) archive¹. Kepler recorded 21 of the 22 transits of planet b that took place in the time span of its observation, the first ten and the last one in long-cadence data (about one point every 29.4 min) and the remaining ten in short-cadence data (about one point per minute). We used the simple aperture photometry (SAP) light curve, which we corrected for the flux contamination (between 0 and 2% depending on the quarter) using the value estimated by the Kepler team. We kept only the data spanning three transit durations around each transit; they were modelled after normalisation using a linear function for each transit. The transit observations are presented in Fig. 1.

2.2 Radial velocities

We observed the star Kepler-419 with SOPHIE (Perruchot et al. 2008; Bouchy et al. 2009) at the 1.93 m telescope of the Observatoire de Haute-Provence (France). Observations were secured in the slow readout mode of the detector and in high-efficiency mode of the spectrograph, with a resolution power of λ∕Δ λ = 39 000. The first optical fibre was used for starlight, whereas the second fibre was placed on the sky to evaluate the sky background pollution, especially from moonlight. Wavelength calibrations were secured approximately every two hours during the night to monitor and correct for the potential spectrograph drifts. Forty-five exposures of Kepler-419 were secured between May 2012 and July 2014. Most of them have 1-hour exposure time and their signal-to-noise ratios (S/Ns) per pixel at 550 nm range from 16 to 43, with a typical value of 30. This translates into a mean radial velocity precision of 28 m s⁻¹, estimated following Boisse et al. (2010).

The spectra were extracted using the SOPHIE pipeline (Bouchy et al. 2009), and cross-correlated with a G2-type numerical mask to produce cross-correlation functions (CCFs). CCFs are fit with Gaussians to derive the radial velocities (Baranne et al. 1996; Pepe et al. 2002). The cross-correlation with other types of numerical masks (e.g. K0 ou K5) does not significantly change the results. A series of corrections were subsequently applied, following Santerne et al. (2016): CCD charge transfer inefficiency was corrected (as exposures have different S/Ns), background pollution was removed (which produced corrections of up to 40 m s⁻¹ on ten affected exposures), and instrumental drifts were subtracted using the monitoring of the constant star HD 185144 (which shows a dispersion of 6.5 m s⁻¹ over our 2.2-year time span). The resulting radial velocities are listed in Table A.1 and plotted in Fig. 2.

The observed CCFs are broad; they have a full width at half maximum of 19.9 ± 0.2 km s⁻¹, which corresponds to a projected rotational velocity v sin i_⋆ = 12.3 ± 1.0 km s⁻¹ (Boisse et al. 2010). From two high-resolution HIRES spectra Dawson et al. (2014) measured v sin i_⋆ = 14.4 ± 1.3 km s⁻¹, in agreement with the SOPHIE estimate. Using these estimates and the amplitude of the photometric variability, the expected stellar jitter produced bystellar spots on this rapidly rotating F7 star can be estimated using a simplified formula for the amplitude of the Rossiter–McLaughlin effect². The Kepler light curve presents a variability of ~2 mmag peak-to-peak amplitude, which leads to an estimate radial velocity jitter of ~ 27 m s⁻¹. Additionally, the time series of the bisector velocity span shows a clear periodicity at 4.58 days (Fig. 3), in agreement with the measured rotational period from Kepler photometry (Mazeh et al. 2015).

For our analyses we also used the 20 HIRES (Vogt et al. 1994) radial velocities presented by Dawson et al. (2014), acquired between April and September 2012, with one additional point on August 2013. Their typical internal precision is of the order of ± 12 m s⁻¹. However, Dawson et al. (2014) report a scatter of 40 m s⁻¹ around their best-fit model, which they attribute to stellar effects. By comparison with the HIRES data set, the new SOPHIE data provide a longer time span with an improved time sampling.

Fig. 1 Transits of Kepler-419b observed by Kepler. Each panel is centred at the linear ephemeris (indicated by the vertical grey lines, and reported in the caption of Fig. 8). For short-cadence data, 29.4-minute binned data is shown in addition to the observed data points. Each panel is labelled with the epoch; zero is the first transit after tref. The black curve is the median oversampled model over 10 000 random MCMC steps. In the lower part of each panel the residuals after subtracting the MAP model to the observed data are shown. The shades of grey represent the 68.3, 95.5, and 99.7% credible intervals, which are hardly distinguishable in the residual panels and show an increased uncertainty at ingress and egress times.

Fig. 2 Radial velocities of Kepler-419 observed by HIRES (upper plot, purple points) and SOPHIE (lower plot, blue points). Theblack curve is the median model from the photodynamical analysis, with the 68.3, 95.5, and 99.7% credible curves.

Fig. 3 Generalised Lomb–Scargle periodogram of the bisector velocity span of the SOPHIE observations. A clear peak is seen at P = 4.58 days, in agreement with the rotational period determination by Mazeh et al. (2015), whose 1 − σ range is shown in grey.

3 Analysis

3.1 First analysis: Keplerian model

We first fitted the radial velocities with a two-planet Keplerian model, i.e. neglecting the mutual gravitational interactions between both planets. The goal here was to have a first idea of the information included in the radial velocity data set.

Figure 4 shows Generalised Lomb–Scargle (GLS) periodograms (Zechmeister & Kürster 2009) of the radial velocities. In each panel the periods of the two planets are indicated with the green and the red lines (70 and 673 days, respectively). The first panel shows the periodogram of the mean-corrected SOPHIE and HIRES velocities together. It reveals no significant signals, but there is a peak close to the rotational period of the star, 4.5 days. The periodogram of the HIRES data is presented in second panel of Fig. 4. No significant power is seen at the periods of any of the planets, and there is a forest of peaks around the rotational period, including a seemingly significant one around 7 days. We cannot explain this signal easily, but as it is only detected in the HIRES data, we are inclined to ascribe it to instrument systematics (see also discussion at the end of the section). The third panel of Fig. 4 presents the periodogram of SOPHIE data only. The p-values of the peaks atthe planet periods, computed by randomly shuffling the model residuals, remain above the customary 1%. However, when the Keplerian orbit is eccentric, the GLS periodogram power is transferred partially to the harmonics of the orbital period, leading to a reduced peak amplitude. More robust methods should be used to assess the question of the significance of these signals in detail, but this is beyond the scope of this paper.

The residuals of a Keplerian fit including Kepler-419b only show a strong signal at the period of Kepler-419c, with a p-value smaller than 1% (Fig. 4, fourth panel). We conclude that the outer, non-transiting planet Kepler-419c, which was only predicted from TTVs, is detected in the SOPHIE radial velocities. In the sameway, if the effect of the outer planet is removed from the SOPHIE data, a significant peak (p-value slightly higher than 1%) appears at the period of the inner transiting planet.

Table 1 shows the results of the two-planet Keplerian model fit (using DACE³; Delisle et al. 2016)to the SOPHIE data alone. The parameters of the Keplerian orbits are both in agreement with the parameters reported by Dawson et al. (2014) and with our photodynamical model presented in Sect. 3.3; in particular, the orbital periods and the phases agree.

The residuals of the two-Keplerian fit show no additional signals (Fig. 4, lower panel). The scatter of the residuals is 38 m s⁻¹, in good agreement with the precision of the observations and the expected stellar jitter⁴. This is similar to the dispersion measured in the HIRES residuals after a similar fit. As the HIRES data have a smaller internal dispersion (12 m s⁻¹), this hints atsome level of systematics in the HIRES data.

Fig. 4 Generalised Lomb–Scargle periodograms of the radial velocities. The green and red vertical lines mark the period of planet b and c, respectively. The grey horizontal lines represent the 50, 10, and 1% p-value levels. From top to bottom: HIRES and SOPHIE data; HIRES velocities; SOPHIE velocities; SOPHIE after substraction of Keplerian curve at the period of planet b; SOPHIE velocities after substraction of a two-Keplerian model with the periods of planets b and c. The periodogram plots are provided by DACE (dace.unige.ch).

3.2 Photodynamical modelling without radial velocities

The photodynamical analysis of the Kepler light curve, including the interaction between the planets, allows the planet-to-star mass ratios to be constrained. The details on the modelling are given in Sect. 3.3. Here we neglected effects related to the time of travel of light, which we confirmed does not change the result significantly (see Sect. 3.3). As a consequence, the model used in this section does not depend on the sizes and masses of the star and planets, only on their densities. In Fig. 5, we present the posterior distributions of the mass ratios of planets b and c relative to the star. Contrary to the commonly accepted notion that TTVs of a given planetare completely insensitive to its own mass, the distribution in Fig. 5 (left panel, black histogram) shows that an upper limit for the mass ratio of planet b can be set (0.010 at 95% confidence level). We study the constraints on the mass of planet b in detail in Sect. 5.1. On the other hand, the planet-to-star mass ratio of planet c is well constrained using Kepler photometry data alone (Fig. 5, right panel).

Fig. 5 Posteriors of planet-to-star mass ratios for planet b (left) and planet c (right). The black histograms correspond to the analysis of the photometric data alone, while the grey histogram includes radial velocities. The filled grey histograms show the PDF for the masses of planets b and c, respectively, for the analysis of the simulated light curve with 54 transits and no radial velocities.

3.3 Finalanalysis: photodynamical model

Finally, we employed a photodynamical model of the observed photometry and radial velocity measurements, accounting for the gravitational interactions of all the known components of the system. Our model is described in detail in Almenara et al. (2015, 2016). In brief, we obtain the positions and velocities of the three bodies in the system in time through numerical integration of the system. The sky-projected positions are used to compute the light curve (Mandel & Agol 2002) using a quadratic limb-darkening law (Manduca et al. 1977). To account for the integration time, the model was oversampled by a factor of 30 and 3, for the long- and short-cadence data, respectively, and then binned back to match the cadence of the data points (Kipping 2010). The line-of-sight projected velocity of the star issued from the integration is used to model the radial velocity measurements (i.e. we do not assume Keplerian motion).

We used the n-body code rebound (Rein & Liu 2012) with the whfast integrator (Rein & Tamayo 2015) and an integration step of 0.01 days, which results in a maximum error of 4 cm s⁻¹ and 1 ppm for the radial velocity and photometric model, respectively (which also takes into account the oversampling factor, Kipping 2010). We included the light-time effect (Irwin 1952), which has an amplitude of ~5 s on the TTVs, corresponding to a displacement of the star by around 2 stellar radii (see Fig. 6). This is small compared to the timing precision of individual transits, and we have checked that the results are not significantly different when the effect is not included, as in Sect. 3.2. The model is parametrised using osculating astrocentric asteroidal orbital elements (Table 2) at the time immediately before the first transit observed by Kepler t_ref = 2 454 958 BJD_TDB, given in Barycentric Dynamical Time. Due to the symmetry of the problem, we fixed the longitude of the ascending node of the interior planet Ω_b at t_ref, we limited the inclination of the outer one i_c < 90°, and rejected models where planet c transits⁵.

Our model has 26 free parameters. In addition to the physical parameters, we considered a radial velocity offset for each instrument with respect to the systemic velocity (assumed to be zero), a global light curve normalisation factor for long- and short-cadence data, and a multiplicative jitter parameter for each data set. A non-informative uniform prior distribution was chosen, and the joint posterior distribution was sampled using the emcee algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013). To minimise correlations, the combinations of parameters listed in Table 2 were used for the Markov chain Monte Carlo (MCMC) algorithm. We ran emcee with 100 walkers from a starting point based on the results of Dawson et al. (2014). We ran 2.4 × 10⁶ steps of the emcee algorithm, and used the last 100 000 steps for the final inference.

Fig. 6 Orbital projections during the Kepler observations for planet b (green), planet c (red), and the star (blue). The origin is the system barycentre, the movement is clockwise, and the orbits are projected in the X-Z reference plane, with the positive Z-axis pointing towards the observer. The two rightmost panels are successive zooms of the panel on the left. A thousand random orbits are drawn form the posterior samples, and the MAP is shown as a black orbit. The red points in the leftmostpanel mark the position of planet c on the MAP orbit at the epochs of transits of planet b, numbered accordingly.

4 Results

In Table 2, we list the maximum a posteriori (MAP) estimate, the median, the 68% credible interval (CI), and the 95% highest density interval (HDI) of the inferred system parameters marginal distributions⁶. The MAP model and credible regions are plotted in Figs. 1 and 2. Figure 7 presents the RV measurements, phase-folded to the best-fit Keplerian period to the MAP photodynamical model.

In Fig. 5, we present the histogram of the marginal distribution of the planet-to-star mass ratios. The distribution for planet c is only slightly affected by the inclusion of the radial velocities. On the other hand, the distribution for planet b is drastically improved by the inclusion of the SOPHIE RVs. By self-consistent modelling of photometry and radial velocity, absolute masses and radii are inferred without resorting to theoretical stellar model. The radii are measured to a precision of 12 and 13% for the star and planet b, respectively. The masses were determined to a precision of 35, 24, and 35% for the star, planet b, and planet c, respectively. The precision on the densities, also derived dynamically, are 19 and 26% for the star and planet b, respectively.

Figure 8 shows the posterior of the TTV, obtained as the mean of the first and fourth contacts interpolated from the sky-projected planet–star separation. The uncertainty in the transit times is not constant. Almenara et al. (2015) showed thatthe TTV posterior distributions are typically narrower around the mid-time of the observations. In addition to this effect, in this case we have a combination of long- and short-cadence transits, giving larger uncertainties for long-cadence epochs 0 to 9 and epoch 21. The uncertainty is reduced for the short-cadence transits (epochs 10 to 20), but there is an increase due to a missed transit at epoch 13. The mean uncertainty of the transit times derived with the photodynamical modelling is two times smaller than any of the timing sets computed on individual transits by Dawson et al. (2014). Clearly, our values rely on the three-body system hypothesis. However, this seems to be a pertinent assumption, as our derived transit times are in agreement with the values presented in Dawson et al. (2014). If more bodies are present in the system, they seem to be dynamically irrelevant at these timescales.

The data favour a model where the planets have low true mutual inclination (i_rel)⁷. We obtained ° (median and 68.3% CI) at t_ref, with the mode at 2.5°, and an upper 99% confidence limit of 23.5°, in agreement with Dawson et al. (2014). These authors reported that the mutual inclination was below 21°, at a 91% confidence level.

Fig. 7 Phase-folded radial velocities for planet b (top) and c (bottom). The points are colour-coded as in Fig. 2 (blue for SOPHIE and purple for HIRES). In the photodynamical model the contribution to RV signal coming from individual planets cannot be separated. Therefore, these figures are constructed by fitting a two-Keplerian model to the photodynamical MAP model radial velocities, and subtracting the individual contribution of each planet from the data. The uncertainties include the MAP jitter value.

Fig. 8 Posterior TTVs of Kepler-419b, computed relative to a linear ephemeris computed using the median values of each epoch during Kepler observations (BJDTDB = 2 454 959.358(23) + 69.728 8(21) × Epoch, where the errors are indicated in parentheses). A thousand random draws from the posterior distribution are used to estimatethe TTV median value and its uncertainty. In the upper panel the median TTV values are labelled with the corresponding epoch number (0 is the first transit observed by Kepler). In the lower panel, the posterior median transit timing value is subtracted to visualise the uncertainty of the distribution. The two sets of transit times derived in Dawson et al. (2014) are shown in light red (TAP) and light blue (GP), slightly offset in the x-direction for clarity. The median transit time was subtracted from each epoch to allow for comparison with our results.

4.1 Secular behaviour

To explore the behaviour of the system at longer timescales, we performed numerical integrations of the system for 10 kyr after the observations. A random sample of size 10 000 from the posterior distribution of the photodynamical modelling was selected as the starting point for the integration. The results for selected parameters are plotted in Fig. A.2. The simulations show that the orbits librate around apsidal anti-alignment with an amplitude of 18°, while the ascending node longitudes librate around alignment with an amplitude of 5°, while the mutual inclination undergoes nutation with an amplitude of 2° about 4° (Mardling 2007, 2010). Furthermore, we put a strong constraint on the mutual inclination of the planets on the longer timescales as well. Our analysis constrains the mutual inclination of the planets to be smaller than 18.9° at 99% confidence level (Fig. A.3).

These results are highly suggestive of gentle relaxation at some point in the system’s history, with the source of dissipation coming from either the protoplanetary disk or tides in planet b (or both). This is, however, apparently not in agreement with the misaligned orbit suggested by the data (Sect. 4.2).

4.2 Stellar models

Although we determined absolute masses and radii dynamically, theoretical stellar models can add further constraints on and improve the precision of some parameters. By doing this, stronger assumptions are introduced in the inference process, which might lead to a degraded accuracy. We interpolated the Dartmouth models (Dotter et al. 2008) using the atmospheric parameters from Dawson et al. (2014; T_eff = 6430 ± 79 K, [Fe/H] = 0.176 ± 0.070), for which we assumed uncorrelated normal distributions (see Fig. 9), and the stellar density from our model, for which we kept only posterior samples that were compatible with the stellar models. In the last column of Table 2 we list the results for the parameters that were improved significantly by this procedure: the semi-major axes, the masses, and the radii. We note that the new posterior distribution of the stellar radius has less mass outside the limit imposed by the non-detection of p-mode oscillations (R_⋆ < 1.9 R_⊙, Dawson et al. 2014) than for the photodynamical determination (R_⋆ = 1.80 ± 0.22 R_⊙).

We obtained an isochronal age of 2.37 ± 0.31 Gyr, which is about half the stellar main sequence lifetime for this mass. We compared this age with that derived via gyrochronology. Coupling the rotational period P_rot = 4.53 ± 0.16 days (Mazeh et al. 2015) with the mass determination using stellar models, we derived a gyrochronological age of 2.59 Gyr (Barnes 2010; Barnes & Kim 2010), where we assumed the zero age main sequence rotational period is between 0.12 and 3.4 days, and we added a systematic 10% error to the statistical error (Meibom et al. 2015). Isochronal and gyrochronological ages agree within the uncertainties.

With the observed v sin i_⋆ = 12.3 ± 1.0 km s⁻¹, and the rotational velocity estimated from the rotational period and the stellar radius, it is possible to derive the inclination of the stellar rotational axis °, in agreement with the more rigorous analysis of Dawson et al. (2014). The orbit of the transiting planet is therefore apparently misaligned with the stellar spin axis. Another star observed to have a coplanar planetary system and a misaligned spin is the red giant Kepler-56 (Huber et al. 2013). While most origin scenarios of spin-orbit misalignment involve misalignment of the orbits themselves (e.g. Fabrycky & Tremaine 2007; Nagasawa et al. 2008), both Kepler-419 and Kepler-56 suggest that the truth is more complex.

Finally, the distance to the system was obtained by modelling the spectral energy distribution of Kepler-419 using the PHOENIX/BT-Settl synthetic spectral library (Allard et al. 2012), and the procedure described in Díaz et al. (2014). Magnitudes from APASS⁸, 2MASS (Skrutskie et al. 2006), and WISE (Wright et al. 2010; Table A.2) were fit to obtain a distance of 993 ± 67 pc, in agreement with the value obtained by Dawson et al. (2012; see Fig. A.4).

Fig. 9 Stellar evolution tracks from Dartmouth models, from the zero age main sequence stage to 25 Gyr age, in the luminosity-effective temperature plane for [Fe/H] = 0.15. The results from the photodynamical modelling using stellar models are denoted by the grey contours corresponding to 39.3%, 86.5%, and 98.9% joint confidence regions. The mass in solar masses is annotated at the beginning of the main sequence of each track.

5 Discussion

With a first detection reported by Holman et al. (2010), the TTVs method is a recent technique. The Kepler satellite continues to be the only facility to have detected them unambiguously. Confirmations are therefore important in order to validate the method, but also to distinguish possibly degenerated TTV predictions or reveal additional companions. The radial velocity technique is the natural method to do this. However, planets detected both with TTVs and radial velocities are rare today. Barros et al. (2014) presented the first radial velocity confirmation of a non-transiting exoplanet discovered by the TTV, and found a mass in agreement with the TTV prediction (Nesvorný et al. 2013). In the Kepler-89 system, Masuda et al. (2013) predicted from TTVs a mass around 50 M_E for Kepler-89d, but Weiss et al. (2013) measured a mass two times greater using radial velocities. It remains unclear where this difference originates.

Non-transiting planets detected by TTVs are also rare (e.g. Nesvorný et al. 2012; Mancini et al. 2016). Radial velocity detections of such planets would be helpful to confirm and refine the TTV predictions. In cases with degenerate predictions (e.g. Ballard et al. 2011), radial velocity measurements would in principle be capable of identifying the correct companion parameters.

In this context, the detection of Kepler-419c in SOPHIE radial velocities is only the second confirmation of a non-transiting exoplanet discovered by TTVs. While Dawson et al. (2014) presented the high-precision radial velocity values of this system, they were not sufficient to detect the planet independently, as we did (see Sect. 3.1). The orbital period and velocity amplitude of Kepler-419c measured with SOPHIE agree with the values predicted from the modelling of the photometry alone. This indicates that the model hypotheses are relevant, and that no additional companions significantly perturb the two detected planets. The dynamical analysis described in Sect. 5.2 also points to this conclusion. Furthermore, we did not detect the signal of additional planets in the residuals of the radial velocity data.

The difference between the radial velocities from the n-body model and those obtained using non-interacting Keplerian curves is relatively small, except close to the periastron passage of planet b where the difference can be up to m s⁻¹ (see Fig. 10). The difference increases as the time span of the observations increases. Here we used the combined time span of HIRES and SOPHIE, i.e. 2.3 years.

The main results presented here come from the combined analysis of the Kepler photometry and the SOPHIE radial velocities using a photodynamical model. Under very simple and general assumptions, the model takes into consideration the mutual gravitational interactions of the planets in a consistent manner. A novel result is the dynamic determination of absolute masses and radii. Our results depend on the validity of Newtonian mechanics alone, as well as some simple model assumptions: sphericity of all bodies in the system, the number of objects in the system, and the chosen limb-darkening law. To date, Kepler-419 is the system with the most accurate masses and radii determined using the photodynamical model for a multiple planetary system with a single host star⁹. While the parameters are determined less precisely than when using theoretical atmosphere and evolutionary stellar models, the limiting factor is certainly the inherent stellar velocity jitter, and to a lesser extent to the photon noise obtained with SOPHIE for this relatively faint star. It is to be expected that more precise results will be achieved for brighter quieter stars, given the same amount of dynamical information.

Under this model, the radius and mass of planet b are constrained to 13% and 24%, respectively. At the same time, the obtained precision on the stellar mass is 35%. That is, the mass of planet b is measured more precisely than the mass of its stellar host. This shows that the constraints provided by our modelling do not imply a measurement relative to the star, as is provided by the Keplerian model of radial velocities. Also, it provides a method for characterising dynamically interacting planets without being limited by knowledge of the star. Ultimately, this technique can be used to test theoretical stellar models, especially for quiet stars.

Fig. 10 Difference between the stellar radial velocity (computed with the n-body integration) and a two-Keplerian fit to the model prediction, computed over 100 random MCMC samples from the full photodynamic model. The black line and the grey shaded region represent the median and the 68.3% confidence interval, respectively. The smaller panel is an enlargementaround the maximum difference, marked in red in the main panel.

5.1 Masscontraints on Kepler-419 b from photometry

Dawson et al. (2014) affirmed that the TTV of the inner planet are “not at all sensitive” to its mass, and that its mass constraint comes exclusively from the RV data. Our analysis with photometry only (Sect. 3.2) shows that an upper limit can be set on the mass ratio of planet b (Fig. 5). In fact, this non-zero signal is a result of the proximity of the period ratio (which is around 9.7) to 9 and 10, together with the significant eccentricities of both orbits, as we show next.

Figure 11 shows the Lomb–Scargle periodogram (plotted against frequency ν in units of the inner orbital frequency ν_b) of the 22-transit signal shown in Fig. 8. Only the Nyquist window 0 < ν∕ν_b ≤ 1∕2 is shown; higher and negative frequencies are aliases of these values and are such that if is the power at scaled frequency , then , where k is an integer¹⁰. While it is true that the TTV signals of near-circular (single-star) planetary systems are dominated by harmonics associated with first- and second-order resonances (and hence require period ratios of around 3 or less to be detected), the TTVs associated with the 22 transits of Kepler-419b are dominated by harmonics with frequencies , n′ = 1, 2, 10, where ν_b and ν_c are the orbital frequencies of planets b and c, respectively. The substantial power in n′ = 10 is due to both eccentricities (see Table 3).

To verify this, we simulated a light curve with 54 transits of planet b using the MAP values in Table 2 with a sampling and white noise amplitude equal to that of the Kepler SC data. The TTV periodogram for the simulated LC is also shown in Fig. 11. It reveals that the dominant harmonic is n′ = 10, and that the powerful n′ = 9 harmonic is not resolved in the 22-transit data set¹¹. This lack of resolution results in a poor constraint on the mass of the transiting planet for the following reasons.

For a coplanar system there are seven unknowns (the masses of the two planets, the two eccentricities and corresponding longitudes of periastron, and the mean longitude of planet c at epoch). Therefore, information from at least four well-resolved harmonics (each with an amplitude and a phase) is needed in order to place reasonable bounds on all seven parameters from TTV data alone. We can show analytically (Mardling, in prep.) that the back effect of the mass of planet b on its own TTVs is to reduce the amplitude of the n′ th harmonic by a factor of approximately 1 − n′(M_b∕M_*). Thus, for the low-order harmonics n′ = 1, 2, 3, typically themost powerful detectable for low-eccentricity systems, the difference is of the order of a fraction of a per cent. On the other hand, for n′ = 9 and n′ = 10 the effect is at the 2% level, and can therefore be detected, provided they can be correctly resolved.

Table 3 lists the normalised Lomb–Scargle power, , for n′ = {10, 9, 1, 2, 3}, for values of the inner and outer eccentricities around the MAP values for 54 transits. The dominant harmonic is highlighted in bold. We note that substantial power is associated with all five harmonics in the observed system (case 1), in contrast with cases 2 and 3 for which the eccentricity of planet c is, respectively, reduced and increased by around 0.1. Similarly, changing the eccentricity of the transiting planet significantly affects the distribution of power in the various harmonics. Thus, we see that the Kepler-419 system is in a serendipitous sweet spot for the determination of system parameters via TTVs (given a sufficient number of transits). This is confirmed by the filled grey histograms in Fig. 5, showing the posterior sample the masses of planets b (left) and c (right), based on the photodynamical analysis of the simulated light curve with 54 transits (and no radial velocities). Not only is the mass of the transiting planet well resolved, the mass of the perturbing planet is also significantly better resolved than it is with just 22 transits, with or without radial velocities.

Fig. 11 Lomb–Scargle periodogram of the TTVs for 22 transits (black curve) and 54 transits (red curve), with dominant frequencies labelled (see text). Normalised power is plotted against frequency in units of the orbital frequency νb of planet b. Only the Nyquist window 0 ≤ ν∕νb ≤ 1∕2 is shown; higher frequencies are aliases of these values.

5.2 Origin of the Kepler-419 system

The coplanar yet highly eccentric nature of Kepler-419 presents a puzzle regarding its origin. The system’s proximity to a stable fixed-point suggests that it is likely to have been brought gently to this relaxed state via disk dissipation or planet (or even stellar) tides or both.

While the existence of low-eccentricity fixed points in coplanar systems is well known (e.g. Mardling 2007), as is the existence of high-eccentricity fixed points in non-coplanar systems (those associated with the Kozai mechanism; e.g. Naoz 2016), the existence of high-eccentricity fixed points in coplanar systems (e.g. Nagasawa et al. 2003) seems to be less well known. Figure 12 shows the existence of low- and high-eccentricity fixed points, with curves produced by integrating the secular equations of motion (in the absence of damping), with the Kepler-419 system parameters as initial conditions except for e_b and ϖ_b − ϖ_c, which are indicated by red dots. In panel a, e_c is held constant, while in panel b, e_c is allowed tovary. The blue cross shows the position of Kepler-419 and its proximity to the fixed point, suggesting that it is highly likely that the system has relaxed towards that state.

In orderto reach planet b’s high eccentricity, a significant source of external torque must be identified, which raises the eccentricity and also maintains the coplanarity of the planets. In the following we discuss the viability of various mechanisms capable of providing such a torque.

Fig. 12 Variation of eb as a functionof ϖb − ϖc in the absence of damping, showing the existence of low- and high-eccentricity fixed points. Curves produced by integrating the secular equations of motion with the Kepler-419 system parameters as initial conditions except for eb and ϖb − ϖc, which are indicated by red dots. In panel a, ec is held constant, while in panel b, ec is allowed to vary. Phase curves which cross have different values of ec at those points. The blue cross shows the position of Kepler-419 and its proximity to the fixed point. We note the different y-axis scales.

5.2.1 External forcing

Dawson et al. (2014) considered Kozai forcing of planet b by planet c, but dismissed it because the mutual inclination (here constrained even further to be no more than a few degrees) is not compatible with the large value normally necessary to achieve eccentricity growth. A potential solution to this quandary is that a significantly inclined external fourth body has brought the system to its current relaxed state. Takeda et al. (2008) showed that under favourable circumstances a pair of planets will respond to the torque from such a body and will remain almost coplanar while undergoing Kozai oscillations. The key to this mechanism is that the resulting precession rates of each planetary orbit should be similar so that their response is in concert, with the difference in the node angles librating around a fixed value. However, this cannot be the case for the Kepler-419 system because the mutual torques are already almost perfectly balanced; the addition of a fourth body would destroy this harmony.

Dawson et al. (2014) also dismissed the “orbit flip” mechanism (Li et al. 2014) because they did not observe this to happen in integrations using their observed orbital parameters. This mechanism initially involves almost coplanar orbits periodically undergoing large excursions in relative inclination and (inner) eccentricity, with periodic 180° flips in the orientation of the orbital angular momentum vector of a body relative to that of an external perturber. We note that the Li et al. (2014) analysis is done in the test particle approximation and therefore does not involve any variations in the perturber’s orbit. It is associated with octopole-level terms in the disturbing function¹², and occurs for favourable orbital configurations which satisfy an analytic condition involving the relative strengths of the octopole and quadrupole terms. In fact, we can study the problem from the point of view of finding the fixed points of the octopole-level equations of motion for the general problem (all three bodies massive) and enquire about their stability (Mardling, in prep.). We find that flip-type solutions are associated with unstable fixed points, and in the presence of damping a system would evolve away from such a state. It therefore seems unlikely that a real system would be observed in a “flip” configuration, and it is clear that the Kepler-419 system is not in this state.

5.2.2 Spin-orbit coupling

Another promising mechanism is spin-orbit coupling between the planet and effectively both orbits. Correia et al. (2012) showed that under favourable circumstances, the eccentricity of the orbit of a spinning planet in a two-planet system can be raised to high values. This mechanism relies on the fact that the spin-synchronisation timescale of a planet is much shorter than the eccentricity-damping timescale, and so is effective for relatively long-period orbits. Since a companion will modulate the eccentricity of the tidally active planet on the secular timescale, and since the spin will tend to synchronise with the orbital frequency at periastron, but will lag behind the eccentricity forcing by an amount that depends on the spin-synchronisation timescale, a small positive drift in the average eccentricity can result if the two timescales are similar. The average eccentricity continues to increase until torques are balanced (the system reaches a fixed point). However, tidal heating at periastron also increases (because secular forcing does not change the semi-major axis and hence decreases the periastron separation as the eccentricity increases), ultimately shrinking the semi-major axis and circularizing the orbit on a timescale which may be longer than the age of the system (see Fig. 1 in Correia et al. 2012). Using a double-averaged code with spin-orbit coupling, relativistic apsidal advance and tides (Mardling & Lin 2002), with realistic structure and damping parameters and the current observed parameters of the Kepler-419 system but with low initial values of the inner eccentricity and arbitrary initial values of the difference in the longitudes of periastron ϖ_b − ϖ_c, some positive drift in the latter is observed for some initial configurations but appears to be quite sensitive to the initial value of ϖ_b − ϖ_c. None of the systems considered experienced eccentricity increases of much more than 0.1 on the several Gyr timescale; however, if some gentle disk migration (weak enough to avoid resonance capture) is included, thereby providing additional torques, as well as the evolution of the planet radius due to gravitational contraction (resulting in stronger tides in the past), it is conceivable that the mechanism could produce the Kepler-419 system as it is currently observed.

5.2.3 Collision and damping

Finally we consider a possible collision and subsequent damping scenario for the origin of Kepler-419, and examine the nature of the high-eccentricity fixed point and its implications for the existence of hot Jupiters.

Consider the scenario in which the Kepler-419 system forms a three-planet 1 : 2 : 4 Laplace configuration via inward convergent migration of three (or initially more) planets, with the current planet c on the outside. Notwithstanding the protective nature of the Laplace resonance, eccentricity growth can result in orbit crossing and eventual collision of the inner two planets for favourable systems and disk conditions. As far as we are aware, fully self-consistent modelling ofsuch an event has not been done, and as such it remains unclear how angular momentum would be distributed between the resultant planet orbits and the ejecta following collision. It is conceivable that the post-collision orbit could have considerablyless angular momentum than that of the precursor orbits, with the excess being returned to the protoplanetary disk or escaping the system.

A less constraining scenario is for the pre-cursor system to consist of a pair of planets in a 2:1 resonance interior to the current planet c (but not participating in a Laplace-type configuration), which again undergo collision following eccentricity pumping and subsequent instability (Goldreich & Schlichting 2014). Either way, we propose a scenario where a collision results in a highly eccentric orbit whose apsidal orientation with respect to that of planet c is such that it is immediatelyor subsequently captured by the high-eccentricity fixed point associated with coplanar two-planet systems.

Depending on the initial value of the angle between the apsidal lines, such a state can in fact prevent a highly eccentric system fromcircularising, as appears to be the case for Kepler-419. Alternatively, it can result in a hot Jupiter, as the following demonstrates.

5.2.4 Formation of hot Jupiters

When damping due to planetary tides is included, the semi-major axis of planet b responds by shrinking at a rate which depends on the minimum value of the periastron separation over a secular cycle and the time spent near that value. Panels a–d of Fig. 13 show the evolution of the eccentricity and semi-major axis for two systems which are identical except for the initial value of ϖ_b − ϖ_c. For panels a and b (Case I), ϖ_b(0) − ϖ_c(0) = 180°, while ϖ_b(0) − ϖ_c(0) = 90° for panels c and d (Case II). The remaining system parameters are as for Kepler-419, except that a_b (0) = 0.4 and e_b (0) = 0.93. The Q-value and Love number of the planet are taken as 10⁵ and 0.3, respectively. The first 10⁸ years of evolution is shown. Very little orbit shrinkage and circularisation has occurred for Case I after 10⁸ years, with the system remaining trapped near the fixed point at (ϖ_b − ϖ_c = π, e_b = 0.83) with an average periastron separation of 0.06 au. In contrast, Case II initially achieves a maximum eccentricity of 0.96 during a short period of libration, after which it escapes the librating region with a non-oscillatory eccentricity, and thus a permanently low periastron separation. With an average periastron separation of 0.02 au, circularisation is rapid. We emphasise that the only difference between the two cases is the inital value of the difference in the apsidal longitudes.

Panels e and f show the dependence of the values of e_b and a_b at 10⁸ yr, 10⁹ yr, and 2.3 × 10⁹ yr (the estimated age of the system) on the initial difference in the apsidal longitudes, suggesting that the long-term state of a system like Kepler-419 is highly dependent on the value of ϖ_b − ϖ_c at the time it was brought to that state (via collision or some other mechanism). For the adopted Q-value the system has circularised after 2.3 × 10⁹ yr, suggesting that the true Q-value is higher or has varied over the lifetime of the system, or that the system is younger, or both. Panel g shows the evolution of the periastron distance of planet b for Cases I and II, while panel h shows its average value, ⟨p_b ⟩, as a function of the initial difference in the apsidal longitudes. Thus while a system is trapped in the libration state it tends to maintain a high eccentricity because ⟨p_b⟩ is relatively high, while escape from libration is associated with a permanently low value of p_b.

This suggests that a system like HAT-P-13 (Bakos et al. 2009), whose periods and planet masses are 2.9 d and 428.5 d, and 0.85 M_J and 15.2 M_J, respectively, and whose star has a mass of 1.22 M_⊙, may have had a similar origin to Kepler-419, but found itself with a value of ϖ_b − ϖ_c conducive to circularisation. Just as for the Kozai migration mechanism, such an evolutionary history would leave little room for intervening planets. However, unlike the Kozai mechanism, the present mechanism produces hot Jupiters without the need of mutual inclination.

Fig. 13 Effect of the initial value of ϖb − ϖc on the subsequent tidal evolution of a system. Panels a and b: ϖb(0) − ϖc(0) = 180° (Case I); panels c and d: ϖb(0) − ϖc(0) = 90° (Case II). The red dot and blue cross in panel a indicate the initial and current state, respectively, of Kepler-419. Panels e and f: valuesof ab and eb at tend = 108 yr (solid black curves), tend = 109 yr (dashed blue curves), and tend = 2.3 × 109 yr (the estimated age of the system; dot-dashed red curves) as functions of the initial value of ϖb − ϖc, which show that for a range of initial relative orientations, planet b can sustain a high value of eb because it spends most of its time with a periastron distance which is too high for significant circularisation to occur. Panel g: evolution of planet b’s periastron distance, pb, for Case I (blue) and Case II (red), the latter showing escape from libration and hence a permanently low value of pb. Panel h: average value of pb, ⟨pb ⟩, as a functionof ϖb(0) − ϖc(0). Angles are in degrees.

Acknowledgements

We thank Rebekah I. Dawson and Geoffrey W. Marcy for helpful discussions. We thank the technical team at the Observatoire de Haute-Provence for their support with the SOPHIE instrument and the 1.93 m telescope and, in particular, for the essential work of the night assistants. Financial support for the SOPHIE observations from the Programme National de Planétologie (PNP) of CNRS/INSU, France, is gratefully acknowledged. We also acknowledge support from the French National Research Agency (ANR-08- JCJC-0102-01). We thank L. Kreidberg for her Mandel & Agol code and Y. Revaz for his assistance with the computing cluster used in this work. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Simulations in this paper made use of the rebound code which can be downloaded freely at http://github.com/hannorein/rebound. These simulations have been run on the Regor cluster kindly provided by the Observatoire de Genève. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange, and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. This work has been carried out within the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. S.C.C.B. acknowledges support by the Fundação para a Ciência e a Tecnologia (FCT) through the Investigador FCT Contract IF/01312/2014/CP1215/CT0004 and also acknowledges support from FCT through national funds and by FEDER through COMPETE2020 by these grants UID/FIS/04434/2013 & POCI-01-0145-FEDER-007672 and PTDC/FIS-AST/1526/2014 & POCI-01-0145-FEDER-016886. A.S.B. acknowledges funding from the European Union Seventh Framework programme (FP7/2007-2013) under grant agreement No. 313014 (ETAEARTH).

Appendix A: Additional figures and tables

Fig. A.1 Two-parameter joint posterior distributions for the most relevant MCMC model parameters. The 39.3, 86.5, and 98.9% two-variable joint confidence regions are denoted by three different grey levels; in the case of a Gaussian posterior, these regions project on to the one-dimensional 1, 2, and 3 σ intervals. The histogram of the marginal distribution for each parameter is shown at the top of each column, except for the parameter on the last line, which is shown at the end of the line. Units are the same as in Table 2.

Fig. A.2 Evolution of the Jacobi orbital elements (ϖ = Ω + ω, is the longitude of the periapsis) over 10 kyr since the beginning of Kepler observations. 10 000 random draws from the posterior distribution are shown. The grey points correspond to the integration using the MAP values.

Fig. A.3 Posterior probability function of the mutual inclination of the planets in Jacobi coordinates, based on integrations over 10 kyr.The shaded regions correspond to the 68.3% [1.20, 7.58], 95.5% [0.26, 14.17], and 99.7% [0.05, 22.02] highest density intervals.

Fig. A.4 Spectral energy distribution of Kepler-419 based on magnitudes from the literature (red circles; see text for details). The best-fit spectrum is plotted as a solid black curve, and the integrated fluxes in the photometric bands are plotted as open circles. The residuals are given in the bottom panel.

References

All Tables

All Figures

Fig. 1 Transits of Kepler-419b observed by Kepler. Each panel is centred at the linear ephemeris (indicated by the vertical grey lines, and reported in the caption of Fig. 8). For short-cadence data, 29.4-minute binned data is shown in addition to the observed data points. Each panel is labelled with the epoch; zero is the first transit after tref. The black curve is the median oversampled model over 10 000 random MCMC steps. In the lower part of each panel the residuals after subtracting the MAP model to the observed data are shown. The shades of grey represent the 68.3, 95.5, and 99.7% credible intervals, which are hardly distinguishable in the residual panels and show an increased uncertainty at ingress and egress times.

In the text

	Fig. 2 Radial velocities of Kepler-419 observed by HIRES (upper plot, purple points) and SOPHIE (lower plot, blue points). Theblack curve is the median model from the photodynamical analysis, with the 68.3, 95.5, and 99.7% credible curves.
In the text

	Fig. 3 Generalised Lomb–Scargle periodogram of the bisector velocity span of the SOPHIE observations. A clear peak is seen at P = 4.58 days, in agreement with the rotational period determination by Mazeh et al. (2015), whose 1 − σ range is shown in grey.
In the text

Fig. 4 Generalised Lomb–Scargle periodograms of the radial velocities. The green and red vertical lines mark the period of planet b and c, respectively. The grey horizontal lines represent the 50, 10, and 1% p-value levels. From top to bottom: HIRES and SOPHIE data; HIRES velocities; SOPHIE velocities; SOPHIE after substraction of Keplerian curve at the period of planet b; SOPHIE velocities after substraction of a two-Keplerian model with the periods of planets b and c. The periodogram plots are provided by DACE (dace.unige.ch).

In the text

	Fig. 5 Posteriors of planet-to-star mass ratios for planet b (left) and planet c (right). The black histograms correspond to the analysis of the photometric data alone, while the grey histogram includes radial velocities. The filled grey histograms show the PDF for the masses of planets b and c, respectively, for the analysis of the simulated light curve with 54 transits and no radial velocities.
In the text

Fig. 6 Orbital projections during the Kepler observations for planet b (green), planet c (red), and the star (blue). The origin is the system barycentre, the movement is clockwise, and the orbits are projected in the X-Z reference plane, with the positive Z-axis pointing towards the observer. The two rightmost panels are successive zooms of the panel on the left. A thousand random orbits are drawn form the posterior samples, and the MAP is shown as a black orbit. The red points in the leftmostpanel mark the position of planet c on the MAP orbit at the epochs of transits of planet b, numbered accordingly.

In the text

Fig. 7 Phase-folded radial velocities for planet b (top) and c (bottom). The points are colour-coded as in Fig. 2 (blue for SOPHIE and purple for HIRES). In the photodynamical model the contribution to RV signal coming from individual planets cannot be separated. Therefore, these figures are constructed by fitting a two-Keplerian model to the photodynamical MAP model radial velocities, and subtracting the individual contribution of each planet from the data. The uncertainties include the MAP jitter value.

In the text

Fig. 8 Posterior TTVs of Kepler-419b, computed relative to a linear ephemeris computed using the median values of each epoch during Kepler observations (BJDTDB = 2 454 959.358(23) + 69.728 8(21) × Epoch, where the errors are indicated in parentheses). A thousand random draws from the posterior distribution are used to estimatethe TTV median value and its uncertainty. In the upper panel the median TTV values are labelled with the corresponding epoch number (0 is the first transit observed by Kepler). In the lower panel, the posterior median transit timing value is subtracted to visualise the uncertainty of the distribution. The two sets of transit times derived in Dawson et al. (2014) are shown in light red (TAP) and light blue (GP), slightly offset in the x-direction for clarity. The median transit time was subtracted from each epoch to allow for comparison with our results.

In the text

Fig. 9 Stellar evolution tracks from Dartmouth models, from the zero age main sequence stage to 25 Gyr age, in the luminosity-effective temperature plane for [Fe/H] = 0.15. The results from the photodynamical modelling using stellar models are denoted by the grey contours corresponding to 39.3%, 86.5%, and 98.9% joint confidence regions. The mass in solar masses is annotated at the beginning of the main sequence of each track.

In the text

Fig. 10 Difference between the stellar radial velocity (computed with the n-body integration) and a two-Keplerian fit to the model prediction, computed over 100 random MCMC samples from the full photodynamic model. The black line and the grey shaded region represent the median and the 68.3% confidence interval, respectively. The smaller panel is an enlargementaround the maximum difference, marked in red in the main panel.

In the text

	Fig. 11 Lomb–Scargle periodogram of the TTVs for 22 transits (black curve) and 54 transits (red curve), with dominant frequencies labelled (see text). Normalised power is plotted against frequency in units of the orbital frequency νb of planet b. Only the Nyquist window 0 ≤ ν∕νb ≤ 1∕2 is shown; higher frequencies are aliases of these values.
In the text

Fig. 12 Variation of eb as a functionof ϖb − ϖc in the absence of damping, showing the existence of low- and high-eccentricity fixed points. Curves produced by integrating the secular equations of motion with the Kepler-419 system parameters as initial conditions except for eb and ϖb − ϖc, which are indicated by red dots. In panel a, ec is held constant, while in panel b, ec is allowed to vary. Phase curves which cross have different values of ec at those points. The blue cross shows the position of Kepler-419 and its proximity to the fixed point. We note the different y-axis scales.

In the text

Fig. 13 Effect of the initial value of ϖb − ϖc on the subsequent tidal evolution of a system. Panels a and b: ϖb(0) − ϖc(0) = 180° (Case I); panels c and d: ϖb(0) − ϖc(0) = 90° (Case II). The red dot and blue cross in panel a indicate the initial and current state, respectively, of Kepler-419. Panels e and f: valuesof ab and eb at tend = 108 yr (solid black curves), tend = 109 yr (dashed blue curves), and tend = 2.3 × 109 yr (the estimated age of the system; dot-dashed red curves) as functions of the initial value of ϖb − ϖc, which show that for a range of initial relative orientations, planet b can sustain a high value of eb because it spends most of its time with a periastron distance which is too high for significant circularisation to occur. Panel g: evolution of planet b’s periastron distance, pb, for Case I (blue) and Case II (red), the latter showing escape from libration and hence a permanently low value of pb. Panel h: average value of pb, ⟨pb ⟩, as a functionof ϖb(0) − ϖc(0). Angles are in degrees.

In the text

Fig. A.1 Two-parameter joint posterior distributions for the most relevant MCMC model parameters. The 39.3, 86.5, and 98.9% two-variable joint confidence regions are denoted by three different grey levels; in the case of a Gaussian posterior, these regions project on to the one-dimensional 1, 2, and 3 σ intervals. The histogram of the marginal distribution for each parameter is shown at the top of each column, except for the parameter on the last line, which is shown at the end of the line. Units are the same as in Table 2.

In the text

	Fig. A.2 Evolution of the Jacobi orbital elements (ϖ = Ω + ω, is the longitude of the periapsis) over 10 kyr since the beginning of Kepler observations. 10 000 random draws from the posterior distribution are shown. The grey points correspond to the integration using the MAP values.
In the text

	Fig. A.3 Posterior probability function of the mutual inclination of the planets in Jacobi coordinates, based on integrations over 10 kyr.The shaded regions correspond to the 68.3% [1.20, 7.58], 95.5% [0.26, 14.17], and 99.7% [0.05, 22.02] highest density intervals.
In the text

	Fig. A.4 Spectral energy distribution of Kepler-419 based on magnitudes from the literature (red circles; see text for details). The best-fit spectrum is plotted as a solid black curve, and the integrated fluxes in the photometric bands are plotted as open circles. The residuals are given in the bottom panel.
In the text