Resonance in the K2-19 system is at odds with its high reported eccentricities

K2-19 hosts a planetary system composed of two outer planets, b and c, with size of $7.0\pm 0.2~R_\oplus$ and $4.1\pm0.2~R_\oplus$ , and an inner planet, d, with a radius of $1.11\pm 0.05 R_\oplus$. A recent analysis of Transit-Timing Variations (TTVs) suggested b and c are close to but not in 3:2 mean motion resonance (MMR) because the classical resonant angles circulate. Such an architecture challenges our understanding of planet formation. Indeed, planet migration through the protoplanetary disk should lead to a capture into the MMR. Here, we show that the planets are in fact, locked into the 3:2 resonance despite circulation of the conventional resonant angles and aligned periapses. However, we show that such an orbital configuration cannot be maintained for more than a few hundred million years due to the tidal dissipation experienced by planet d. The tidal dissipation remains efficient because of a secular forcing of the innermost planet eccentricity by planets b and c. While the observations strongly rule out an orbital solution where the three planets are on close to circular orbits, it remains possible that a fourth planet is affecting the TTVs such that the four planet system is consistent with the tidal constraints.


INTRODUCTION
The numerous planet discoveries over the past decades have revealed the large diversity in sizes and orbital architecture of exoplanetary systems (Winn & Fabrycky 2015). The discovery of Hot Jupiters (Mayor & Queloz 1995) and of super-Earths have profoundly changed how we see planets and their formation outside of the Solar System narrative. Since then, planet formation theories have been adapted to take this diversity into account. The models now prefer a fast formation within the protoplanetary disk lifetime involving a migration process (e.g. Bitsch et al. 2019;Izidoro et al. 2019;Lambrechts et al. 2019, and references therein). In particular, the migration process can lead to the capture of planet in mean motion resonant (MMR) chains (Cresswell & Nelson 2008). Most of the resonant chains are expected to break once the protoplanetary disk is dissipated as most of the systems are observed outside of MMR (Izidoro et al. 2017. Nevertheless, the period ratio distribution still shows signpost of this past history (Fabrycky et al. 2014). antoine.petit@astro.lu.se The study of the few remaining resonant chains is thus of particular interest to unravel the very early life of planetary systems.
The capture into MMR is a complex process that depends on the migration parameters, the planet masses and the particular resonance where capture happen (e.g. Mustill & Wyatt 2011;Batygin 2015). Nevertheless, we expect planets in resonant chains to have close to circular orbits due to the eccentricity damping during the migration (Cresswell & Nelson 2008). Such configurations have been observed for various systems observed through both radial velocities (RV) and with the analysis of Transit-Timing Variations (TTV, see e.g., Agol et al. 2005;Holman & Murray 2005). Systems where observations challenge this theoretical picture are of particular interest as they are the only way to probe the validity and/or the generality of theories.
K2-19 hosts three known transiting planets. The first two were reported (Armstrong et al. 2015) based on the photometry collected by the Kepler Space Telescope during the K2 operations (Howell et al. 2014). K2-19 b and c appeared to be close to the 3:2 MMR but it was not possible to conclude whether or not the pair was indeed in resonance.
Planets b and c sizes lies between Uranus and Saturn with respective radii of 7.0 ± 0.2 R ⊕ and 4.1 ± 0.2 R ⊕ and they orbit in respectively 7.9 days and 11.9 days. A third inner planet was detected by Sinukoff et al. (2016). Planet d orbits in 2.5 days and has a size similar to Earth with a radius of 1.11 ± 0.05R ⊕ . The recent observations by Petigura et al. (2020) of TTV for both planets b and c, as well as Radial Velocities (RV) measurements showing the reflex motion due to planet b, have given a precise set of orbital elements for the system. The 5% fractionnal uncertainties, are among the smallest for sub-Jovian exoplanets. Using a photodynamical model, they obtained a determination of the planet masses and orbital element down to a few percent. One puzzling aspect of this system is the moderate eccentricities of 0.2 with well-aligned apsides (∆ = 2 ± 2 • ) of planets b and c. Petigura et al. (2020) conclude from these observations that the system is very close to the 3:2 commensurability while not being resonant based on the analysis of the classical resonant angles. K2-19 system's architecture is thus puzzling from a dynamical point of view as no clear mechanism is identified to explain its stability. Indeed, Petigura et al. (2020) report that the system is stable in numerical simulations but is strongly AMD-unstable  and not protected by the resonance. Moreover, such a configuration is in tension with our current understanding of planet formation. Indeed, convergent migration within the protoplanetary disk leads to eccentricity damping and capture into MMR (Cresswell & Nelson 2008). The planets trapped into such state have low eccentricity (comparable to the protoplanetary disk aspect ratio 0.05) and have their periapses antialigned (Batygin & Morbidelli 2013).
The dynamics of the system and its origin thus necessitate an in-depth study. From a dynamical point of view, the proximity to the resonance and the eccentric and aligned orbits require one to go beyond the simple study of the resonant angles as they may not be representative of the nature of the dynamics. Given that all the planets are within 0.1 au, tidal effects must also be considered.
In this paper, we revisit the dynamical study of the system orbiting K2-19. We show in section 2 that the two outer planets are indeed trapped into the 3:2 mean motion resonance despite being apsidally aligned and the resonant angles circulation. The first order resonant model explains all the dynamics properties discussed by Petigura et al. (2020). We also show that the inner planet is secularly coupled to the b-c pair. We then study in section 3 the effect of tidal dissipation onto the inner planet d. We show that due to the eccentricity forcing from the outer planets, planet d's orbit tends to decay while the outer planets circularize. The time scale for the system evolution is shorter than its lifetime rendering the configuration unlikely. We finally discuss in section 4 the tension between the observations and our theoretical understanding of the system history. In particular we highlight the constraints on the three planet best fit and discuss whether the TTVs might be affected by an unmodelled effect, including the presence of a planet not yet detected. Table 1. K2-19 system parameters from (Petigura et al. 2020).

Parameter
Value In this section, we reanalyse the dynamics of the best fit three planet solution given by Petigura et al. (2020). We show that the outer planets are indeed inside the MMR and that they are coupled secularly to the inner planet. We partially reproduce in Table 1, the orbital elements and planet characteristics from the best photodynamical fit from Petigura et al. (2020).

The 3:2 mean motion resonance
The analysis of the K2 photometric data makes it clear that K2-19 b and c are close to the 3:2 MMR. Being close to the 3:2 MMR means that the planet mean motions n k = 2π/P k (P k being the planets' orbital periods) satisfy the arithmetic relation 2n b − 3n c 0 . As a result, the motion of the planets are coupled and one cannot average over the fast motions to study the longterm orbital evolution. Instead, the classical approach to analyse resonant motions consists of averaging the planet interactions over the non-resonant angles to reduce the problem to a one degree of freedom problem that is integrable. For first order MMR such as the 3:2 resonance, d'Alembert relations (see Morbidelli 2002) impose that at first order, the resonant terms in the development of the perturbation depend on the angles where λ k and k are respectively the mean longitude and the longitude of the periapsis of planet k. There are two different combination of angles related to the 3:2 resonance. In principle, the interaction between the two terms should make the system not integrable. In reality, the system can be reduced to a one degree of freedom resonant system thanks to a constant of motion that appears after a canonical transformation (Sessin & Ferraz-Mello 1984;Henrard et al. 1986). The integrable approximation for first order MMR has been called the second fundamental model of resonance 1 (Henrard & Lemaitre 1983). The analytical derivation of the integrable model for two massive planets has been carried out by several authors (Henrard et al. 1986;Batygin & Morbidelli 2013;Deck et al. 2013;Delisle et al. 2014;Petit et al. 2017;Hadden 2019). It is obtained by an expansion to first order in eccentricity, averaging over the fast angle and a rotation of the two classical resonant coordinates is the angular momentum deficit (AMD, Laskar 1997) of the planet k, Λ k = m k √ Gm S a k , G is the gravitational constant and θ res = 3λ c − 2λ b . Following ), we also definẽ The rotation, first described in Sessin & Ferraz-Mello (1984), transforms the coordinates x b , x c into two complex coordinates y 1 and y 2 (we follow here the notations from Petit et al. 2017). The norm of y 2 is a constant of motion and the dynamics of y 1 are described by the second fundamental model of resonance (see Ferraz-Mello 2007, for a complete description of the dynamics). It is also worth pointing out that the total AMD of the system is given by C = I 1 + I 2 where I k = y kȳk . For the 3:2 MMR, the expressions of y 1 and y 2 are where γ = m b /m c 3.0 for K2-19 and the numerical coefficients come from the expansion of the resonant terms of the 3:2 resonance. More precisely, the coefficients are linear combinations of different Laplace coefficients evaluated at the exact Keplerian resonance. We give the analytical expression in appendix A and refer to Batygin & Morbidelli (2013); Petit et al. (2017) or Hadden (2019) for recent complete derivations of the Hamiltonian.
Without loss of generality, one may rescale y 1 and y 2 They are linear combinations of the eccentricities where Y 1 is roughly the eccentricity vectors difference and Y 2 the mass weighted sum.
In the case of a system with two massive planets, the real resonant angle corresponds to the argument of Y 1 . Hence, to determine if the system is indeed in resonance one should in principle verify if the variable Y 1 evolves within the resonant island shown in figure 1. In reality, for most of the resonant chains observed in exoplanet systems, the resonant angles ϕ k (eq. 1) are good proxies for the actual resonant angle and such an analysis is not needed (see below).
On the other hand, the actual value of Y 2 is less critical to determine whether or not the system is resonant because it does not affect directly the shape of the resonance (fig 1). It can also be shown (see appendix A) that within the limit of the first order model, Y 2 precesses at the same frequency as θ res . Hence Y 2 e ιθ res is almost a constant that we will notẽ In reality, secular terms are neglected in this approximation. Nevertheless, the evolution ofỸ 2 happens on a much longer time scale. The phase space dynamics of Y 1 are shown in Figure  1. We compute the phase portrait using the expression of the Hamiltonian given in the appendix (eq. A16). We also plot the distribution of Y 1 using the posterior distribution of Petigura et al. (2020). It should be noted that while the uncertainties on the eccentricities are of the order of 0.03 in this paper, Y 1 has a much more restricted spread, reinforcing the argument that the system is in a resonant state i.e., the data require that Y 1 is well-within the resonant island shown in Figure 1. Finally we plot the result of a numerical integration of the two planet system. We integrated multiple draws from the posterior distribution from Petigura et al. (2020) with REBOUND (Rein & Liu 2012) and the high order integrator SABA(10,6,4) (Blanes et al. 2013;Rein et al. 2019) and for 5000 orbits of planet b (roughly 108 years). We show a sample tragectory in Figure 1. All draws behaved in qualitatively the same way, i.e., showing the libration of the resonant variable Y 1 .
It is clear that the dynamics are resonant, moreover, we can see that all the trajectories lie very close to the center of the resonance. It suggests that the mechanism that led to the capture must have been gentle and dissipative (Batygin 2015). The most favored mechanism for the formation is through migration within a protoplanetary disk (e.g. Cresswell & Nelson 2008). However, such a scenario is incompatible with the large eccentricities as well as the apsidal alignment of planet b and c as discussed in Petigura et al. (2020).

Posterior distribution
Example of trajectory Figure 1. Phase portrait of the integrable approximate Hamiltonian in the Y 1 (eq. 6) space. The black lines correspond to circulating orbits, the green lines to resonant orbits and the red line is the separatrix. The red set represents the posterior probability distribution of Y 1 at time t = 2020 BJD − 2454833. The grey points corresponds to the result of one N -body integration of the two planet system for 5000 planet b orbits.
over a period of about 6 yrs. We show that this oscillation is well explained by the resonant dynamics.
In most of the systems that have been observed in resonant configurations, Y 2 is usually negligible with respect to Y 1 . When |Y 2 | |Y 1 |, the libration of the angles ϕ k (eq. 1) is a good proxy to test whether or not a system is in MMR. However, in cases where Y 2 cannot be neglected with respect to Y 1 , the transformation from this set of resonant variables to the classical orbital elements is not straightforward. This means that the angles ϕ k can circulate while the system is actually very well described by the resonant first order integrable model. Indeed, the inverse transformation from Y k to the complex eccentricities can be written as Since Y 1 oscillates around the resonant center, its argument is librating. In a very rough approximation, we can consider it as constant. On the other hand Y 2 has a constant norm and is rotating, with a frequency comparable to the one of θ res that can be approximated as where δ = 2n b /(3n c ) − 1, represents the distance to the exact Keplerian resonance. In this approximation,ẽ b e ι( b −θ res ) describes a circle centred on Y 1 / √ 1 + 1.31γ and of radius 1.22Y 2 / √ 1 + 1.31γ (a similar analysis can be done for planet c). The angle ϕ b = b − θ res librates if the complex plane  Figure 2. Dynamical evolution of the eccentricities and resonant variables of the numerical solution described in section 2.1. The quantities Y k are defined by eq. (6) and (7) and e app. k by eq. (11).
origin lies outside of this circle. From eq. (9), we see that if In the case of the best fit studied here, both of these conditions are fulfilled. As a result, we cannot use the classical angles to probe the resonance. To our knowledge, this is the first system observed where the resonance cannot be characterized thanks to the classical angles.
The 3:2 resonance can also explain the eccentricity dynamics and in particular the apsidal aligned configuration. Let us denote with Y 0 1 and Y 0 2 the initial conditions for Y 1 and Y 2 at t 0 , the evolution of the eccentricities can be approximated as where γ was replaced by its average value and we used the approximationẽ k e k . From eq. (11), we see that the oscillations of e c are about three times as large as the one of e b . The period of oscillations should be of order 6.3 yr. We plot in Figure 2 the evolution of the eccentricities of planet b and c, of |Y 1 | and |Y 2 | as well as the approximation from eq. (11). While there is a small discrepancy on the frequency (the error is of the order of 10%), the amplitude of the motion is well reproduced. Moreover, the agreement in the plane e k -∆ is very good.
The main observational evidence for the resonance comes from the TTV. While the eccentricities evolve with a period of about 6 yr, the main frequency in the observed TTVs from Petigura et al. (2020) is about 2 years, which is the period of the oscillation of the resonant variable Y 1 .

Secular evolution
On short timescales (i.e. comparable to the TTV baseline), the simple model described above gives a good description of the system dynamics. On longer timescales (more than a few kyr), the system is subject to orbital precession due to secular interactions. We integrate the same initial condition as in the previous section, but this time we add planet d to the system. The initial condition was drawn from the posterior distribution and the mass of planet d in this particular realisation 2 is 5.9 M ⊕ . The simulation is run for 10,000 years, general relativity and stellar oblateness slightly change the precession rate but are not included in the example shown.
We plot the eccentricities as well as |Y 1 | and |Y 2 | in Figure 3. We see that as in the case with only planets b and c, the eccentricities e b and e c evolve very rapidly with the roughly 6 yr period seen in section 2.2 while |Y 1 | is almost constant. We checked that the resonance is indeed preserved in this case during the whole integration. However, |Y 2 | is no longer a constant and there are large AMD exchanges between planet d and the b-c pair. Due to the smaller planet d mass and semi-major axis, its eccentricity rises to values around 0.37 and its mean value is 0.24. We conclude that even starting with a circular orbit, planet d is largely coupled to the two outer planets and therefore cannot be considered in isolation.

TIDAL DECAY DURING LONG TERM EVOLUTION
From the last section we have seen that K2-19's system is stable over long timescales. But until now we have only taken into account the purely N-body gravitational interactions. However, due to their eccentric orbits, the three planets are subject to tidal effects. Indeed, the change of orbital distance leads to friction inside the planet and thus energy dissipation. Tidal effects conserve the total angular momentum 3 (Goldreich & Soter 1966). The energy loss results in a circularization of the orbits and a decay of the semi-major axes.
2 The value is close to the average of the posterior distribution but results were qualitatively similar for other realizations (not shown here).
3 In this analysis, we neglect the planets spin as well as tides raised on the star, that influences the stellar's spin. The total orbital angular momentum is thus conserved Dissipatitve effects are of particular importance due to the estimated age of the system. Indeed, based on K2-19's rotation period, the star is older than 1 Gyr. Besides, an age of a few Gyrs is compatible with the rotation rate and effective temperature (David Trevor, private communication: the star is more slowly rotating so more likely older than similar temperature 1 Gyr old stars).

Low-eccentricity tidal migration
Pu & Lai (2019) proposed a formation mechanism for ultra short period planets (with an orbital period of around 1 day), that they called the low-eccentricity migration scenario (see also Mardling 2007;Laskar et al. 2012). They showed that, in a multiplanetary system with slightly eccentric outer planets, the inner planet migrates very efficiently inward up until the point where it is decoupled due to the precession induced by the star oblateness and general relativity or if the outer planets run out of AMD.
Since planet b and c interact secularly with planet d only through the variation of Y 2 , we will assume that the interaction can be reduced to a two planet case: planet d and an outer planet. This simplification does not affect the results since the resonance is not affected by a variation of Y 2 (see figure 3). We use the two planet model presented by Pu & Lai (2019) to compute the effect of tides on planet d's orbit. It should be noted that considering the two outer planet as a single one is a conservative assumption as their interactions could lead to a faster migration (Pu & Lai 2019).
Following the weak friction theory of equilibrium tides (Darwin 1880;Alexander 1973;Hut 1981;Pu & Lai 2019), the evolution of the planets' semi-major axes in presence of tides is given by the equation a k a k = − 1.9 × 10 −9 k 2,k ∆t L,k 100 s e k 0.02 where k 2,k is the tidal Love number, ∆t L,k is the tidal lag time and R k is the radius for planet k. We note that the decay timescale depends on the eccentricity squared. The decay is fast at moderate eccentricity and slows down as the orbit becomes more and more circular. Moreover, the tides become less important for the farthest planet because of the steep a −8 k dependency. In principle tidal dissipation in the two outer planets should be considered. However, the tidal dissipation in large planets is not well constrained (Ogilvie 2014). We thus follow the conservative assumption made in Pu & Lai (2019) to neglect tides affecting planets b and c.
The decay rate also depends on the two coefficients k 2,d , of order unity, and ∆t L,d . Using this parametrized formalism allows us to study the dissipation while remaining agnostic on the actual physical mechanisms at play. We show below that the results are compatible with a large range of values for the coefficients. As Pu & Lai (2019), we take k 2,d = 1. The tidal lag ∆t L,d is inversely proportional to the planet's quality factor Planet d orbital period (days) Figure 4. Planet d period evolution due to low-eccentricity tidal migration for 1,000 initial conditions drawn from the three planet best fit. The blue curve is the average value. The spread is explained by the range of tidal lag times explored, the uncertainties on planet d's mass and the system's initial total AMD. to a quality factor close to 170. Using values close to the Solar system terrestrial planets is common in the field and motivated by studies of the viscoelastic response of planets to tidal deformations . In our analysis, we choose to draw planet d tidal lag time from a log-uniform distribution with boundaries 50 s ≤ ∆t L,d ≤ 500 s.
We first detail our model for the long term evolution of the system and then discuss the possible outcomes. As shown in section 2.3, the eccentricity of planet d is driven by its secular coupling with the outer planets. For planet d, we replace in eq. (12) the eccentricity e d by a forced eccentricity e d,forced given by the secular coupling. The forced eccentricity can be estimated as a function of planets b and c eccentricities. We use eq. (40) from Pu & Lai (2019), where ν d,b and ω b,d correspond to secular interaction terms between b and d, and ω d,gr and ω d,tide are respectively the apsidal precession of planet d due to general relativity and tides. While planet d is close to its current position, a d 0.035 au, the ratio e d,forced /e b is close to 0.5. For shorter orbits (a d < ∼ 0.01 au), the ratio sharply decays as general relativity and tidal precession become significant and decouple planet d from the outer planets, effectively stopping the tidal migration.
Finally, the outer planets eccentricities evolve as planet d migrates because tides raised on planets conserve the total orbital angular momentum 4 . One can estimate e b as a function of the new semi-major axis of planet d, and the initial values for e b , a d and a b where C 0 is the total initial AMD of the system and Λ k,0 = m k Gm S a k,0 is the initial circular angular momentum of planet k. To obtain expression (15), we assumed that planets b and c see no variation of their semi-major axis due to tidal migration and that e c and e b were equal, which is reasonable in first approximation due to the resonance. By inserting (14) and (15) into (12), we obtain a differential equation for a d that gives results comparable to a secular complete integration (Pu & Lai 2019). We draw 1,000 initial conditions from the posterior distribution of Petigura et al. (2020). We plot in figure 4 the evolution of the period of planet d over 10 Gyr. The individual evolutions are in red while the thick blue line corresponds to the averaged value.
The final orbital periods extend from 0 days (where planet d would be consumed by the star) to almost 2 days. It should be noted that for periods smaller than 0.4 days (semi-major axis of about 0.01 au), the decay is expected to be faster due to stellar tides that are neglected in this analysis. The typical outcome is the formation of an ultra short period planet with an orbit of about a day. The final orbit of planet d is mainly determined by the mass of planet d and the initial AMD. For larger AMD or smaller mass m d , the final orbit is shorter.
We plot in figure 5, the resulting distribution of planet b eccentricity 5 and of planet d period at t = 0 Myr, t = 500 Myr and at t = 2 Gyr. We see that even after a few hundreds of Myr, the decay of the eccentricity and orbital period are general and significant. After 2 Gyr, the eccentricity of planet b is smaller than 0.1 in 85% of the simulations. In the simulations where planet d does not migrate up to the star, the orbital decay is stopped because the AMD reservoir has been emptied, i.e. the outer planets' orbits have been circularized.
More importantly, the final state is reached within a few hundred Myr. We define the half-decay time T hf as the time such that planet d has undergone half of its orbital decay over 10 Gyr. The median half decay time is T hf = 472 Myr and 80% of the initial conditions have T hf < 1 Gyr.

Tidal decay in the past history of the system
We also considered the case where planet d started on a wider orbit and is currently experiencing tidal decay. We run the same model but with planet d starting with a period of 3.5 days, correcting the total AMD such that the system keeps the same total angular momentum. This initial period is the largest one before planet d's orbit crosses planet b's. After 2 Gyrs, only 30% of the systems have planet d with a period larger than 2.4 days and e b > 0.1. The systems compatible with today's observations give a constrain on planet d tidal lag. The 1σ upper limit is ∆t 86% L,d = 130 s (which corresponds to a quality factor of 180 at the current orbit). The constraints become stronger if the system is older than 2 Gyr.
It results that it may be possible that the observed system is on its way to circularize in the next billion of years. Such a scenario necessitates a very particular initial configuration where planet d is originally on an orbit at the limit of instability.

TRULY ECCENTRIC? TENSIONS BETWEEN THE OBSERVATIONS AND THE THEORY
The observations for this system are very precise and from multiple sources. The planets were detected thanks to the photometry from the K2 campaign (Armstrong et al. 2015;Sinukoff et al. 2016). The orbital elements and masses are also constrained by RV data and TTV obtained 2 years after the K2 campaign. Petigura et al. (2020) performs a photodynamical fit that forward model the lightcurve. The eccentricities and periapsis argument are parametrized as { √ e cos ω, √ e sin ω} such that the prior is uniform in eccentricity (Eastman et al. 2013). Given the posterior distribution obtained, a close to circular, three planet model is strongly ruled out. Petigura et al. (2020) already noticed that the architecture of the system was puzzling depsite it being stable in numerical simulations.
However, it appears from the previous section that the tidal dissipation make the current system's architecture harder to explain. Besides, the formation of the system remains unexplained. So one needs to explain how to fit the observations while ensuring that the system configuration can be observed after a few Gyr.

The formation challenge
We showed in section 2.1 that K2-19 b and c are trapped in the 3:2 MMR with a small libration amplitude. Capture into MMR generally emerges from dissipative effects leading to convergent migration (Batygin 2015). The most common mechanism is migration within the protoplanetary disk (e.g. Cresswell & Nelson 2008;Pichierri et al. 2019). Capture can also occurs due to convergent tidal migration (Papaloizou et al. 2018). However, both disk and tidal migration show shorter timescales for eccentricity damping than for change in orbital period. Systems are thus capture close to circular orbits, an increase of the eccentricity while the system is in the resonance typically leads to an anti-aligned configuration as pointed out by Petigura et al. (2020).
In order to explain the present configuration, one has to imagine a mechanism where the planet's eccentricity vectors are not damped to zero but to a common value. In this case, the capture can occurs in the aligned configuration since the resonant dynamics and the capture mainly depend on Y 1 (Batygin 2015, eq. 22 and fig. 7). Migration in eccentric disk has been studied theoretically (Papaloizou 2002), but eccentric disks arise only in the presence of large planets (Teyssandier & Ogilvie 2016) or for circumbinary disks. Moreover, eccentricities close to 0.2 leads to an outward migration (D'Angelo et al. 2006), which seems at odds with the short period of planet b and c. In anyway, it is clear that the formation of the three planets around K2-19 remains a challenge that does not fit the classical scenarios.

New constraints from tidal dissipation
We have shown in section 3 that taking into account tides is critical to understand the long term evolution of this system. Indeed, the best fit orbital solution is stable in the presence of purely gravitational interactions over long timescales. However, the secular coupling between planet d and the resonant pair (see 2.3) leads to a strong tidal dissipation in the inner planet. Over a few hundreds Myr, the outer planets' AMD is depleted and planet d experiences a period decay. If the system had truly formed as we see it today, we would expect to observe the outer planets on circular orbits and the inner one on a shorter orbit.
Yet, the current system configuration might be explained in presence of tidal decay. However, it implies that planet d orbital period was originally larger ( up to 3.5 days) and that the three planets started on eccentric orbits with planet d and b at the limit of the orbit crossing. It also requires that the tidal lag time of planet d is smaller than 130s (or that the quality factor is larger than 180), a value that is not incompatible with our understanding of tides in rocky bodies but that still gives a constraint on the dissipation rate. Lower quality factors are strongly ruled out. Such a fine tuning of the original configuration is necessary to explain the observed architecture if the system only hosts three planets.
In this work, we have not taken into account the dissipation in planets b and c because of the poor constraints on tidal dissipation in gaseous planets. However, works on on tidal dissipation for systems in MMR (e.g. Delisle et al. 2012Delisle et al. , 2014 have shown that the resonance is often broken before the planets circularization. Such studies could give an additional constraint on the system.

A fourth planet?
When it is hard to reconcile the observations to the theory, it is common to explore the possibility that the planets motion can be perturbed by an unseen companion. From Le Verrier's work on Neptune to the recent Planet 9 hypothesis , this approach is historically tied to the progresses in our understanding of the Solar System because of the very precise constraints on the planets' motion. However, the method have been applied with success also in the context of exoplanet dynamics. One can cite the example of Kepler-56 where the 40 deg obliquity of the two transiting planets is due to a distant planetary companion (Huber et al. 2013;Otor et al. 2016). But the most obvious application to exoplanets is the TTV method itself that allows us to find non-transiting planets, most of the time in resonance with a transiting one (Kepler-88, the "King of TTVs" is an example of a system where a non-transiting planet perturbing the motion of Kepler-88b, Nesvorný et al. 2013).
The photodynamical fit strongly reject solutions with three planet on circular orbits. Nevertheless, it remains possible that the model is misspecified, which would be the case if, for instance, the TTVs are affected by a fourth planet in the system. In this case, the orbits of the observed planets b, c and d might be close to circular at the expense of the addition of another planet in the resonant chain. Indeed, a few Earth masses planet trapped in an inner resonance with planet b can have significant TTVs contribution while being non transiting or even not being detected due to its small radius. We also point out that the TTV coverage of Petigura et al. (2020) is sparse, especially compared to TTV datasets from the prime Kepler mission. This dataset had timing measurements over three distinct epochs (the K2 campaign and two sets from Spitzer ). We expect sparse datasets to be more susceptible to model misspecification errors and encourage additional transit time measurements.
While it will necessitate more transit observations in the future to verify such a claim, it should be noted that it is possible to add a planet between planet d and b without destabilizing the system. Indeed, we show the positions of the planets alongside the resonances 2:1 and 3:2 with planet b in figure 6. Another planet in 3:2 or 2:1 resonance with planet b is consistent with the 'peas in a pod' pattern observed in the architecture of the Kepler systems (Weiss et al. 2018). Moreover, assuming this unseen planet is in a 2:1 resonance with planet b leads to a period ratio of 1.58 with planet d, which is just wide of the 3:2 resonance. In other words, the whole system could have been placed into a four planet resonant chain during its formation before tidal effects broke the resonance with the innermost planet as it has been proposed by Millholland & Laughlin (2019) or Pichierri et al. (2019).

CONCLUSIONS AND DISCUSSIONS
In systems such as K2-19, the precise orbital parameters obtained by TTVs allow for rich dynamical studies that can 0 2 4 6 8 10 12 Orbital period (day) d b c 3:2 2:1 Figure 6. Architecture of the system K2-19 including the inner 2:1 and 3:2 MMR with planet b. The size of the planets is proportional to their radius.
provide a lot of insight onto the system history and configuration. Such systems act as laboratories to test theories of formation and dynamical evolution.
Following their observations, Petigura et al. (2020) conclude that K2-19's two outer planets are very close to the 3:2 MMR but due to their large aligned eccentricity, the classical resonant angles are not librating. A system with such a configuration is a challenge for classical planet formation scenarios. Indeed, planets close to MMR tend to be captured due to type I migration during the disk lifetime (Cresswell & Nelson 2008). Besides, migration is very effective to damp planet's eccentricities.
We have shown in section 2.1 that by considering the true resonant variables of the system that is a combination of the two complex eccentricities, the system rapid dynamics are well explained by the integrable first order model. Besides, the very small libration amplitude of the resonant variable ( fig. 1) strengthens the classical scenario of a smooth capture through disk migration (Batygin 2015). We also see that the resonant variable is more constrained than the two eccentricities as it is expected for systems presenting TTVs (Hadden & Lithwick 2016). In particular, we want to highlight that in the context of orbital fit with such large eccentricities, the classical methods to spot a MMR such as monitoring the classical resonant angles (eq. 1) is not reliable. A resonant configuration can also exist even if the two orbits are not anti-aligned.
We then studied the entire system rather than only the two outer more massive planets. The fitted configuration leads to a very strong secular coupling between the inner terrestrial planet d and the pair b and c. The coupling does not disrupt the MMR but leads to eccentricities of the order of 0.35 for the inner planet. In particular, this planetary system's architecture is long-lived in the presence of purely N-body interactions.
Because planet d orbits in about 2.5 days around its host star, it is subject to large tidal effects. The tides raised on planet d tend to circularize its orbit at the expense of a period decay. Following the low-eccentricity migration model presented in (Pu & Lai 2019), we show that the two outer planets transfer their AMD to the inner planet through the secular coupling. As a result the inner planet continues to decay up until the point where tidal and general relativity precession decouple it from the outer planets or if the outer planets run out of AMD. The typical timescale of the process is less than 500 Myr whereas the system is expected to be billions of years old. The system's current architecture remains compatible with the tidal decay if planet d started in a longer orbit, at the limit of orbit crossing with planet b. Note that in the absence of planet d, it would be much harder to rule out the aligned eccentric resonant configuration due to the poor constraint we have on the dissipation onto gaseous planets.
Even though the photodynamical fit gives a configuration that can only be maintained for a short amount of time with respect to the system lifetime, the observations have to be explained. While biases in eccentricity determination are known and have been quantified in RV observations (Anglada-Escudé et al. 2010;Hara et al. 2019), no such study has been carried for TTVs systems. One could speculate on the fact that another undiscovered planet in the system can affect the outer planets periods. Indeed, in the picture of the "peas in a pod" systems (Weiss et al. 2018), the system is compatible with the presence of another planet between planet d and b. If such a non-transiting planet could also be in resonance with planet b and c and lead to TTV unaccounted for. However, the presence of an eventual fourth planet would need to be confirmed by more measurements of planet b and c transits.
The detailed analysis of the K2-19 system has revealed even richer dynamics than originally reported. Taking into account the dissipative effects also appears to be crucial for the understanding of the system history. Nevertheless, the system formation remains mostly unexplained. Future photometric or RV monitoring will be crucial to unveil the nature of this system. Γ has been called in the literature the scaling factor Sessin & Ferraz-Mello (1984), G is the total angular momentum. Because of the d'Alembert relations, the Hamiltonian does not depend explicitly on θ res (i.e. the angular momentum is conserved). The two resonant angles are k − θ res (see eq. 1). We then expand the perturbation εH 1 in power series of x k and Fourier series of θ Γ and only keep the first order terms. The next step consists in averaging the motion over the fast angle θ Γ . This operation is also a canonical transformations and the new coordinates are ε close to the old ones. They correspond to the average coordinates over a Keplerian orbit. The resulting Hamiltonian no longer depends on θ Γ . As a result, Γ is a constant of the averaged system. Note that the inverse transformation of eq. (A4) allow to express Λ k as a function of the total AMD, C and the constants of motion Γ and G.