The Six-planet Resonant Chain of HD 110067

HD 110067 is the brightest star known to have six transiting planets. Each adjacent pair of planets has a period ratio that is nearly equal to a ratio of small integers, suggesting the planets are in a chain of mean-motion resonances, but the limited time span of the available data has prevented firm conclusions. Here, we show that the requirement of long-term dynamical stability implies that all six planets are very likely to form a resonant chain. Dynamical simulations of nonresonant systems with initial conditions compatible with the available data almost always suffer an instability within 25 Myr (∼0.3% of the system’s age). Assuming the system is in resonance, we place upper limits on the planets’ eccentricities and lower limits on the masses of the planets that have not yet been measured. We also predict the characteristics of transit timing variations and the values of the three-body libration centers.


INTRODUCTION
A "resonant chain" is a planetary system in which multiple adjacent planet pairs are locked in meanmotion resonances (MMRs).Although resonant chains are believed to be natural outcomes of slow migration in protoplanetary disks (e.g., Henrard 1982;Lee & Peale 2002;Kley et al. 2005;Terquem & Papaloizou 2007), only ∼ 1 % of known exoplanetary systems are in a resonant chain configuration (Fabrycky et al. 2014).Because external perturbations tend to break resonances, resonant chains are believed to have avoided major disturbances since their formation and thereby provide a view of the immediate outcomes of the planet formation process.
Gliese 876 was the first known example of an exoplanetary resonant chain, with three planets -two giant planets and a super-Earth-sized planet -having periods nearly in the ratios 4:2:1 (Rivera et al. 2010).Later, Kepler-223 was found to host a four-planet resonant chain consisting of sub-Neptune-sized planets (Mills et al. 2016).Starting with the innermost planet and moving outward, the ratios of the periods of adjacent planets are nearly 4:3, 3:2, and 4:3.Nearly commensurable orbital periods do not guarantee, however, that a planetary system is trapped in a resonant configuration.Even planet pairs located within 1 % of a first-order MMR (i.e., |∆| = | j − k j Pout Pin − 1| < 0.01, where j and k are the index and order, respectively, of a j:j − k MMR) typically have circulating two-body resonant angles (Hadden & Lithwick 2017).Mills et al. (2016) showed that the planets of Kepler-223 form a resonant chain based on the analysis of four years of observed transit timing variations (TTVs).Specifically, they demonstrated that the systems' threebody resonant angles (also termed "Laplace angles") are librating rather than circulating.TRAPPIST-1 is another well-characterized resonant chain system, which consists of seven approximately Earth-sized planets, the largest number of planets known in any resonant chain (Gillon et al. 2017;Luger et al. 2017).It is also possible for a resonant chain to consist of only a subset of a system's planets; for instance, the outer five planets of TOI-178 form a resonant chain, but the innermost planet is not part of the chain (Leleu et al. 2021).
HD 110067 is a K0 dwarf star with mass M * ≈ 0.8 M ⊙ and age 8 ± 4 Gyr.With a V magnitude of 8.4, it is the brightest star known to have more than four transiting planets.Luque et al. (2023) discovered six transiting planets with sizes ranging from 1.9 to 2.9 R ⊕ and orbital periods ranging from 9.1 to 55 days.Starting with the innermost planet, the period ratios between adjacent planetary orbits are nearly equal to 3:2, 3:2, 3:2, 4:3, and 4:3, suggesting that the system forms a resonant chain.If true, HD 110067 would be one of only three known systems with six or more planets trapped in a resonant chain, the other two being TRAPPIST-1 (Gillon et al. 2017;Luger et al. 2017) and TOI-1136 Table 1.Properties of the HD 110067 planetary system, as reported in Luque et al. (2023).Planets c, e, and g have only 3 σ upper limits on their masses, and all six eccentricities are unconstrained.Our N -body simulations were initialized based on the observed posteriors on the orbital parameters of the six planets (see Section 2 for more details).

Planet
Period a Periods determined based on two transits and the prior that the system is in a specific resonant chain configuration.(Dai et al. 2023).Under the assumption that the system forms a resonant chain, Luque et al. (2023) successfully scheduled follow-up observations of the transits of the two outer planets before their periods had been well established.However, because the number of observed transits is still relatively small, and because TTVs have not yet been detected, the case for a resonant chain in HD 110067 has not been established as securely as it was for TRAPPIST-1 and TOI-1136.Due to the potential importance of the HD 110067 planetary system, dynamical characterization is warranted (see, e.g., MacDonald et al. 2016;Mills et al. 2016;Tamayo et al. 2017;Siegel & Fabrycky 2021;Mac-Donald et al. 2022;Quinn & MacDonald 2023 for dynamical analyzes of other resonant chain systems).This Letter describes long-term N -body simulations that were undertaken to explore the dynamics of HD 110067.By assuming that the system is dynamically stable on timescales comparable to the age of the system, we wanted to know how confidently we may conclude that the system is a resonant chain.Since only three of the planets' masses have been measured (using the Doppler technique), and none of the eccentricities have been measured, we also wanted to know whether dynamical arguments can be used to place meaningful constraints on the masses and eccentricities.

SIMULATION SETUPS
We have carried out N -body simulations of HD 110067-like systems initialized in two different ways: with parameters that are consistent with the available observations and no other constraints; and with parameters that emerge from a simple model of disk-driven migration into the observed configuration.Below, we describe the setup of our migration simulations.
The theory of resonant capture in two-planet systems is fairly well understood analytically (Henrard 1982;Borderies & Goldreich 1984;Mustill & Wyatt 2011;Batygin et al. 2015), but the capture of highermultiplicity systems into resonance is more complex and is typically modeled numerically.In this work, we adopt the migration model of Tamayo et al. (2017), which allows us to generate many physically plausible resonant chains that are consistent with the observed properties of HD 110067.By initializing planets with random orbital orientations and sampling different choices of parameters governing the timescale and turbulence of migration, this procedure produces resonant chain systems that span a range of plausible eccentricities, deviations from exact resonant period ratios, and libration amplitudes.Similar approaches have been adopted in subsequent works (e.g., MacDonald & Dawson 2018;Siegel & Fabrycky 2021;MacDonald et al. 2022).Below, we describe our implementation of this migration model; more background on the method can be found in the works cited above.
At the outset of our simulations, the five outer planets were assigned periods that placed them 2 % wide of the exact MMRs they are found near today (i.e., initial periods were chosen so that ∆ = 0.02 for each adjacent planet pair).Each planet was initialized on a circular orbit, with an initial inclination and a mass drawn from independent normal distributions based on the posteriors reported in Table 1.We set the stellar mass to M * = 0.8 M ⊙ and drew the initial orbital orientations (i.e., true longitudes and longitudes of the ascending node) from a uniform distribution spanning [0, 2π].
To capture the planets into a sequence of MMRs, we damped the semi-major axis of the outermost planet exponentially, with a damping timescale of τ a = 5 × 10 6 P 1 = 1.2 × 10 5 years. 1 The eccentricity of each planet was exponentially damped with a timescale τ e = τ a /K, where K is a was drawn from a log-uniform distribution from 10 to 10 3 .We integrated each system for one τ a timescale and then adiabatically removed the 1 τa was chosen by manually increasing its value until the migration process was slow enough to successfully capture the six planets into a 3:2, 3:2, 3:2, 4:3, 4:3 resonance chain with a ∼ 50 % success rate.eccentricity and semi-major axis damping over 5 × τ e .
To produce systems with a range of libration amplitudes (as expected due to turbulent perturbations during migration; e.g., Rein & Papaloizou 2009), we then applied a parameterized kick to each of the planets' eccentricities.Specifically, we drew a random number η from a log-uniform distribution between 10 −3 and 1, and then either increased or decreased (with equal probability) each planet's eccentricity by a factor of η.
Finally, configurations that did not end up being captured into a 3:2, 3:2, 3:2, 4:3, 4:3 resonance chain were discarded.Configurations that successfully locked into the desired resonant chain were rescaled so that P 1 = 9.11678 days.Fig. 1 shows an illustrative example of our migration procedure.Resonant capture in our simulations proceeds outside-in.First, planets f and g are captured into the 4:3 MMR, followed by planets e and f, and then the inner three planet pairs are sequentially captured into 3:2 MMRs.Over the course of this process, each planet acquires a new equilibrium eccentricity, which can be influenced by the capture of subsequent planets.After ∼ 3 × 10 6 P 1 , the system is trapped in the resonant chain configuration, with librating twoand three-body resonant angles.
We then carried out long-term N -body integrations using the WHFast integrator (Rein & Tamayo 2015) from the open-source REBOUND package (Rein & Liu 2012) with a time step of P 1 /20, where P 1 is the initial period of the planet b.Integrations were stopped when the Hill radii [r H, i = a i (m i /M * ) 1/3 ] of two planets intersected, or after 10 9 P 1 ≈ 25 Myr, whichever came earlier.Damping forces were implemented using the modify orbit forces functionality from REBOUNDx (Tamayo et al. 2020).

LONG-TERM STABILITY
To simulate systems that are consistent with observations but are not necessarily in resonance, we integrated 350 systems initialized by drawing parameters directly from the observational posteriors.In detail, we drew the initial periods from independent normal distributions with mean and standard deviation taken from Table 1.2Similarly, we drew the masses of the planets with radial-velocity (RV) measured masses (planets b, d, and e) from independent normal distributions, and we drew masses for the remaining three planets independently from [0, m max ], where m max is the 3σ RV mass upper limit from Luque et al. (2023).Inclinations were also drawn independently from the posteriors (taking into account the discrete degeneracy in inclination inherent to transit observations) and eccentricities for each planet were drawn uniformly from [0, 0.05], motivated by the eccentricity distribution that has been inferred for nonresonant multiplanet systems (Van Eylen & Albrecht 2015; Hadden & Lithwick 2017; the initial eccentricities of our resonant chain systems span a similar range).Each longitude of pericenter and longitude of the ascending node was drawn randomly from [0, 2π], and true longitudes were determined from the reported transit times.Mixed systems, in which the inner planets are resonant but the outer planets are not, were constructed by carrying out the migration process described above with a Figure 2. Cumulative fraction of 350 systems that destabilized within 25 Myr.99 % of the nonresonant systems with properties consistent with the observations were unstable, as opposed to 11 % of systems initialized in a six-planet resonant chain (via a disk migration simulation).Systems in which only a subset of the planets are captured into a resonant chain were also mostly unstable.
system of fewer planets, then adding the outer nonresonant planets to the system by drawing their properties from the observational posteriors.
Figure 2 shows the cumulative fraction of systems that proved to be unstable as a function of time.Systems with parameters drawn from the posteriors of Luque et al. (2023), subject to no other constraints, were overwhelmingly unstable, with 345 of the 350 systems (99 %) destabilizing within 25 Myr (∼ 0.3 % of the system's age).In contrast, only 38 of the 350 systems (11 %) that underwent convergent migration into a six-planet resonance chain destabilized within 25 Myr.Systems in which the inner five planets formed a resonant chain but the outermost planet was decoupled from the resonant chain were also largely unstable, with 65 % of such systems going unstable within 25 Myr.Decoupling more of the outermost planets from the resonant chain increased the fraction of systems that destabilized.
As a further experiment, we carried out N -body simulations in which a subset of the outer planets formed a resonant chain, and the inner planets were decoupled.This was achieved by carrying out migration simulations with a reduced number of planets, then adding the inner planets to the system by drawing their properties from the posteriors.With the innermost one, two, or three planets initialized directly from the observed posteriors, 41 %, 58 %, and 75 % of systems destabilized within 25 Myr, respectively.We consider these results to be strong evidence that the six planets of HD 110067 form a six-planet resonant chain.

CONSTRAINTS ON ORBITAL ARCHITECTURE
Next, we explored the possibility of constraining the properties of the planets in HD 110067 by requiring that the six-planet resonant chain systems are long-term stable.The small sample size of the ensemble of simulations presented above was a barrier: there were only 38 unstable systems out of the total of 350 resonant chain systems, limiting our ability to find trends among the unstable systems.To expand the sample, we carried out 1050 additional N -body simulations of six-planet resonant chain systems (of which 111 destabilized), increasing the total number of resonant chain simulations to 1400.
First, we considered the influence of the planets' orbital eccentricities on the long-term stability of our synthetic HD 110067-like systems.Unsurprisingly, systems with larger eccentricities were more likely to destabilize, but the initial eccentricities of the six planets affected the stability of the systems to varying degrees.For our estimate of the upper limit on each planet's eccentricity e i , we found the value of e max for which more than 33 % systems with e i > e max destabilized within 25 Myr.We report our upper limits on planet eccentricities in Table 2.The weakest constraint is for planet c, probably because planet c is the only planet in the system in a 3:2 MMR (as opposed to 4:3) with both its inner and outer neighbors.The tightest constraint is on the eccentricity of planet g, probably because of the lack of the stabilizing influence of an outer resonant planet and the fairly tight (4:3) MMR with planet f. Figure 3 shows the cumulative unstable fraction of resonant systems that violate the eccentricity constraints given in Table 2.The sharp upturn and rising trend in the unstable fraction over 25 Myr makes it unlikely that many of these configurations would survive for the full ∼ 8 Gyr age of HD 110067.
Next, we turned to the influence of the planets' masses on long-term stability.Planets a, c, and e have RVmeasured masses, but planets c, e, and g only have RVbased upper limits (see Table 2).We did not find the systems for which m c , m e , or m g approached the RV mass upper limits to be especially unstable, indicating that the RV-based upper limits are more strict than the upper limits set by the requirement of long-term stability (although we did see an uptick in instability when m c , m e , or m g approached about 20 M ⊕ ).On the other hand, systems in which a very low mass was assigned to planet c, e, or g were significantly more likely to destabilize on long timescales.This is probably because The dashed and dotted curves, respectively, show the cumulative unstable fraction of the 278 systems that contain a low-mass planet and the 206 systems that contain a higheccentricity orbit (i.e., systems with masses or eccentricities outside of the limits reported in Table 2).Systems that contain a low-mass or high-eccentricity planet are significantly more likely to destabilize on long timescales, allowing us to place approximate lower limits on planet masses and upper limits on orbital eccentricities for the HD 110067 system.low-mass planets are more strongly perturbed by the gravitational interactions with the other planets in the system, making them more vulnerable to instability (see also Section 6).
Lower limits on the masses of planets c, e, and g, determined analogously to the eccentricity upper limits, are given in Table 2.As was the case for the orbital eccentricities, the mass constraints are most stringent for planet g and least stringent for planet c.The cumulative fraction of unstable resonant chain systems that violate these mass lower limits increases rapidly at late times (see Fig. 3).Systems that violate any of the eccentricity/mass constraints reported in Table 2 make up 86 % of the unstable resonant chain simulations, whereas they account for 30 % of all resonant chain systems in our sample.There is only moderate overlap between the high-eccentricity and low-mass samples; 67 systems belong to both samples.
In addition to the properties of the planets, we attempted to constrain the values of K and η, inspired by the observation that the unstable systems were initialized with preferentially low K values and preferentially high η values when compared with the stable systems.We found that requiring long-term stability leads to weak constraints on K and η.Specifically, by requiring that fewer than 33 % of systems should destabilize within 10 9 P 1 , we found that K ≳ 13 and η ≲ 0.98 (in other words, there are stable resonant chain configurations initialized with K ≈ 10 -1000 and η ≈ 0.001 -1.00).These relatively weak constraints prevent us from drawing strong conclusions about the convergent migration of HD 110067.

OBSERVATIONAL COMPARISONS AND FUTURE PROSPECTS
Unlike two-body resonant angles, which depend on the planets' longitudes of pericenter, three-body resonant angles depend only on the planets' mean longitudes and can therefore be determined directly from transit timings (e.g., ϕ bcd = 2λ b − 5λ c + 3λ d ). Figure 4 shows the distribution of the three-body resonant angles across the stable resonant chain systems at the beginning of our Nbody integrations.Because the resonant state of the system depends on how migration proceeds, different initial conditions can result in different libration centers and amplitudes.Assuming low orbital eccentricities, the observed transit times of the six planets (which have negligible uncertainty) can be used to determine the value of the mean longitudes, and therefore the three-body resonant angles, at one instant in time.The distributions of ϕ cde , ϕ def , and ϕ ef g are consistent with the observed values, although the observed value of ϕ bcd is higher than 95 % of the ϕ bcd values in our simulated systems (a "tension" of 1.6 σ; see Section 6).Additional observations of HD 110067 are required to determine the true three-body libration centers and amplitudes.
Shifting our attention to the future, we expect that TTVs will be detected and will provide more precise constraints on the planets' masses and orbits.The TTVs of a planet pair that lie near a two-body j:j − k MMR are expected to oscillate over the superperiod P s = 1 j Pout ∆ (Lithwick et al. 2012).Due to the close proximity of planets to first-order MMRs in HD 110067 (∆ ∼ 10 −4 ), superperiod oscillations would occur over long timescales (≳ 15,000 days), much longer than the  2022), for HD 110067, we expect P ℓ ≈ 2,000 -6,000 days, depending on the planet pair.
The amplitude of the TTV signals depends on the planets' proximity to resonance, masses, eccentricities, and orbital orientations.As a result, the predicted TTVs vary widely across our synthetic resonant chain systems.Commonly, the predicted TTVs show ∼ 10 hr oscillations over thousands of days, in agreement with our analytical estimate of P ℓ .The median peak-to-peak amplitude and 1 σ spread of the TTVs in our simulated systems is 7 +15 −6.7 hr.To illustrate the diversity of potential TTVs for HD 110067, Fig. 5 shows the TTVs exhibited by three simulated resonant chain systems in our sample.Although there exist a few systems in our sample with negligible short-term TTVs, the typical TTV signal is large enough to be easily detectable with a several-year observational baseline.Through such an investigation, we have concluded that HD 110067 is very likely to be trapped in a six-planet resonant chain.Assuming a resonant configuration, and requiring that the system be stable over timescales com-parable to its age, we also derived constraints on orbital eccentricities and planet masses.
HD 110067 is notable among resonant chain systems in that when planetary parameters are drawn randomly from the distributions allowed by the available data, the resulting system is unstable at least 99 % of the time.Previously, Hadden & Lithwick (2017) noted that > 10 % of N -body simulations of resonant chain systems Kepler-60 and Kepler-223 go unstable over 10 6 P 1 .Similarly, Tamayo et al. (2017) found that ∼ 40 % systems initialized based on the observed posteriors of TRAPPIST-1 destabilized within 10 9 P 1 , whereas systems initialized via disk migration simulations were largely stable over 10 9 P 1 .HD 110067 is an even more extreme example of this trend.Part of the explanation is that HD 110067 was discovered recently and the current observational constraints are relatively weak.The dynamical fragility of the system is likely also related to the system's high multiplicity and compactness, with all adjacent planet period ratios Pout Pin ≲ 1.5.The proximity of adjacent planets in the system to first-order MMRs (∆ ≈ 10 −4 ) is also a factor, as instability times are known to shorten by several orders of magnitude in the vicinity of such MMRs (e.g., Obertas et al. 2017;Petit et al. 2020;Lammers et al. 2024).
Before carrying out our N -body simulations, we did not anticipate being able to place lower limits on any of the planet masses.Intuitively, one might think that since lowering a planet's mass weakens its mutual gravitational interactions, low planet masses promote longterm stability.However, this may not be true if there are other planets in the system with known masses, as is the case for HD 110067.Low-mass planets are more strongly perturbed by the other planets, potentially disturbing the resonant chain configuration.However, the exact dynamical mechanism for the instability is somewhat unclear.Theoretical works have suggested that resonant chain systems are destabilized by a resonance between a fast libration frequency and a slow difference between two synodic frequencies (Pichierri & Morbidelli 2020;Goldberg et al. 2022).It is unclear why, in this dynamical picture, the inclusion of a low-mass planet would cause systems to destabilize; possibly, overstable librations play a role (Goldreich & Schlichting 2014).
The observed value of the three-body angle ϕ bcd = 2λ b − 5λ c + 3λ d ≈ 224 • is larger than 95 % of the ϕ bcd values in our simulated six-planet resonant chain systems, which cluster with fairly low amplitudes (∼ 20 • ) about the 180 • equilibrium.Although the statistical significance of this "tension" is modest, we wondered about the possible reasons, given that previous works were able to reproduce the observed three-body resonant angles of other resonant chains using similar methods (e.g., Mills et al. 2016;Tamayo et al. 2017).
Through experimentation, we found that the discrepancy is not significantly reduced by modifying the range of K and η values in our migration simulations.The large observed value of ϕ bcd can be explained if ϕ bcd is circulating, rather than librating.However, if the planets b, c, and d are not trapped in a three-body resonance, our N -body simulations indicate that the system would likely be unstable.With a shortest period of P ≈ 9 days, it is also unlikely that tides from the host star are affecting the dynamics of the system.Could the discrepancy be caused by an undetected inner planet in the HD 110067 system?To explore this possibility, we carried out migration simulations with a 1 -10 M ⊕ innermost planet in a first-/second-order MMR with planet b.We found that including an additional inner planet in the system does widen the distribution of ϕ bcd values, and relaxes the tension with the observed value from 1.6 σ to about 1.0 σ.Obviously, this is not yet a compelling hypothesis, but it should be kept in mind as additional data are collected.Given the system's observational accessibility and rich dynamical structure, HD 110067 promises to be a popular system for observers and theorists alike.

Figure 1 .
Figure1.Evolution of adjacent planet period ratios, orbital eccentricities, and three-body resonant angles during an example convergent migration simulation (in this case, K = 114).Adjacent planet pairs are captured into two-body MMRs outside-in, beginning with planets f-g and ending with planets b-c.Each planet obtains a non-zero eccentricity through this process, and three-body resonant angles librate with small amplitudes (two-body resonant angles also librate).Before carrying out our long-term N -body simulations, we adiabatically remove the eccentricity and semi-major axis damping and provide a random kick to the planets' eccentricities.

Figure 3 .
Figure3.The solid curve shows the cumulative unstable fraction of the 1400 HD 110067-like resonant chain systems.The dashed and dotted curves, respectively, show the cumulative unstable fraction of the 278 systems that contain a low-mass planet and the 206 systems that contain a higheccentricity orbit (i.e., systems with masses or eccentricities outside of the limits reported in Table2).Systems that contain a low-mass or high-eccentricity planet are significantly more likely to destabilize on long timescales, allowing us to place approximate lower limits on planet masses and upper limits on orbital eccentricities for the HD 110067 system.

Figure 4 .(
Figure 4. Distributions of three-body resonant angle values in our 1251 stable HD 110067-like resonant chain systems at the outset of the N -body simulations.Horizontal dashed lines mark the observed values, based on measured transit times.For all of the three-body resonant angles except ϕ bcd , the observed value falls comfortably within the range of resonant angles in our simulated resonant systems.

6.
DISCUSSION AND CONCLUSION Thanks to the gigayear ages of typical observed exoplanet systems, investigations of long-term dynamical stability often help to characterize exoplanetary systems beyond what current observations can provide.

Figure 5 .
Figure5.Transit timing variations (TTVs) over 3000 days for three HD 110067-like resonant chain systems.The predicted TTVs typically evolve over timescales of thousands of days.The predicted amplitudes vary widely across our simulated systems, with a median peak-to-peak amplitude of 7 +15 −6.7 hr.

Table 2 .
Constraints on the orbital eccentricities and planet masses of HD 110067 from the requirement that < 33 % of six-planet resonant chains destabilize within 25 Myr.The masses of planets b, d, and f have been measured with RV data (see Table1).