Stability Constrained Characterization of the 23 Myr Old V1298 Tau System: Do Young Planets Form in Mean Motion Resonance Chains?

Planets c and d transited nine and five times during the K2 observations, and four and three times when observed by TESS, roughly 6.5 yr (≈200 orbits) later, respectively. This effectively provides two "instantaneous" measurements of their physically meaningful period ratio: 1.503 ± 0.0001 as measured in 2015 with K2 (David et al. 2019a), and 1.503 ± 0.0005 as measured in 2021 with TESS (Feinstein et al. 2021).

The simplest explanation for the consistent period ratio measurements in the two epochs is that the orbital periods are approximately constant in time, though more general oscillating solutions one might expect near MMRs (see Chapter 8 of Murray & Dermott 1999) are also possible. Most of the techniques for speeding up Markov Chain Monte Carlo parameter estimation methods (e.g., Foreman-Mackey et al. 2013) do not work for sampling the wide and disconnected sets of masses and orbital parameters consistent with two observed period ratios. We therefore opt for a simple rejection sampling approach (e.g., Price-Whelan et al. 2017) where we accept a set of masses and orbital parameters based on the likelihood L of the period ratio passing through the observations: ²

$$\begin{eqnarray}&&\mathrm{log}L=-\displaystyle \frac{1}{2}\sum _{i}\displaystyle \frac{{\left[{p}_{\mathrm{obs}}({t}_{i})-{p}_{\mathrm{model}}({t}_{i})\right]}^{2}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 1 }$$

where the t_i are the times of the (two) observations, σ_i the observational uncertainties (assumed Gaussian), p_obs are the observed period ratios, and p_model are the period ratios as determined using the WHFast integrator (Rein & Tamayo 2015) in the REBOUND N-body package (Rein & Liu 2012).

In order to sample a wide range of resonant and near-resonant configurations consistent with the observed period ratios, while reducing the number of parameters and aiding in their interpretation, we exploit analytical models of the planetary dynamics in and near MMRs. In particular, we sample four parameters: the planet pair's total planet–star mass ratio μ = (m_c + m_d )/m_⋆, the equilibrium combined eccentricity e_forced parameterizing the strength of the 3:2 MMR, the combined free eccentricity e_free, measuring how far the planet pair is from the equilibrium value, and the resonant angle ϕ, which measures the azimuthal angle at which conjunctions between the two planets occur, and sets the phase of the period ratio oscillations. The fact that two masses and twelve orbital parameters for the two planets can be approximately distilled into four parameters determining the period ratio evolution is not obvious, so we provide more details in the Appendix. We perform these and all other calculations of resonant parameters in this paper using the open-source celmech package. ³

We sampled ${e}_{\mathrm{forced}}^{\mathrm{cd}}$ and ${e}_{\mathrm{free}}^{\mathrm{cd}}$ uniformly from 0 to 0.2 (i.e., nearly orbit-crossing), and the total planet–star mass ratio μ log-uniformly. These planets may be significantly inflated (David et al. 2019a), so we adopt a low-mass bound of μ at 3 M_⊕ based on measurements of the least dense exoplanet known—Kepler-51b (Steffen et al. 2013; Masuda 2014). For the upper bound, we chose 1 M_Jup to conservatively account for the RV upper bounds (Suárez Mascareño et al. 2021).

We sampled 60 million configurations, accepting 9693. This posterior size is adequate to explore our parameter space, and is shown in Figure 1, where the dashed histograms are the rejection posteriors. The color-coding and solid histograms incorporate stability constraints and are discussed in Section 3.

Zoom In Zoom Out Reset image size

We divide our configurations into ones that are inside the 3:2 MMR on the left half of Figure 1, and ones that are not (≈78%) on the right (see Appendix A.3 for details). We see that in both the resonant and nonresonant cases, the observed period ratios for planets c and d significantly favor lower masses. In the resonant cases, the forced eccentricities parameterizing the strength of the MMR are higher than in the nonresonance configurations. This reflects the fact that for most of the configurations on the right side of Figure 1, the MMR is weak enough that there is no resonant region at all (Appendix A.3). The oscillations around the equilibrium, parameterized by e_free are small in both cases.

We also show the period ratio evolution for a random selection of 100 solutions in the bottom panels with the two observed values in blue. Many of the nonresonant configurations on the right exhibit flat period ratios because the masses and orbital parameters render the MMR weak. For the resonant cases where the MMR is strong, we get finite-amplitude oscillations in the period ratio (the bottom left panel of Figure 1) even if we restrict e_free = 0, which by definition should correspond to a configuration at the resonant equilibrium with a constant period ratio. This discrepancy is due to the fact that the analytical MMR model, and its prediction for the location of the resonant equilibrium, is only approximate. Given that the two observations in blue (the bottom panels of Figure 1) are consistent within very narrow error bars (plotted), improved or numerical models would provide a stronger preference for nearly flat period ratio solutions with e_free ≈ 0. Since our goal is simply to explore a representative range of resonant and near-resonant configurations for planets c and d in the presence of a massive planet b, we do not pursue this further. A third set of transits for planets c and d would provide a strong and valuable constraint on their resonant configuration, and for whether the period ratio is constant with time.

For the resonant cases (the left side of Figure 1), the configurations most vulnerable to instabilities due to perturbations from other planets are those near the boundaries of the MMR, i.e., near the separatrix (e.g., Rath et al. 2021). The separatrix corresponds to a particular value of e_free, or equivalently to the black "cat's-eye" boundary in the plane spanned by ϕ and the period ratio deviation from 1.5 shown in the left panel of Figure 2 (see the Appendix for details). To account for the fact that the size of this resonant region varies with μ and e_forced, for each of the resonant configurations, we use celmech to normalize the period ratio deviation by the value at the separatrix so that the MMR boundary for all the resonant configurations plotted in Figure 2 extends from [−1, 1]. We see that while a diffuse distribution of configurations fill the "cat's eye" resonant region, most configurations cluster closer to the equilibrium (equivalently at low e_free), with higher mass solutions closer to the equilibrium.

Zoom In Zoom Out Reset image size

Our approach requires closer examination for this particularly young system, given the possibility that planetary systems may generically be unstable, undergoing dynamical instabilities throughout their lives (e.g., Laskar 1996; Pu & Wu 2015; Volk & Gladman 2015; Izidoro et al. 2017; Bitsch et al. 2019; Izidoro et al. 2021). Nevertheless, in that picture, one might expect a system observed at an arbitrary snapshot in time to always have a time to instability comparable to its present lifetime (Laskar 1996). Thus, even if planetary systems are ultimately unstable, one can still require stability on timescales much shorter than their age, since otherwise instabilities would be happening all the time.

The SPOCK FeatureClassifier (Tamayo et al. 2020) provides a probability of stability over 10⁹ orbits of the innermost planet, corresponding to ≈22 Myr for V1298 Tau. For this young system, that timescale is close to its age, and thus too stringent a timescale on which to require stability.

However, most instabilities occur early, and we do not expect this to significantly bias our results. About 96% of our tested configurations had SPOCK probabilities ≤0.2. We checked the median instability times for these low-probability configurations with SPOCK's deep learning classifier (Cranmer et al. 2021), and found 97% had predicted lifetimes ≤10⁷ orbits (≲0.2 Myr). This preference toward short lifetimes also matches expectations from the analytic investigation of these instabilities by Tamayo et al. (2021b).

For the nonresonant cases on the right side of Figure 1, stability sets a 99.7th percentile upper limit for the combined mass of planets c and d of μ < 0.35 M_Jup, slightly more constraining than the upper limits from the RV observations (Suárez Mascareño et al. 2021). Stability of resonant configurations (the left side of Figure 1), on the other hand, allows for values up to 0.95 M_Jup (99.7th percentile). Moreover, SPOCK shows a preference toward particular values of μ (peaks in the solid histogram for μ, the lower right triangle plot panel on the left side of Figure 1).

These modes are apparent in Figure 2, where in the right panel we set the opaqueness of each configuration according to its probability of stability according to SPOCK. As expected, stability preferentially eliminates configurations near the separatrix, and selects for setups closer to the resonant equilibrium. While we do not fully understand why stability is clearing out the annulus of configurations between the turquoise and purple rings in the left panel of Figure 2, we speculate that this may be due to chaos induced by secondary resonances, i.e., with the frequency of oscillation around the equilibrium point (see Rath et al. 2021). A third measurement of the period ratio between planets c and d would strongly constrain the amplitude of the period ratio oscillations, and provide a valuable constraint on these planets' resonant configuration.

Figure 3 plots stability constraints on planet b, taking our 10⁵ samples drawn from the mass and orbital eccentricity distributions reported by Suárez Mascareño et al. (2021). As with Figure 1, red points are configurations that went unstable in N-body integrations within 10⁴ orbits, while the surviving black points have their opacity set according to their SPOCK probability of surviving for the next ≈20 Myr.

Zoom In Zoom Out Reset image size

We obtain a 99.7th percentile upper limit for m_b of 1.07 M_Jup. Even though this upper limit is only ∼0.1 M_Jup lower than the Gaussian 3σ RV limit, we find that even if we extend our mass priors to higher values, masses above this limit are unstable at any eccentricity. Therefore, coincidentally, stability and RV observations yield very similar upper limits on planet b's mass. We also note that the median mass estimated through RV observations of 0.64 M_Jup (Suárez Mascareño et al. 2021) is within a factor of 2 from the stability limit even at zero eccentricity (Figure 3).

This implies strong stability constraints on planet b's orbital eccentricity, e_b . We obtain a 99.7th percentile upper limit of 0.17 for the e_b , which is half of the corresponding RV limit. Finally, we find that the constraints on planet b in Figure 3 do not depend strongly on the parameters for planets c and d (e.g., the bottom two rows of Table 2 comparing constraints for resonant and nonresonant configurations for c and d). This renders our approach of examining the configuration of planets c and d separately from those on planet b particularly valuable.

Table 2. 99.7th Percentile Upper Limits

Parameter	Inside Resonance	Outside Resonance	RV (3σ) Reported
μ (MJup)	0.94	0.35	0.55
${e}_{\mathrm{forced}}^{\mathrm{cd}}$	0.17	0.06	⋯
${e}_{\mathrm{free}}^{\mathrm{cd}}$	0.02	0.02	⋯
mb (MJup)	1.07	1.07	1.21
eb	0.17	0.17	0.34

Note. 99.7th percentile upper limits. The second and third columns correspond to stability limits in the cases where planets c and d are inside and outside the 3:2 MMR, respectively (the left and right sides of Figure 1). The last column are the Gaussian 3σ RV limits reported by Suárez Mascareño et al. (2021). Bootstrap resampling of the stability posteriors gives errors on the order of ∼10⁻³ for each parameter.

Download table as:  ASCIITypeset image

A challenge to establishing whether or not V1298 Tau is in an MMR chain is that the theory for MMR chains (e.g., Delisle 2017; Siegel & Fabrycky 2021) is less developed than for MMRs between a single pair of planets. We therefore take two approximate approaches.

First, we ask how many stable configurations place planets d and b (observed period ratio ≈1.946) inside the 2:1 MMR. In particular, when counting the fraction of configurations in the 2:1 MMR, we weight each one by its corresponding SPOCK probability of stability p_i (Tamayo et al. 2021a):

$$\begin{eqnarray}&&P(\mathrm{MMR}\,\mathrm{chain})=\sum _{i\ {in}\ 2:1}{p}_{i}/\sum _{i}{p}_{i}.\end{eqnarray} \tag{ 2 }$$

This yields P(MMR chain) ≈ 1%. Even if future observations yield lower masses and eccentricities for planet b, the 2:1 MMR region lies at higher masses, making the system unstable. The fact that P_b /P_d ≈ 1.946, far from 2, is only in the 2:1 resonance region for the highest masses and eccentricities makes it even less likely that planets d and b are in the 2:1 resonance.

However, this analytical two-planet model ignores the gravitational effects of the third planet c, which can be important. We therefore instead check for libration of the three-body angle ϕ_chain = λ_c − 2λ_d + λ_b over timescales of 10³ orbits of planet c, where the λ are the planets' mean longitudes. Using this libration condition to select configurations in the 2:1 MMR for the first sum of Equation (2), we obtain P(MMR chain) ≈0.02%. Both tests, therefore, rule out an MMR chain configuration for V1298 Tau at ≳99% confidence.

Similar to the stability step performed for planet b, we sampled the eccentricity and mass distributions of planet e from Suárez Mascareño et al. (2021) and the period distribution from Feinstein et al. (2021). We drew 10⁵ samples from the stable configurations already with planet b, and we conducted uninformative and RV prior analysis similar to those we performed with planet b. In both cases, we found the inclusion of planet e does not affect the results of the inner planets, nor does stability constrain the period, mass, or eccentricity any better than what is already reported. The fact that plausible parameters for planet e do not affect our results justifies our choice to simplify the parameter space by excluding it.

The dynamics near first-order MMRs depend on the orbital period ratio, masses m_i , and eccentricity vectors e _i (pointing in the direction toward the pericenter with a magnitude given by the eccentricity), where the subscript i indexes the planet. Two approximations simplify how the period ratio evolution depends on system parameters. First, Deck et al. (2013) show that the period ratio dynamics depend approximately only on the total mass of the planet pair divided by the stellar mass, μ = (m₁ + m₂)/m_⋆ (and not the mass ratio m₁/m₂). For concreteness, we set planets c and d to have equal masses, but we have checked that varying the mass ratio does not significantly affect our results.

Second, the period ratio evolution depends primarily on a single linear combination of the eccentricity vectors e ₋ (Sessin 1983; Henrard 1984; Hadden 2019),

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{-}=\sqrt{2}\ \displaystyle \frac{{{fz}}_{1}+{{gz}}_{2}}{\sqrt{{f}^{2}+{g}^{2}}}\approx {{\boldsymbol{e}}}_{2}-{{\boldsymbol{e}}}_{1},\end{eqnarray} \tag{ A1 }$$

where f and g are Fourier coefficients in the disturbing function expansion of the interplanetary potential (e.g., Hadden 2019). The approximation that e ₋ is nearly the vector difference e ₂ − e ₁ is excellent for all first-order MMRs except the 2:1 (in which the relationship between f and g is altered by additional indirect terms that need to be considered), and can thus be thought of as a relative or antialigned eccentricity. This transformation additionally defines a combined pericenter ϖ₋, which approximately specifies the longitude at which the two orbits come closest to one another.

The remaining linear combination of eccentricities, which can be expressed as an approximately conserved center-of-mass eccentricity (e.g., Hadden 2019), to excellent approximation does not affect the period ratio evolution for a single pair of planets on short timescales. However, larger values of this center-of-mass eccentricity do affect stability in multiplanet configurations, by leading to larger secular oscillations in e₋ on longer timescales (Tamayo et al. 2021b). In order to limit the parameter space, we set the center-of-mass eccentricity to zero as a best-case scenario for stability. Our derived upper limits on masses and orbital parameters in this best case are thus conservative and reliable.

For a period ratio near an MMR commensurability, different masses, eccentricities, and pericenter orientations will lead to oscillations in the period. However, there always exists a resonant equilibrium corresponding to a particular combination of the above parameters, where the period ratio remains constant.

To a good approximation, for a given j: j − 1 MMR, the equilibrium, or forced, eccentricity e_forced is only a function of the equilibrium period ratio and μ (e.g., Hadden 2019). In addition the equilibrium resonant angle ϕ = j λ₂ − (j −1)λ₁ − ϖ₋, which approximately measures the location at which conjunctions occur, is always at ϕ = π. This corresponds to conjunctions occurring at the location where the two orbits are farthest from one another.

Initial conditions near, but not at, the resonant equilibrium will oscillate around it. The antialigned eccentricity e₋ can thus be profitably decomposed into free and forced components (e.g., Murray & Dermott 1999),

$$\begin{eqnarray}&&{e}_{-}={e}_{\mathrm{forced}}+{e}_{\mathrm{free}}\end{eqnarray} \tag{ A2 }$$

In this picture, e_forced is the equilibrium value at which the antialigned eccentricity and period ratio would remain constant (if ϕ = π), and the free eccentricity e_free is how far away the eccentricity is from the fixed point. Additionally, the forced eccentricity is a conserved quantity involving both e₋ and the period ratio's deviation from the resonant value (in this case P_d /P_c = 3/2), which lets one switch back and forth between the two variables.

Because various configurations with different resonance strengths can have very different values of e_forced, but the equilibrium period ratio is always approximately 3/2, we choose in Figure 2 to plot the period ratio deviation δ P. The center at (ϕ, δ P) ≈ (π, 0) is the equilibrium, corresponding to different nonzero values of e_forced for each configuration (point). Points near the equilibrium approximately execute circles around the fixed point, allowing for a simple scalar distance from the fixed point to parameterize the oscillation amplitude. At larger distances from the resonant fixed point, the trajectories become more deformed, and can even oscillate around a different fixed point. Since the distance from the equilibrium therefore varies for general, noncircular trajectories, we define e_free to be the value for the trajectory (through Equation (A2)) when the resonant angle ϕ crosses through π, and we convert it to a period deviation with celmech to make Figure 2.

In order to account for the phase of oscillations to match the observed period ratios (the bottom panels of Figure 1), we initialize the system with total mass ratio μ, and e_forced and e_free at ϕ = π (i.e., along a vertical line in Figure 2). We then integrate the orbits using the REBOUND N-body integrator (Rein & Liu 2012) for a uniformly drawn time ΔT, which we marginalize over when presenting the posterior distributions in Figure 1.

While libration of the resonant angle ϕ is often used as a test of whether a pair of planets is in resonance, this distinction is not physically meaningful when the MMR is weak (i.e., in these cases there is no dynamically meaningful difference in behavior between a configuration where ϕ circulates, and one where it librates with a large amplitude). A more meaningful distinction is to separate this question into two parts.

First, there is a threshold strength for first-order MMRs (parameterized by μ and e_forced), beyond which a separatrix appears (e.g., Henrard 1984; Deck et al. 2013). This separatrix trajectory is (in the analytic approximation of the MMR model) an infinite-period trajectory (black "cat's-eye" curve in Figure 2). Configurations inside the separatrix (small e_free) are resonant, and cleanly separated from nonresonant trajectories (large e_free) on the outside. This is a physically important boundary because trajectories near the separatrix are most susceptible to chaos under perturbation (e.g., Lichtenberg & Lieberman 1992; Rath et al. 2021). When the MMR is weak, i.e., below the threshold MMR strength, there is no separatrix, and there is no dynamically meaningful boundary between resonant and nonresonant configurations.

We therefore use the celmech package to calculate whether a separatrix exists or not for each of the orbital configurations in Figure 1. If it does, we use the fact that the numerical value of the (conserved) Hamiltonian for a resonant trajectory (inside the separatrix) is always smaller than the value on the separatrix. On the left side of Figure 1, we then plot all the configurations where (a) a separatrix exists and (b) the value of the Hamiltonian places it inside the resonant region. All remaining configurations are then plotted on the right side of Figure 1, i.e., (1) cases where the MMR is so weak there is no separatrix and (2) cases where a separatrix does exist, but the system is outside the resonant region.