The star K2-19 hosts a pair of Neptunian planets deep inside the 3:2 resonance. They induce strong transit-timing variations with two incommensurate frequencies. Previous photodynamical modeling of 3.3 years of transit and radial velocity data produced mass estimates of \(32.4\pm 1.7\,M_\oplus\) and \(10.8\pm 0.6\,M_\oplus\) for planets b and c, respectively, and corresponding eccentricity estimates of \(0.20\pm 0.03\) and \(0.21\pm 0.03\). These high eccentricities raise questions about the formation origin of the system, and this motivated us to extend the observing baseline in an attempt to better constrain their values. We present a photodynamical analysis of 10 years of transit data that confirms the previous mass estimates (\(30.8\pm 1.3\,M_\oplus\) and \(11.1\pm 0.4\,M_\oplus\)), but reduces the median eccentricities to \(0.04\pm 0.02\) and \(0.07\pm 0.02\) for b and c, respectively. These values are more consistent with standard formation models, but still involve nonzero free eccentricity. The previously reported high eccentricities appear to be due to a single transit for which measurements taken at twilight mimicked ingress. This resulted in a 12-minute error in the midtransit time. The data that covered 1.3 and 5 so-called super and resonant periods were used to match a Fourier analysis of the transit-timing variation signal with simple analytic expressions for the frequencies and amplitudes to obtain planet mass estimates within 2% of the median photodynamical values, regardless of the eccentricities. Theoretical details of the analysis are presented in a companion paper. Additionally, we identified a possible planet candidate situated exterior to the b–c pair. Finally, in contrast to a previous study, our internal structure modeling of K2-19 b yields a metal mass fraction that is consistent with core accretion.

1 Introduction↩︎

The migration process is a complex mechanism that takes place in protoplanetary disks. It can result in the capture of planets in mean motion resonance (MMR) chains [1]. When the protoplanetary disk dissipates, most of the resonant chains are expected to break [2], [3] because most systems are observed outside the MMR [4]. It is particularly interesting to study the few remaining resonant chains because it can help us to understand the early life of planetary systems. K2-19 [5]–[7] is one of these systems. It is composed of three known transiting planets: an inner Earth-size planet [8], and an outer pair of Neptunian planets (K2-19 b and c) inside the 3:2 MMR (with periods 7.9 and 11.9 days) that present significant transit-timing variations [9], [10]. These TTVs allowed the first measurements of the masses of K2-19 b and c [11], [12]. The latest characterization of the system by [13], who used a mix of TTVs and a photodynamical modeling [14], found eccentricities of \(0.20 \pm 0.03\) and \(0.21 \pm 0.03\) for K2-19 b and c, respectively. Disk migration does not favor significant orbital eccentricity excitation, however, because planet-disk interactions tend to damp eccentricities [15]. [16] proposed secular angular momentum exchange with the interior planet d at 2.5 days, as well as an additional planet in the system. When this additional planet is included, the eccentricities of K2-19 b and c are reduced.

The precision of the radial velocities in the K2-19 system is moderate because the target is faint, and significant variability is caused by stellar activity [13], [17], [18]. Therefore, transit photometry (notably through TTVs) remains the main technique for investigating this system, and it alone currently has the sensitivity to accurately measure orbital eccentricities and mutual inclinations. These measurements might help us to unveil the formation history of the system.

Motivated by the work of [16], we obtained additional transit observations of K2-19 b and c with several facilities from the ground and space. They extend the observing baseline of [13] from 3 to 10 years and cover more than a complete superperiod [19]. With 18 and 9 additional transits of planets b and c, respectively, we find that the system is considerably closer to the stationary state of the 3:2 resonance (see Fig. 8). Our posterior included nonaligned systems with very little free eccentricity (consistent with the standard model of disk migration) as well as aligned systems with a free eccentricity no higher than 0.1. The forced eccentricities of planets b and c are about 0.015 and 0.045 for all samples in the posterior (see Sect. 5.3).

The article is organized as follows: We present the new transit photometry of K2-19 b and c in Sect. 2. We revise the stellar parameters in Sect. 3. We describe our analysis in Sect. 4, and we present our results in Sect. 5. Finally, we discuss our results in Sect. 6.

2 Observations↩︎

5pt

Table 1 summarizes all transit observations we used. We detail the new observations by telescope below together with the revision of two datasets from the literature.

2.1 K2↩︎

K2-19 was discovered during Campaign 1 of the K2 mission of the Kepler space telescope [20]. The light curve presents significant systematics because the fine-pointing capability of the spacecraft was lost. We used EVEREST[21], [22] to correct for these systematics by masking the transits of the known transiting planets (including candidate K2-19 e; Sect. 4.1) to minimize changes in the transits during the detrending process.

2.2 LCOGT↩︎

[13] presented one partial transit of K2-19 b observed with the Las Cumbres Observatory Global Telescope [23] on June 25, 2017², using a 1 m telescope at the South African Astronomical Observatory. We reanalyzed these transit data because we found that the transit timing reported by [13] deviated from the model posterior (Sect. 4.2) by \(3.5\,\sigma\) (equivalent to \(12.4\pm3.6\) minutes³). We note that the observations they interpreted as the transit ingress were acquired during twilight. We performed photometry with prose[24] on images calibrated with the standard LCOGT pipeline [25].

2.3 ExTrA↩︎

We started our follow-up of K2-19 b and c in 2021 with the facility called Exoplanets in Transits and their Atmospheres [26], which is located at La Silla Observatory in Chile. ExTrA is a low-resolution near-infrared (0.85 to 1.55 \(\mu\)m) multi-object spectrograph fed by three 60 cm telescopes. We used 8\(\arcsec\) aperture fibers and the lowest-resolution mode (\(R\) \(\sim\)​20) of the spectrograph with an exposure time of 60 seconds. During the 2021 and 2022 seasons, we targeted several transit windows of K2-19 b and c using the transit forecast from [13], which had uncertainties between 35 and 65 minutes for K2-19 b and between 2 and 4 hours for K2-19 c. We were able to observe one full transit and two partial transits of K2-19 b and no transit of K2-19 c. The observed transits of K2-19 b differed between \(2.5\,\sigma\) to \(4.5\,\sigma\) from the predictions by [13]. From season 2023 onward, we used a TTV modeling of the available timings to forecast transits for ExTrA and other facilities. As K2-19 b and c have detected gravitational interactions, the observations of K2-19 b also improved the transit forecast of K2-19 c. This allowed us to finally observe a transit of K2-19 c on April 19, 2023 with ExTrA, which to our knowledge is the first time that the transit of this 6-hour TTV planet has been observed from the ground.

2.4 TESS↩︎

The Transiting Exoplanet Survey Satellite [27] observed K2-19 in sectors 45, 46, and 72 with a cadence of 2 minutes. For the analysis, we used the presearch data-conditioning simple aperture photometry (PDCSAP; [28]; [29], [30]) light curve of K2-19 that was produced by the TESS Science Processing Operations Center [31].

2.5 CHEOPS↩︎

The CHaracterising ExOPlanet Satellite [32] observed two transits of K2-19 b on February 29, 2024 and May 10, 2024, and one transit of K2-19 c on March 29, 2024 (Table 2). The raw data were automatically processed by the CHEOPS data reduction pipeline [33]. We performed point-spread function photometry with the PIPE package⁴ [34]. Unfortunately, even with an observation efficiency of 93.7%⁵, no first contact of the transit of K2-19 c was observed.

2.6 EulerCam↩︎

We observed two transits of K2-19 b on March 8 and 16, 2024, and one transit of K2-19 c on May 4, 2024 with the EulerCam (ECAM) at the Swiss 1.2 m Euler telescope located at La Silla Observatory. The image reduction was carried out following [35] and [36]. Aperture photometry was performed using prose[24], which relies on astropy[37]–[39] and photutils[40]. The optimal differential photometry followed [41].

3 Stellar parameters↩︎

To determine the stellar parameters of K2-19, we fit its spectral energy distribution (SED) using stellar atmosphere and evolution models, along with the parallax determination from Gaia [42], [43]. We constructed the SED of K2-19 using the magnitudes from Gaia data release 3 [44], the 2-Micron All-Sky Survey [45], [46], and the Wide-field Infrared Survey Explorer [47], [48]. The measurements are listed in Table 3. We used the PHOENIX/BT-Settl [49] stellar atmosphere model and the two stellar evolution models Dartmouth [50] and PARSEC[51]. We modeled the SED using the procedure described by [52], with informative priors for the effective temperature (\(T_{\mathrm{eff}}\)), surface gravity (\(\log g\)), and metallicity (\([\rm{Fe/H}]_\star\)) from [53], and for the distance from the Gaia zeropoint-corrected parallax [54]. We added an uncertainty of 2% quadratically to the \(T_{\mathrm{eff}}\) following [55]. We used uninformative priors for the remaining parameters. We used a jitter [56] for each set of photometric bands (Gaia, 2MASS, and WISE). The offset in the posterior median between the two tested stellar evolution models is 0.008 \(R_{\odot}\) and 0.017 \(M_{\odot}\). We merged the results from the two different stellar evolution models assuming an equal probability for each. The parameters priors and posteriors are listed in Table 6. The data with the maximum a posteriori (MAP) stellar atmosphere model are shown in Fig. 11. For the stellar radius, we used a systematic uncertainty floor of 4% following [55] (we added an uncertainty of 4% quadratically to the radius estimated with the SED analysis). We confirmed the mass values with the four stellar evolution models implemented by kiauhoku[57] and determined a mean value of 0.882 \(M_{\odot}\), a standard deviation between models of 0.012 \(M_{\odot}\), and a maximum offset of 0.026 \(M_{\odot}\). Our mass value and uncertainty already encompass the differences between the models in kiauhoku. We compare the stellar parameters derived in this section (Table 3) with literature values listed in Table 7. The stellar mass and radius were used as priors in the photodynamical modeling (Sect. 4.2).

We used the stellar mass and the stellar rotation period determined in Sect. 4.2 (\(P_{\rm rot} = 18.36 \pm 0.51\) days) to derive a gyrochronological age, neglecting the influence of the planets. Based on the formulation of [58] and [59] with initial periods P\(_0\) between 0.12 and 3.4 days, the age is \(1.70 \pm 0.25\) Gyr [60].

4pt

4 Analysis↩︎

4.1 Search for additional transiting planets↩︎

The K2 light curve of K2-19 is dominated by rotational modulation, with an amplitude variability (peak to peak) of about 1.2% and a periodicity of \(18.36 \pm 0.51\) days (Sect. 4.2). We searched for additional transiting planets in the K2 data using nuance[64], a transit search algorithm that accounts for correlated noise (stellar variability in this case), which we modeled with the rotation kernel Gaussian process (GP).

We found the already known planets with an S/N of 54 (K2-19 b), 33 (K2-19 c), and 13 (K2-19 d). Then, we detected an additional signal, the candidate planet e, at 19.97 days with an S/N of 7.1 (Fig. 1). We fit this planet candidate with juliet[65], [66] using batman[67] for the transit model and the quasi-periodic kernel GP included in celerite[68] to model the stellar variability. We modeled the K2 light curve without the transits of K2-19 b, c, and d. We used a normal prior for the stellar density and uninformative priors for the rest of the parameters. We obtained an orbital period of \(19.9657_{-0.0067}^{+0.0047}\) days, a transit duration of \(4.28 \pm 0.31\) hours, and a planet size of \(1.26 \pm 0.17\) \(\mathrm{R_E}\). We repeated the analysis after removing the transits of the candidate planet e from the model. A model comparison strongly favors [69] the model with planet e over the model without the planet (log-Bayes factor \(\mathcal{Z}_{\rm planet~e}/\mathcal{Z}_{\rm no~planet~e}\) of \(5.60 \pm 0.37\)). We repeated the analysis with juliet with the timing of each transit as a free parameter. The individual timings have a precision of about 20 minutes and are compatible with a linear ephemeris. Interestingly, the candidate planet e is close to the 5:3 MMR with K2-19 c. We included this planet candidate (K2-19 e) in the analysis of the system.

Subsequent searches yielded low S/N candidates that are not favored by the model comparison. One interesting candidate that emerged is a \(0.591 \pm 0.025\) \(\mathrm{R_E}\) planet with a period of 1.89 days (close to the 4:3 MMR with K2-19 d), however. When this candidate was added to the analysis described in Sect. 4.2, its mass was \(0.089 \pm 0.035\) \(\mathrm{M_E}\).

4.2 Photodynamical modeling↩︎

We fit the observed photometry accounting for the gravitational interactions between the five bodies assumed in the system using a photodynamical model. Its positions and velocities in time were obtained through an n-body integration. The sky-projected positions were used to compute the light curve [70] with batman[67] using a quadratic limb-darkening law [71], which we parameterized following [72]. To account for the integration time of 30 minutes of the K2 light curve, the model was oversampled by a factor of 30, and then binned back to match the cadence of the data points [73]. We used the n-body code REBOUND[74] with the WHFast integrator [75] and an integration step of 0.01 days, which resulted in a maximum error [76] of \(\sim\)​5 ppm for the photometric model. The light-time effect [77] was not included; with an expected amplitude of \(\sim\)​6 ms in the transit timing, it is a negligible effect for this system. We included some transit timings from the literature, assuming that the reported time is the calculated time of minimum sky-projected planet–star separation. Table 1 specifies whether we used the time-series photometry or only the reported timing for each observation. The model was parameterized using the stellar mass and radius, limb-darkening coefficients, planet-to-star mass and radius ratios, and Jacobi orbital elements (Table 4) at the reference time, \(t_{\mathrm{ref}} = 2456829.222\) BJD\(_{\mathrm{TDB}}\), during the K2 observations. Because the problem is symmetrical, we fixed the longitude of the ascending node of K2-19 b \(\Omega_{\mathrm{b}}\) at \(t_{\mathrm{ref}}\) to 180, and we limited the inclination of K2-19 c to \(i_c<90\). In order to reduce the correlation between the jump parameters, we used the [78] parameterization for the eccentricity and the argument of pericenter of K2-19 b and c.

We used a GP regression [68] for the model of the error terms of the transit light curves. We used the rotational kernel for K2 and an approximate Matern kernel for the rest of the datasets. For CHEOPS observations, we used the satellite roll angle as a regressor instead of the time for the rest, and we added a linear model with \(\log\)-background (for the transit on May 10, 2024, we also added a quadratic term). For the ECAM transits, we added linear models with the point spread function centroid shifts, full width at half maximum, and peak flux value, and a quadratic term for the airmass.

In total, the model had 149 free parameters. We used normal priors for the stellar mass and radius from Sect. 3 and uninformative prior distributions for the rest of the parameters. The joint posterior distribution was sampled using the emcee algorithm [79], [80].

Our parameterization of the eccentricity and argument of pericenter for planets b and c leads to a nonuniform prior in eccentricity. To account for this, we applied the algorithm called sampling importance resampling [81] as an a posteriori correction.

5 Results↩︎

In Table 4 we list the median and the 68% credible interval (CI) of the marginal distribution of the inferred system parameters. The one- and two-dimensional projections of the posterior sample of selected system parameters are shown in Fig. 12. The MAP model of the transit photometry is plotted in Fig. 2. The noise-model-corrected transits for each planet are plotted in Figs. 13 to 16. Figure 17 shows the posterior of the planet orbits. The value of \(P_{\rm rot}\) derived from the K2 data with the celerite rotational kernel is \(18.36 \pm 0.51\) days, which is lower (\(3.7\,\sigma\)) than the \(20.54 \pm 0.30\) days of [18] that were obtained by applying the autocorrelation function [82] to the K2 light curve.

The current transit-timing \(1\,\sigma\) uncertainty for K2-19 d and e is about 7 and 33 hours, respectively. With a transit depth of about 190 and 270 ppm for K2-19 d and e, respectively, these planets are difficult to follow-up (or confirm). Notably, their shallow transit depths place them below the detection threshold of TESS, which otherwise would have been well suited for recovering planets with currently poorly constrained ephemerides. In addition, K2-19 e might not always transit in the long term (Sect. 5.2). If the candidate K2-19 e is not real, the data provide an upper limit of 2.0 \(\mathrm{M_E}\) at a confidence of 99% for a body in 5:3 MMR with K2-19 c. We verified the results for planets b and c did not change significantly when planet e was excluded from the model. In the subsequent sections, we concentrate on K2-19 b and c because we have achieved a more detailed characterization of these planets.

5.1 Evolution of the orbital parameters during the observations↩︎

We studied the changes in the orbital parameters of the planets during the period covered by the observations (10 years) from 1000 random draws from the posterior distribution. The eccentricities of K2-19 b and c vary significantly during the observations for most posterior samples and are correlated. The periodicity is generally dominated by the superperiod (Fig. 3).

The correlation between the eccentricities and difference in the longitudes of periastron exhibits three types of behavior (see Fig. 4): Libration around \(\varpi_{\mathrm{b}}-\varpi_{\mathrm{c}} = 0°\) (53% of the posterior samples), circulation (41%), and libration around \(\varpi_{\mathrm{b}}-\varpi_{\mathrm{c}} = 180°\) (6%). This behavior is discussed in detail in the companion paper (hereafter C2).

The period ratio librates around its mean value \(3/2+\Delta\sigma\) with a small amplitude \(A_\sigma\Delta\sigma\) and with a period equal to the resonant period (see Sect. 5.3), each determined by conditions at the time of formation, which generally result in two-planet systems having period ratios slightly larger than exact commensurability [4]. The posterior median and 68.3% CI values for \(\Delta\sigma\) and \(A_\sigma\) are \(0.00291^{+0.00001}_{-0.00002}\) and \(0.40\pm 0.02\), respectively, suggesting that the system is resonant and close to its equilibrium state (Sect. 5.3).

5.2 Secular behavior↩︎

Over longer timescales, the precession of the orbital planes exhibits a period of approximately 165 years (Fig. 20), which results in oscillations in both the difference in orbital inclinations and the relative longitudes of the ascending nodes. During certain epochs, this configuration permits mutual eclipses between planets b and c.

5.3 Transit-timing variations: Mass estimates↩︎

Figure 5

shows the posterior TTVs of K2-19 b and c, obtained from the times of minimum projected separation of the star and planet. Two well-separated periodicities are clearly evident, with the ratio of the associated super and (shorter, noncommensurate) resonant periods, \(P_S\) and \(P_R\), respectively, which immediately indicates that the system is formally resonant and in the librating state, even before any analysis is performed. Moreover, the magnitudes of the amplitudes of each component are similar, which reveals (again without a detailed analysis) that the system is close to (but not at) the equilibrium state of the resonance (cyan star in Fig. 8): This is the primary fixed point of the phase space in which the system resides. An algorithm for drawing these conclusions by simple inspection is given in C2.

When a system is in this state and the quality of the timing data is high enough, the Bayesian analysis process is expected to find little correlation between the masses and eccentricities, in contrast to the mass-eccentricity degeneracy that is understood to plague many systems exhibiting TTVs. Figure 6

shows selected posterior projections (see also Figure 12), demonstrating that this is indeed the case for the K2-19 system, while the planet masses themselves are correlated, as are the two eccentricities (panels (e) and (f), respectively). A full analysis of the nature of these correlations is presented in C2 [83], but we note that the masses can be determined independently when the resonant and superfrequencies are distinct and their associated Fourier contributions to the TTV signal can be accurately separated and modeled. Moreover, while in general the interference between the resonant and super harmonics is nonlinear, for a system close to the primary fixed point such as K2-19 b,c there is little interference, and the masses can be directly determined from the amplitudes of the supercomponent together with the ratio of the superperiod and resonant period. This result is summarized below, and a more general treatment for systems that are not necessarily close to the primary fixed point (using the information in the resonant harmonic) is presented in C2.

The analysis in C2 is based on a new formalism for the study of the dynamics of systems near first-order resonances (Mardling, in prep.). A key result we used here is that for a system that is close to the primary fixed point of the forced part of the associated phase space,⁶ a quantity \(R\) called the resonance parameter relates the ratio of the superperiod to the resonant period to the planet-to-star mass ratios such that \[\begin{align} R &=& (P_S/P_R)^2-1\nonumber\\ &=& 3\left(\frac{\sigma}{\Delta\sigma}\right)^2\left(\overline{m}_b+ \alpha\, \overline{m}_c\right) \left(\mathrm{v}_1(\alpha)\,e_b^{(+)}+\mathrm{v}_2(\alpha)\,e_c^{(+)}\right), \label{PR} \end{align}\tag{1}\] for a system close to the \(n+1:n\) commensurability. Here, \(\sigma=(n+1)/n\), \(\alpha=\sigma^{-2/3}\) is the ratio of the semimajor axes at exact commensurability, \(\Delta\sigma\) is the distance of the primary fixed point to exact commensurability,⁷ \(\overline{m}_{b,c}=m_{b,c}/m_*\) are the planet-to-star mass ratios, \(\mathrm{v}_1(\alpha)\) and \(\mathrm{v}_2(\alpha)\) are constants related to Laplace coefficients, such that \(\mathrm{v}_1=-2.025\) and \(\mathrm{v}_2=2.484\) for \(\sigma=3/2\), and \[e_b^{(+)}=\frac{\alpha^{-1/2} |\mathrm{v}_1|\,\overline{m}_c}{n\Delta\sigma} {\rm and} e_c^{(+)}=\frac{\mathrm{v}_2\,\overline{m}_b}{n\Delta\sigma} \label{ebc}\tag{2}\] are the values of the forced eccentricities at the primary fixed point. Equation (1 ) represents one constraint on the three quantities \(\overline{m}_b\), \(\overline{m}_c\), and \(\Delta\sigma\). The other two come from the amplitudes of the supercomponents of the TTV signals, as described below.

As discussed in C2, an accurate model for the TTVs of a system in the librating state is such that the \(jth\) TTV of planet b is given by \[\begin{align} {\rm TTV}_b(j) &=& A_S^{(b)}\cos[\nu_S (t_j-t_{\rm ref})+\beta_S]\nonumber\\ && +A_R^{(b)}\cos[\nu_R (t_j-t_{\rm ref})+\beta_R]+c_1^{(b)}t_j+c_0^{(b)} \label{TTVb} \end{align}\tag{3}\] with a similar expression for the \(kth\) TTV of planet c. Here, \(t_j\) is the \(jth\) midtransit time, \(t_{\rm ref}\) is the reference time, \(\bar P_b\) is the long-term average value of the inner orbital period, which is well approximated by its linear ephemeris value for a long enough baseline, \(A_S^{(b)}\) and \(A_R^{(b)}\) are the amplitudes of the supercomponent and resonant component, \(\beta_S\) and \(\beta_R\) are the associated phases, and \(c_0^{(b)}\) and \(c_1^{(b)}\) are constants that tend to zero in the case of an unlimited observing baseline (or more generally, model any long-term secular trends that are not accounted for in the model). For our purposes, we have in particular that \[A_S^{(b)}\simeq\bar P_b\,e_b^{(+)}/\pi\propto\overline{m}_c {\rm and} A_S^{(c)}\simeq\bar P_c\,e_c^{(+)}/\pi\propto\overline{m}_b, \label{ASbc}\tag{4}\] where \(e_b^{(+)}\) and \(e_c^{(+)}\) are given in (2 ), while \(A_R^{(b)}\propto\overline{m}_c\) and \(A_R^{(c)}\propto\overline{m}_b\) and the phases \(\beta_S\) and \(\beta_R\) depend on the forced eccentricities away from their fixed-point values (C2). If the K2-19 b,c system were precisely at the fixed point, we would have \(A_R^{(b)}=A_R^{(c)}=0\) and \[\beta_S=n\lambda_b(t_{\rm ref})-(n+1)\lambda_b(t_{\rm ref}),\] the latter being minus the longitude of conjunctions at the reference time, while if the system were farther from the fixed point, the resonant contribution would dominate the TTV signal. Figure 7

plots the best fits to the median TTV values shown in Fig. 5, and the individual contributing harmonics and residual linear trends, verifying that the superamplitude and resonant amplitude are similar, while \(P_S/P_R\simeq 3.8\). The values for the Fourier parameters are listed in Table [table:fourier],

      \tiny
\centering
\caption{Fourier and photodynamical mass estimates}\label{table:fourier}
\begin{tabular}{lccccccc}
\hline
Parameter  & Units & Median and 68.3\% CI  & MAP & \phantom{+}$\nu_R/\nu_S$, $A_S$ $^{(1)}$ 
&  \phantom{-}Petigura$^{(2)}$ & From Radii$^{(3)}$ \\
&& Photodyn & Photodyn & Fourier & Photodyn\phantom{-}  \\
&& 10 yr & 10 yr & 10 yr    & 3.3 yr \phantom{-} &  \\
\hline
\vspace{-2mm}
&\\
$m_b$  & $M_\oplus$ &  ${\bf 30.8\pm 1.3}$  & \bf{29.62}  & {\bf 30.03}  &  { $\bf 32.4\pm 1.7$}  & 39.6 \\
$m_c$  & $M_\oplus$ &  {$\bf 11.12\pm 0.44$}  & {\bf 10.94} & {\bf 11.08}   &  { $\bf 10.8\pm 0.6$} & 18.0\\
\vspace{1mm}
$\Delta m_{b,c}/m_{b,c}$ &  &    & $-0.04$, $-0.02$  & $-0.02$, $-0.004$      & $+0.05$, $+0.03$   & $+0.29$, $+0.62$           \\
$\rho_b$ & ${\rm g\,cm}^{-3}$ & 0.59 & 0.57 & 0.58  &  0.52 & 0.76  \\
\vspace{1mm}
$\rho_c$ & ${\rm g\,cm}^{-3}$  & 1.05 & 1.03 & 1.04  &  0.86 & 1.69  \\
$P_S$  ($\overline P_S$) & d &  &  & 2711.2    & & (2751.9) \\
\vspace{1mm}
$P_R$ & d &  & & 717.4     \\
\vspace{1mm}
$A_S^{(b)}$, $A_R^{(b)}$ &  min &  & & 53.4, (41.7)    \\
\vspace{1mm}
$A_S^{(c)}$, $A_R^{(c)}$  &  min &  & & 233.2, (148.3)   \\
\vspace{0.5mm}
$\beta_S$, \, $\beta_R$ & deg & $[78.1]^{+1.5}_{-1.4}$ \hspace{1.0mm} $[137.3]^{+2.2}_{-1.9}$ & [76.1], [138.1]& 
(76.5, 144.3)  \\
\vspace{0.5mm}
$\lambda_n$ & deg & $85.8^{+1.7}_{-1.9}$  &  82.8 &    & \phantom{-}37.5$^{(4)}$ \\
$R$ & & $[14.81]^{+0.37}_{-0.26}$ & [15.49] & 13.28   & [15.7]\\
\vspace{1mm}
$\Delta\sigma$ & & $[0.0029]^{+1.e-5}_{-2.e-5}$ &  [0.0029]  &  0.0030    &  [0.0029]\\
\vspace{1mm}
$e_b^{(+)}$ ($e_b$)    &  & $[0.015]^{+0.0002}_{-0.0002}$ ($0.043^{+0.024}_{-0.024}$)  & [0.015] (0.069) & 0.015   &  [0.015]  ($0.20^{+0.03}_{-0.03}$) &   \\
\vspace{1mm}
$e_c^{(+)}$  ($e_c$)  &  & $[0.044]^{+0.0005}_{-0.0004}$ ($0.067^{+0.017}_{-0.017}$) & [0.045] (0.078)& 0.043    &  [0.047]  ($0.21^{+0.03}_{-0.03}$)&  \\
\hline
\end{tabular}
\tablefoot{\tiny 
\phantom{.}\\
(1): $\nu_R/\nu_S$, $A_S$ / Fourier A uses the frequency ratio and superamplitudes only to
estimate masses;\\
(2): Petigura / Photodyn uses system parameters from the photodynamical analysis of \cite{Petigura2020};\\
(3): Estimates using mass-radius relation $m_p/M_E=(R_p/0.56\,R_E)^{1/0.67}$ \cite{Mueller2024};\\
(4): \cite{Petigura2020} value minus $180^{\rm o}$ to account for the fact that  the authors assumed $\Omega_b(t_{\rm ref})=0$;\\
$[${\it Square brackets}$]$ show the quantity calculated analytically from elements and masses using the formalism of Mardling (2025a);\\
$\Delta m_{b,c}/m_{b,c}=$ relative difference in masses compared to median value;\\
$\overline P_S$ (defined from ephemeris periods such that $\bar P_S^{-1}= n\bar P_b^{-1}-(n+1)\bar P_c^{-1}$) listed in the last column for comparison with the Fourier superperiod;\\
{\textit{\textbf {Boldface entries}}}: highlighted for comparison;\\
Amplitudes and phases in brackets were not used directly in mass estimates.  \\
Last two rows: Entries in brackets are the median and 68.3\% CI eccentricities for comparison with the fixed-point values $e_b^{(+)}$ and  $e_b^{(+)}$.\\
The reference mean densities are $\rho_{\rm Saturn}=0.69\,{\rm g\,cm}^{-3}$, $\rho_{\rm Uranus}=1.15\,{\rm g\,cm}^{-3}$.
}

together with \(R\), \(\Delta\sigma\), \(e_b^{(+)}\) and \(e_c^{(+)}\). The planet mass estimates listed in the column labeled “Fourier” were obtained by solving (1 ), (2 ) and (4 ) for \(\overline{m}_b\), \(\overline{m}_c\) and \(\Delta\sigma\), yielding values that are 2.5% and 0.4% lower than the median values obtained from the photodynamical analysis for planets b and c, respectively, and are well inside our error bars. We also list our median and MAP values for the planet masses (taken from Table 4) and those of [13].

The accuracy of the Fourier mass estimates can be attributed to a good estimate of \(\Delta\sigma\) via the resonance parameter \(R\) (which is vital because the superamplitudes go like \(m_{b,c}/\Delta\sigma\)); if we had used \(\Delta S=\bar P_c/\bar P_b-(n+1)/n\) as an estimate for \(\Delta\sigma\) in (4 ), the planet masses would both have been underestimated by about 30%.

5.4 Transit-timing variations: Eccentricity constraints↩︎

The high eccentricities of K2-19 b and c reported by [13] (both about 0.2) have inspired much discussion about the formation of the system, including the influence of additional planets [16], [84]–[86] and stochastic stirring by turbulent eddies [86]. Based on seven additional years of transit data, however, our posterior is significantly better constrained, with median and 68.3% CI values at \(e_b=0.043\pm0.024\) and \(e_c=0.067\pm 0.017\). While substantially lower than 0.2, these values still imply some free eccentricity, with values lower than 0.005 unlikely at a confidence of 99%. Panel (f) of Figure 6 again shows the posterior coverage of the \(e_b-e_c\) plane (see the figure caption for the meaning of the colors and the discussion below). The shape of this posterior projection is a result of a third integral (in addition to the total energy and angular momentum), which exists at first order in eccentricity for systems that are close to a first-order commensurability. Because of this integral, the free components of the Runge-Lenz vectors for each planet orbit are approximately constant,⁸ and as a result, their contribution to the TTVs is limited because their variation in time is second order in eccentricity. This leads to the degeneracy that is visible in the figure and is confined to the upper half of the allowed region (the latter being the entire colored part) as a result of second-order effects and of the sensitivity to the average transit duration time (see the next section).

The colors in Figure 6 (f) are associated with the libration and/or circulation behavior of the standard resonance angles \(\phi_1=n\lambda_b-(n+1)\lambda_c+\varpi_b\) and \(\phi_2=n\lambda_b-(n+1)\lambda_c+\varpi_c\) and the difference in the longitudes of periastron \(\varpi_b-\varpi_c\); this can be derived from the new definitions of forced and free eccentricity given in C2 together with the complex functions⁹ \[u=-3\left(\frac{\sigma}{\Delta\sigma}\right)^2\left(\overline{m}_b+ \alpha\, \overline{m}_c\right) \left(\mathrm{v}_1\,z_b+\mathrm{v}_2\,z_c\right) \label{u}\tag{5}\] and \[\upsilon=-3\left(\frac{\sigma}{\Delta\sigma}\right)^2\left(\overline{m}_b+ \alpha\, \overline{m}_c\right) \left(g^{-1}\mathrm{v}_2\,z_b-g\mathrm{v}_1\,z_c\right), \label{v}\tag{6}\] posterior values that are plotted in Figure 8.

Here, \(g=(m_c/m_b)^{1/2}\alpha^{-1/4}\), \(z_1=e_b\,{\rm e}^{i\phi_1}\) and \(z_2=e_c\,{\rm e}^{i\phi_2}\). We note from (1 ) that \(u=R\) at the primary fixed point (at which \(\phi_1=0\) and \(\phi_2=\pi\)). We also note that the true resonance angle \(\psi\), defined such that \(u=|u|{\rm e}^{i\psi}\), librates with a small amplitude for the whole posterior. Finally, we note that the third integral mentioned above is simply \(\upsilon\upsilon^*={\rm constant}\) (so that the radius of the black curve in the right panel of Figure 8 is constant). The details are provided in C2.

5.5 The transit duration: Further eccentricity constraints↩︎

The transit duration was computed as the difference between the first and fourth contact using the sky-projected planet-star separation. The transit durations of K2-19 b and c are shown in Fig. 9

and exhibit transit duration variations (TDVs). The transit duration of K2-19 b varies by about 6 minutes peak to peak over the 10-year baseline, and that of K2-19 c varies by about 20 minutes. This is 3% and 11%, respectively, of the median transit durations. As we show in C2 the variation period is the superperiod, and the amplitudes are proportional to the primary fixed-point eccentricities \(e_b^{(+)}\) and \(e_c^{(+)}\) and not the true eccentricities, as is commonly assumed. The average transit durations of planets b and c do provide constraints via the quantities \(e_b\sin\omega_b\) and \(e_c\sin\omega_c\), however, which further confines the \(e_b-e_c\) posterior projection in the upper half of the colored region in panel (f) of Figure 6.

5.6 Comparison with previous works↩︎

The first photodynamical analysis of K2-19 with a time span of 285 days was presented by [11]. They found \(e_{\rm b} = 0.119^{+0.082}_{-0.035}\) and \(e_{\rm c} = 0.095^{+0.073}_{-0.035}\) at 2456813 BJD\(_{\rm TDB}\). These values are different by \(1.8\,\sigma\) (\(e_{\rm b}\)) and \(0.7\,\sigma\) (\(e_{\rm c}\)) from our results. Their mass ratios \(m_b/m_\star = (13.4^{+4.3}_{-2.9})\times10^{-5}\), \(m_c/m_\star = (5.5^{+2.0}_{-1.3})\times10^{-5}\), and \(m_c/m_b = 0.42 \pm 0.12\) agree with our results with differences of \(1.0\,\sigma\), \(1.3\,\sigma\), and \(0.5\,\sigma\), respectively.

The next photodynamical analysis of K2-19 was performed by [13] with a time span of 3.3 years. They obtained osculating \(e_{\rm b} = 0.20 \pm 0.03\), \(e_{\rm c} = 0.21 \pm 0.03\), and well-aligned apsides (\(\omega_{\rm c} - \omega_{\rm b} = 2 \pm 2°\)). We did not find the reference time for the orbital parameters in [13]. Assuming that their reference time is close to ours, the difference from our result is \(4.1\,\sigma\) for both planets, which cannot be explained by orbital evolution if the reference time was significantly different from ours. Despite the differences in eccentricities, their masses of K2-19 b and c agree with our results (for the reasons explained in Sect. 5.3). Based on the difference between the transit times predicted by [13] and those observed (Fig. 25), which are up to \(6.7\,\sigma\) and 14 hours, we suspect an issue in their analysis. As shown in Sect. 9, K2-19 d and e cannot have caused this discrepancy, and neither can the shorter time span of the observations. Perhaps the discrepant individual timing they found for the LCOGT partial transit (\(3.5\,\sigma\) off based on our results) has driven their solution. As shown in Fig. 6 of [13], their model already struggled to fit this transit timing. When [16] enforced a low-eccentricity solution, the model for the LCOGT timing fell at the same value as we found, which corresponds to \(12.4 \pm 0.9\) minutes in the residuals of Fig. 7 in [16].

We briefly compare our results with the analyses that only used radial velocities from [17] and [18]. [17] found radial velocity amplitudes of \(K_{\rm b} = 9.6^{+1.8}_{-1.6}\) m s\(^{-1}\)and \(K_{\rm c} = 7.5 \pm 2.1\) m s\(^{-1}\), which differ by \(0.6\,\sigma\) and \(2.0\,\sigma\), respectively, from our results. On the other hand, [18] found \(K_{\rm b} = 18.8 \pm 2.4\) m s\(^{-1}\), which differs by \(3.4\,\sigma\) from our results, and did not significantly detect K2-19 c.

5.7 Internal structure↩︎

Internal structure retrieval is especially degenerate for gas giants. [13] used the models of [87] and found that K2-19 b has a hydrogen-helium envelope mass fraction of \(f_{\rm env} = 44 \pm 3\)% by mass, which challenges the core-accretion theory. Runaway accretion begins at about \(f_{\rm env} \approx 50\)%, and fine-tuning would be required for the disk to disappear at just the moment when K2-19 b approached the runaway phase.

We followed [88] and investigated the internal structure of K2-19 b using GASTLI[89], [90], which is a one-dimensional coupled interior-atmosphere model for warm (\(T_{\rm eq} < 1000\) K) gas giants (between 17 \(\mathrm{M_E}\) and 6 \(\mathrm{M}_{\rm J}\)). The interior model is composed of two layers, that is, the core and the envelope. The core is assumed to be a 1:1 mixture of silicate rock and water, whereas the envelope is a mixture of H/He and water. The model is parameterized with two values: the core mass fraction (CMF) and the mass fraction of water (representative of metals) in the envelope (\(Z_{\rm env}\)).

We generated a grid of GASTLI models for \(T_{\rm eq}=857\) K, C/O=0.55 (solar value), masses from 24 to 38 \(\mathrm{M_E}\) with increments of 0.25 \(\mathrm{M_E}\), a CMF from 0 to 0.99 with increments of 0.1, log(Fe/H) from -2 to 2.4 with increments of 0.5, and \(T_{\rm int}\) from 50 to 250 K with increments of 25 K. For the retrieval, we used the planetary mass and radius and the stellar age as observable parameters. We sampled from the posterior distribution with emcee[79], [80]. The interior parameter priors and posteriors are listed in Table 5. The one- and two-dimensional projections of the posterior sample [91] of selected system parameters are shown in Fig. 21. The envelope mass fraction [87] is \(0.283^{+0.38}_{-0.049}\). The total water (metals) mass fraction is \(Z_{\rm planet} = 0.729^{+0.041}_{-0.14}\). We compared our results with those by [92]. We computed the stellar metal mass fraction as \(Z_{\star} = 0.014 \times 10^{\rm [Fe/H]_\star}\) to infer the heavy-element enrichment relative to the host star: \(Z_{\rm planet}/Z_{\star} = 48.1^{+4.1}_{-8.8}\). This value is compatible with formation via core accretion for a planet mass of \(0.0968 \pm 0.0039\) \(\mathrm{M}_{\rm J}\)[92].

3pt

6 Discussion↩︎

Figure 10 illustrates that K2-19 b and c lie within the Neptune savanna [97], which is a moderately populated region in the period–radius diagram that is found at orbital periods longer than those that define the Neptune desert [98], [99] and the Neptunian ridge [96]. Both planets exhibit low densities (\(0.590 \pm 0.053~\mathrm{g\;cm^{-3}}\) and \(1.046 \pm 0.097~\mathrm{g\;cm^{-3}}\) for K2-19 b and K2-19 c, respectively), as anticipated for the planets in the savanna [88]. K2-19 c satisfies the selection criteria of the sample defined by [100], which supports their finding that resonant sub-Neptunes tend to be puffier.

The 3:2 mean-motion resonance is one of the most common configurations observed in known multiplanetary systems [3]. This resonance is thought to be particularly robust, to often act as the first capture point for migrating planets, and to provide long-term dynamical stability. The updated orbital parameters of the K2-19 system are broadly consistent with a disk-driven migration scenario. The nonzero free eccentricity suggests, however, that an additional excitation mechanism must have operated during or after migration. The system is part of the ATREIDES spin-orbit angle survey [101], which will bring further constraints on the 3D orbital architecture of the planets and their dynamical origin.

The precision on the planet-to-star mass ratio is 0.8% for planet b and 1.5% for planet c, and the masses have a relative precision of 4.2% for planet b and 4.0% for planet c, which are dominated by the 3.8% uncertainty on the stellar mass. The two well-separated periodicities in the TTVs (super and resonant) enables a robust determination of the planetary masses independent of the eccentricities. In contrast, constraints on the free eccentricities are weak because an inherent degeneracy is associated with low-eccentricity two-planet systems, and the libration and/or circulation behavior of the usual resonance angles and difference in the longitudes of periastron is similarly degenerate. Thus, while future transit observations aimed at refining the system architecture should prioritize epochs with the largest predicted timing uncertainties, these observations are unlikely to constrain the free eccentricities much more than they are at present.

Although faint (J=11.6 mag), K2-19 b has a transmission spectroscopy metric [102] of \(64.8 \pm 4.5\), which is among the ten most favorable planets of this size (6-7 \(\mathrm{R_E}\)) with an equilibrium temperature below 1000 K. The James Webb Space Telescope [103] might be able to measure the metallicity of the atmosphere of K2-19 b and place an additional constraint on the planet interior.

[13] used photodynamical modeling for K2 data, but only transit timing for the other observations. As we showed in Sect. 9, the likely cause of their high retrieved eccentricities is a single timing from a partial transit near the end of the transit-timing dataset that was miscalculated (\(3.5\,\sigma\), or \(12.4\pm3.6\) minutes, off). This highlights the susceptibility of the results to single outliers when the available data are barely enough to constrain the solution. In this case, using a photodynamical model for all observations might have enabled a better error propagation and avoided an overestimation of the eccentricities in the past.

Data availability↩︎

Light curves, samples from the posterior, and transit times forecast up to the year 2035 can be found at https://zenodo.org/records/17084030.

7 Additional figures and tables↩︎

2pt

Table 6: Modeling of the SED.

Parameter		Prior	Posterior median
			and 68.3% CI
Effective temperature, \(T_{\mathrm{eff}}\)	(K)	\(N\)(5350, 110)	\(5363 \pm 74\)
Surface gravity, \(\log g\)	(cgs)	\(N\)(4.42, 0.06)	\(4.527 \pm 0.029\)
Metallicity, \([\rm{Fe/H}]_\star\)	(dex)	\(N\)(0.03, 0.02)	\(0.026 \pm 0.019\)
Distance	(pc)	\(N\)(296.9, 1.7)	\(297.2 \pm 1.7\)
\(E_{\mathrm{(B-V)}}\)	(mag)	\(U\)(0, 3)	\(0.030^{+0.052}_{-0.023}\)
Jitter Gaia	(mag)	\(U\)(0, 1)	\(0.169^{+0.24}_{-0.092}\)
Jitter 2MASS	(mag)	\(U\)(0, 1)	\(0.057^{+0.090}_{-0.037}\)
Jitter WISE	(mag)	\(U\)(0, 1)	\(0.042^{+0.10}_{-0.032}\)
Mass, \(m_\star\)	(M\(_\odot\))		\(0.885 \pm 0.034\)
Radius, \(R_\star\)	(R\(_\odot\))		\(0.845 \pm 0.017\)
Density, \(\rho_\star\)	(\(\mathrm{g\;cm^{-3}}\))		\(2.06 \pm 0.16\)
Isochronal age	(Gyr)		\(4.8 \pm 4.3\)
Luminosity	(L\(_\odot\))		\(0.531 \pm 0.030\)

Table 7: Parameters of K2-19 from the literature.

Work	\(T_{\rm eff}\) (K)	\(\log g\)	[Fe/H]\(_\star\)	[M/H]\(_\star\)	\(m_{\star}\) (\(M_{\odot}\))	\(R_{\star}\) (\(R_{\odot}\))	d (pc)
[6]	\(5519^{+49}_{-82}\)		\(-0.27 \pm 0.10\)		\(0.84 \pm 0.04\)	\(0.81^{+0.09}_{-0.05}\)	\(291^{+33}_{-20}\)
[7]	\(5230 \pm 420\)	\(4.39 \pm 0.79\)	\(0.38 \pm 0.23\)		\(0.92 \pm 0.14\)	\(1.03 \pm 0.20\)
[11]	\(5390 \pm 180\)	\(4.42 \pm 0.34\)	\(0.19 \pm 0.12\)		\(0.918^{+0.086}_{-0.070}\)	\(0.926^{+0.19}_{-0.069}\)
[12]	\(5345 \pm 17\)	\(4.394 \pm 0.050\)	\(0.07 \pm 0.03\)		\(0.902 \pm 0.011\)	\(0.914 \pm 0.027\)
[61]	\(5340 \pm 160\)	\(4.520 \pm 0.040\)	\(0.019^{+0.12}_{-0.080}\)		\(0.901^{+0.055}_{-0.031}\)	\(0.857^{+0.069}_{-0.054}\)	\(296^{+32}_{-29}\)
[8]	\(5391 \pm 50\)	\(4.60 \pm 0.10\)		\(-0.03 \pm 0.08\)
[108]	\(5430 \pm 60\)	\(4.63 \pm 0.07\)	\(0.10 \pm 0.04\)		\(0.93 \pm 0.05\)	\(0.86 \pm 0.04\)
[109]	\(5322 \pm 25\)	\(4.510 \pm 0.028\)		\(0.05 \pm 0.01\)
[18]	\(5250 \pm 70\)	\(4.50 \pm 0.10\)	\(0.10 \pm 0.05\)		\(0.918 \pm 0.064\)	\(0.88 \pm 0.11\)
[110]	\(5355 \pm 35\)	\(4.46 \pm 0.06\)	\(0.05 \pm 0.02\)		\(0.86 \pm 0.06\)
[111]	\(5320 \pm 100\)	\(4.553 \pm 0.077\)		\(0.05 \pm 0.10\)	\(0.92 \pm 0.11\)	\(0.840 \pm 0.043\)	\(289.8 \pm 5.6\)
[112]	\(5418 \pm 21\)	\(4.54 \pm 0.04\)	\(0.03 \pm 0.02\)
[113]	\(5481 \pm 31\)	\(4.726 \pm 0.051\)	\(0.216 \pm 0.029\)		\(1.22^{+0.17}_{-0.15}\)	\(0.793^{+0.021}_{-0.020}\)	\(289.8^{+5.7}_{-5.5}\)
[13]	\(5320 \pm 100\)	\(4.51 \pm 0.08\)	\(0.06 \pm 0.05\)		\(0.88 \pm 0.03\)	\(0.82 \pm 0.03\)
[53]	\(5346 \pm 27\)	\(4.42 \pm 0.06\)	\(0.03 \pm 0.02\)		\(0.846 \pm 0.006\)	\(0.864 \pm 0.028\)	299.3
[114]	\(5285 \pm 82\)	\(4.53 \pm 0.15\)	\(0.09 \pm 0.09\)		\(0.89^{+0.07}_{-0.11}\)
[115]	\(5378 \pm 82\)	\(4.712 \pm 0.049\)	\(0.146 \pm 0.029\)
[116]	\(5360 \pm 110\)	\(4.50 \pm 0.16\)	\(0.055 \pm 0.054\)
[117]	\(5348^{+16}_{-20}\)	\(4.509^{+0.013}_{-0.015}\)		\(0.026^{+0.015}_{-0.018}\)			\(295.7^{+1.7}_{-1.6}\)
[118]	\(5247 \pm 38\)	\(4.511 \pm 0.027\)	\(0.032 \pm 0.055\)
[119]	\(5420 \pm 100\)	\(4.55^{+0.01}_{-0.02}\)	\(0.14 \pm 0.06\)		\(0.91 \pm 0.02\)	\(0.84 \pm 0.01\)
This work (prior)		\(4.531 \pm 0.044\)			\(0.885 \pm 0.034\)	\(0.845 \pm 0.039\)	\(296.9 \pm 1.7\)
This work (posterior)		\(4.556 \pm 0.024\)			\(0.892 \pm 0.034\)	\(0.824 \pm 0.024\)

8 Determination of transit times↩︎

Transit times of multiple transits observations can be obtained by fitting them simultaneously while assuming constant transit parameters. This approach can obtain more precise values than when each transit is analyzed individually, be less sensitive to stellar activity or instrumental systematics, and it is the safest option if there are partial transits among the observations. Depending on the photometric precision, the transit timing of partial transits presents a degeneracy with the transit duration. This degeneracy is broken if the the signal of limb darkening is detected. However, when the transits present TDVs, assuming one transit shape for all transits can bias the derived transit timings. To obtain the transit timings in K2-19, we found a compromise by fitting together groups of transits observed close in time. But we also fit individual transit times for the most precise observations. Simultaneous observations of the same transit were analyzed together. For planet b, we analyzed epochs [311-362] as a group and epochs [398-438] as another group. Epochs 449, 450, 451, and 458 were analyzed individually. For planet c, epochs [228-232] and [290-291] were analyzed as groups, while epochs 272, 301, and 304 were analyzed individually. The partial transit of planet b at epoch 141 (observed with LCOGT), for which the timing presented in [13] represents an outlier, was analyzed assuming the transit parameters of K2 observations. The timings of planets d, only detected in K2 observations, were analyzed simultaneously. And the same for planet e candidate (Sect. 4.1). We modeled the datasets with juliet[65]–[68] using the same noise model as in Sect. 4.2. Table 8 presents literature and the new transit timings we computed. These timings, derived independently of the assumptions used in the photodynamical modeling (Sect. 4.2), were utilized to validate its results by comparing them with the posterior TTVs (Fig. 5).

2pt

Table 8: Transit times of the observations.

Epoch	Posterior median	Telescope	Source
	and 68.3% CI [BJD\(_{\mathrm{TDB}}\)]	Instrument
Planet d
0	\(2456811.426_{-0.016}^{+0.015}\)	K2	This work
1	\(2456813.945_{-0.023}^{+0.031}\)	K2	This work
2	\(2456816.433_{-0.014}^{+0.018}\)	K2	This work
3	\(2456818.965_{-0.019}^{+0.023}\)	K2	This work
4	\(2456821.468_{-0.021}^{+0.017}\)	K2	This work
5	\(2456823.990_{-0.014}^{+0.019}\)	K2	This work
6	\(2456826.490_{-0.018}^{+0.014}\)	K2	This work
7	\(2456828.954_{-0.011}^{+0.013}\)	K2	This work
8	\(2456831.486_{-0.015}^{+0.013}\)	K2	This work
9	\(2456834.007_{-0.011}^{+0.012}\)	K2	This work
10	\(2456836.5443_{-0.0084}^{+0.0093}\)	K2	This work
11	\(2456839.025_{-0.016}^{+0.012}\)	K2	This work
12	\(2456841.514_{-0.051}^{+0.033}\)	K2	This work
13	\(2456844.041_{-0.016}^{+0.015}\)	K2	This work
14	\(2456846.582_{-0.018}^{+0.020}\)	K2	This work
16	\(2456851.523_{-0.016}^{+0.026}\)	K2	This work
17	\(2456854.069_{-0.013}^{+0.014}\)	K2	This work
18	\(2456856.571_{-0.016}^{+0.019}\)	K2	This work
19	\(2456859.082_{-0.015}^{+0.014}\)	K2	This work
20	\(2456861.6036_{-0.0088}^{+0.019}\)	K2	This work
21	\(2456864.105_{-0.019}^{+0.012}\)	K2	This work
22	\(2456866.610_{-0.022}^{+0.024}\)	K2	This work
23	\(2456869.124_{-0.022}^{+0.042}\)	K2	This work
24	\(2456871.641_{-0.018}^{+0.014}\)	K2	This work
25	\(2456874.144_{-0.029}^{+0.018}\)	K2	This work
26	\(2456876.651_{-0.016}^{+0.038}\)	K2	This work
27	\(2456879.163_{-0.014}^{+0.028}\)	K2	This work
28	\(2456881.677_{-0.022}^{+0.018}\)	K2	This work
29	\(2456884.195_{-0.080}^{+0.032}\)	K2	This work
30	\(2456886.677_{-0.025}^{+0.022}\)	K2	This work
31	\(2456889.193_{-0.025}^{+0.028}\)	K2	This work

2pt

Table 9: No caption

Epoch	Posterior median	Telescope	Source
	and 68.3% CI [BJD\(_{\mathrm{TDB}}\)]	Instrument
Planet b
0	\(2456813.38390 \pm 0.00056\)	K2	[11]
1	\(2456821.30410 \pm 0.00088\)	K2	[11]
2	\(2456829.22194 \pm 0.00025\)	K2	[11]
3	\(2456837.13900 \pm 0.00064\)	K2	[11]
4	\(2456845.06156 \pm 0.00053\)	K2	[11]
5	\(2456852.9794 \pm 0.0028\)	K2	[11]
6	\(2456860.90005 \pm 0.00079\)	K2	[11]
7	\(2456868.82080 \pm 0.00052\)	K2	[11]
8	\(2456876.73856 \pm 0.00064\)	K2	[11]
9	\(2456884.65834 \pm 0.00080\)	K2	[11]
30	\(2457051.00413^{+0.00218}_{-0.00225}\)	FLWO 1.2 m	[12]
34	\(2457082.69550^{+0.00140}_{-0.00116}\)	TRAPPIST-South	[12]
34	\(2457082.6892 \pm 0.0019\)	NITES	[11]
35	\(2457090.6157 \pm 0.0014\)	1-m C2PU Omicron	[11]
36	\(2457098.5368 \pm 0.0022\)	Belesta 82-cm	[11]
41	\(2457138.15047^{+0.00145}_{-0.00176}\)	OAO 1.88 m / MuSCAT	[12]
133	\(2457866.8604 \pm 0.0009\)	Spitzer	[13]
141	\(2457930.2403^{+0.0015}_{-0.0013}\)	LCOGT-SAAO 1 m	This work
150	\(2458001.5368 \pm 0.0014\)	Spitzer	[13]
311	\(2459276.8341^{+0.0019}_{-0.0020}\)	ExTrA T23	This work
343	\(2459530.2882^{+0.0029}_{-0.0032}\)	TESS s45	This work
345	\(2459546.1269 \pm 0.0023\)	TESS s45	This work
346	\(2459554.0421^{+0.0021}_{-0.0023}\)	TESS s46	This work
347	\(2459561.9601 \pm 0.0023\)	TESS s46	This work
348	\(2459569.8789^{+0.0033}_{-0.0029}\)	TESS s46	This work
349	\(2459577.7967^{+0.0024}_{-0.0025}\)	TESS s46	This work
361	\(2459672.8198^{+0.0016}_{-0.0017}\)	ExTrA T23	This work
362	\(2459680.7365^{+0.0022}_{-0.0023}\)	ExTrA T23	This work
398	\(2459965.8909^{+0.0021}_{-0.0019}\)	ExTrA T23	This work
399	\(2459973.8158^{+0.0014}_{-0.0017}\)	ExTrA T123	This work
400	\(2459981.7360 \pm 0.0019\)	ExTrA T123	This work
436	\(2460266.8799^{+0.0033}_{-0.0032}\)	TESS s72	This work
438	\(2460282.7201^{+0.0021}_{-0.0022}\)	TESS s72	This work
449	\(2460369.8272 \pm 0.0012\)	CHEOPS, ExTrA T12	This work
450	\(2460377.74640^{+0.00057}_{-0.00055}\)	Euler / ECAM, ExTrA T12	This work
451	\(2460385.66556^{+0.00070}_{-0.00072}\)	Euler / ECAM, ExTrA T12	This work
458	\(2460441.10519^{+0.00064}_{-0.00066}\)	CHEOPS	This work
Planet c
0	\(2456817.2730 \pm 0.0015\)	K2	[11]
1	\(2456829.1835 \pm 0.0013\)	K2	[11]
2	\(2456841.0922 \pm 0.0014\)	K2	[11]
3	\(2456853.0091 \pm 0.0077\)	K2	[11]
4	\(2456864.9083 \pm 0.0019\)	K2	[11]
5	\(2456876.8150 \pm 0.0014\)	K2	[11]
6	\(2456888.7161 \pm 0.0013\)	K2	[11]
87	\(2457852.4774 \pm 0.0074\)	Spitzer	[13]
102	\(2458030.8645 \pm 0.0059\)	Spitzer	[13]
228	\(2459529.950^{+0.019}_{-0.040}\)	TESS s45	This work
229	\(2459541.884^{+0.015}_{-0.022}\)	TESS s45	This work
230	\(2459553.822^{+0.034}_{-0.017}\)	TESS s46	This work
232	\(2459577.602^{+0.027}_{-0.018}\)	TESS s46	This work
272	\(2460053.6499^{+0.0038}_{-0.0028}\)	ExTrA T23	This work
290	\(2460267.868^{+0.011}_{-0.013}\)	TESS s72	This work
291	\(2460279.798^{+0.017}_{-0.014}\)	TESS s72	This work
301	\(2460398.8787^{+0.0016}_{-0.0017}\)	CHEOPS	This work
304	\(2460434.5866^{+0.0023}_{-0.0022}\)	Euler / ECAM, ExTrA T23	This work
Candidate planet e
0	\(2456812.250^{+0.019}_{-0.014}\)	K2	This work
1	\(2456832.220^{+0.016}_{-0.011}\)	K2	This work
2	\(2456852.185^{+0.011}_{-0.013}\)	K2	This work
3	\(2456872.143^{+0.013}_{-0.016}\)	K2	This work

9 TTVs analysis↩︎

In this section, we analyze only the transit times of planets b and c (listed in Table 8). We exclude the partial transit observed by LCOGT, which was analyzed assuming the transit parameters of K2 observations. However, there is a significant time span between both observations. According to the photodynamical modeling, the difference between transits durations is about 6 minutes, so the determination can be biased. For epoch 34 of planet b, for which two timings can be found in the literature, we keep the one by [12] for the analysis in this section because it is more precise.

To investigate the effect of the coverage of the TTVs signal, we separate the timings into seven datasets: 2014 (only K2 data), 2015 [11], 2017 [13], 2021, 2022, 2023, and 2024. For each dataset, only transit times up to the end of the named year were considered.

The model only includes the star and planets b and c. Transit times were modeled as the moments of minimum sky-projected planet–star separation, calculated using Eq. 1 from [120], with n-body integrations performed in REBOUND, employing WHFast and an integration step of 0.01 days. We fixed the orbital inclination of both planets at 90(which biases the mutual inclination) and included a jitter term [56] for the dataset of each planet. The eccentricities and argument of pericenter were parameterized using \(\sqrt{e}\cos{\omega}\) and \(\sqrt{e}\sin{\omega}\). We used uniform priors for all parameters. The joint posterior distribution was sampled using emcee[79], [80]. The marginal posterior distribution of selected parameters are shown in Fig. 22, and the TTVs plots are in Fig. 23.

[121] asserts that the period modulation of the TTV needs to be fully mapped to avoid degeneracy in the dynamical retrieval. When comparing the 2015 dataset to [11], the photodynamical modeling appears to be less impacted by this degeneracy than the TTV analysis. The distinction between modeling just transit times and employing photodynamical modeling is more noticeable in 2015, which has less coverage of the TTV signal (see Fig. 24 and the 2015 panel in Fig. 23), than in 2024 (see Fig. 5 and the 2024 panel in Fig. 23). The posteriors from the photodynamical modeling (Sect. 4.2) are slightly more precise than those from the transit-time analysis of the 2024 dataset for the shared parameters, although both agree within \(1\,\sigma\). With the exception of the 2015 dataset (and 2024, which uses all available data), the remaining datasets predict future observed transits within \(1.6\,\sigma\) for planet b and within \(2.6\,\sigma\) for planet c (Fig. 23).

We repeated the analysis of the 2017 dataset, including the discrepant LCOGT transit timing from [13], to investigate if it could be responsible for their retrieved higher eccentricities and poor transit timing forecasts (Fig. 25). The TTV model is able to fit the discrepant timing and produces transit timing forecasts similar to those of the 2017 dataset without the LCOGT timing (see Figs. 23 and 26). However, the posterior eccentricities exhibit notable differences (see Fig. 22). The posteriors for the dataset that includes the LCOGT transit timing from [13] are narrower and centered around a higher eccentricity. Therefore, this discrepant timing could be at least partially responsible for the higher eccentricity value found by these authors. The reasoning is as follows: The [13] analysis used a photodynamical modeling for the K2 data and modeled only the transit timings for the other datasets. The discrepant LCOGT transit timing from [13] pushes to moderate eccentricity, but the photodynamical modeling, which has an additional constrain on the eccentricity from the transit shape, pushes to lower eccentricities. In the end, their posterior converges to an intermediate value between these two. If this is correct, performing photodynamical modeling of all datasets would have avoided this problem.

References↩︎