Spin Dynamics of Planets in Resonant Chains

About a dozen exoplanetary systems have been discovered with three or more planets participating in a sequence of mean-motion resonances. The unique and complex architectures of these so-called “resonant chains” motivate efforts to characterize their planets holistically. In this work, we perform a comprehensive exploration of the spin-axis dynamics of planets in resonant chains. Planetary spin states are closely linked with atmospheric dynamics and habitability and are thus especially relevant to resonant chains like TRAPPIST-1, which hosts several temperate planets. Considering a set of observed resonant chains, we calculate the equilibrium states of the planetary axial tilts (“obliquities”). We show that high-obliquity states exist for ∼60% of planets in our sample, and many of these states can be stable in the presence of tidal dissipation. Using case studies of two observed systems (Kepler-223 and TOI-1136), we demonstrate how these high-obliquity states could have been attained during the initial epoch of disk-driven orbital migration that established the resonant orbital architectures. We show that the TRAPPIST-1 planets most likely have zero obliquities, with the possible exception of planet d. Overall, our results highlight that both the orbital and spin states of resonant chains are valuable relics of the early stages of planet formation and evolution.


INTRODUCTION
One of the most intriguing types of exoplanetary systems are resonant chains such as TRAPPIST-1 (Gillon et al. 2017).Resonant chains contain three or more planets participating in a sequence of mean-motion resonances (MMRs), in which the planetary orbital periods are related by ratios of simple integers.Such configurations likely arise as a consequence of convergent migration driven by planet-disk interactions (e.g. Lee & Peale 2002;Cresswell & Nelson 2006;Terquem & Papaloizou 2007).If disk migration is ubiquitous, one might expect that resonant chains would be very common.However, this is not the case; resonant chains (and MMRs in general) are rare throughout the exoplanetary sample (Lissauer et al. 2011;Fabrycky et al. 2014).The paucity of resonances has been a topic of extensive research (e.g.Adams et al. 2008;Goldreich & Schlichting 2014;Batygin & Adams 2017;Izidoro et al. 2017).In this paper, we will not explore why resonances are rare, but rather we will focus on developing a deeper understanding of the small number (∼ 10) of known resonant chains.These systems are valuable relics of the planet formation and sarah.millholland@mit.edumigration process, and they also may host potentially habitable planets.It is thus beneficial to characterize them thoroughly and holistically.
One underexplored aspect of resonant chains is the dynamical behavior of the planetary spin vectors.Planetary rotation rates and axial tilts ("obliquities") are closely linked with fundamental geologic and atmospheric features, including climate stability and habitability (e.g.Atobe et al. 2004;Spiegel et al. 2009;Kilic et al. 2017;Shan & Li 2018;Shakespeare & Steffen 2023).For close-in planets, spin dynamics can even affect a planet's orbital and bulk physical properties by way of the connection between spin and tidal dynamics.Planets with non-zero obliquities experience a dissipative tidal torque that generates heat in the planetary interior at the expense of orbital energy.These "obliquity tides" can drive orbital decay and planetary radius inflation, which have been theorized to explain a range of exoplanetary observations, ranging from features in the period ratio and radius ratio distributions (Millholland & Laughlin 2019;Millholland 2019;Goldberg & Batygin 2021;Li 2021;Hong et al. 2021;Su & Lai 2022a) to the existence of ultra-short period planets (Millholland & Spalding 2020) and puffy sub-Saturns (Millholland et al. 2020).
Although the impacts of planetary spin states are wide-ranging, we are currently unable to directly observe them for the majority of planets, with the exceptions being some wide-orbiting giant planets and brown dwarfs (Bryan et al. 2020(Bryan et al. , 2021)).Obliquities of close-in planets may soon be detected in high-precision photometry through the distortion of the transit light curve created by planetary oblateness (Carter & Winn 2010a,b;Berardo & de Wit 2022).For now though, we must rely on theory to explore planetary spin states and their past histories.This process, in turn, requires theoretical investigations of the orbital architectures, given the tight coupling between the spin and orbital dynamics (e.g.Correia & Delisle 2019).Compared to typical planetary systems which leave little trace of their formation history, the orbital histories of resonant chains are more well-understood since they can only arise through convergent orbital migration (e.g.Terquem & Papaloizou 2007).The equilibrium positions and libration amplitudes of the critical resonant angles may even encode clues to the disk conditions (Adams et al. 2008;Rein & Papaloizou 2009;Delisle 2017) (though there are some caveats; see Jensen & Millholland 2022).
The orbital migration process can strongly affect the planetary spin states (e.g.Vokrouhlický & Nesvorný 2015;Brasser & Lee 2015;Millholland & Laughlin 2019).As the semi-major axes change, the rates of orbital precession induced by planet-planet interactions change up to orders of magnitude in the case of long-range migration.The evolving orbital precession rates can cross through commensurabilities with the rates of precession of the planetary spin axes.Specifically, commensurabilities between the precession rates of the longitude of ascending node and the spin vector are called "secular spin-orbit resonances" (e.g.Peale 1969;Ward 1973;Ward & Hamilton 2004).When they are encountered, they can force planetary spin axes to tilt to high obliquities, which the planets can then maintain for billions of years, even in the presence of tidal dissipation.Secular spin-orbit resonances are common in the Solar System.Most notably, the mechanism is thought to be responsible for Saturn's 27 • obliquity (Ward & Hamilton 2004;Hamilton & Ward 2004;Wisdom et al. 2022;Saillenfest et al. 2021;Wisdom et al. 2022).They are also thought to be common in extrasolar systems (Millholland & Laughlin 2019;Su & Lai 2020;Li 2021;Su & Lai 2022b).
High-multiplicity systems are particularly susceptible to secular spin-orbit resonances.A system of N p planets contains N p − 1 nodal precession frequencies, and the total nodal evolution is a superposition of the various modes (Murray & Dermott 1999).In a secular spin-orbit resonance, the planet's spin-axis precession frequency is resonant with just one of the nodal frequency modes.Many of the observed resonant chains have high multiplicities (e.g.≳ 5 observed planets), meaning that there is a multitude of nodal frequency modes that the planetary spin axes could be resonant with.If the frequencies are closely-separated, it can lead to chaotic spin dynamics (Shan & Li 2018;Saillenfest et al. 2019), as it does for the terrestrial planets in the Solar System (Laskar & Robutel 1993;Touma & Wisdom 1993).
Altogether, the combination of a history of orbital migration and the existence of many orbital frequency modes suggests that secular spin-orbit resonances could be a common occurrence in resonant chains.If this is the case, high obliquities could be prevalent.However, this conjecture has not yet been tested.In this work, we conduct a systematic analysis of the spin dynamics of planets in observed resonant chains.We will approach the problem in two parts: (1) We will first identify the locations of the spin equilibria and numerically assess their stability, remaining agnostic to how these states may have been reached in the first place.(2) Second, we will explore the spin dynamical history and consider how disk migration may have simultaneously created the resonant orbital architecture and the present-day spin states.These steps will allow us to answer whether high obliquity states are common for planets in resonant chains.
This paper is organized as follows.We begin by reviewing the dynamics of secular spin-orbit resonances (Section 2).Considering our sample of resonant chains (Section 3), we first identify the spin equilibria for each planet and assess their stability (Section 4).We then explore the formation and migration history of two case study systems (Section 5).We show that Kepler-223 is a particularly strong candidate for stable high obliquities.We discuss further implications in Section 6 and conclude in Section 7.

SPIN DYNAMICS AND CASSINI STATES
Planetary spin vectors are subject to complex dynamical evolution as a result of the combined effects of spin-axis precession and orbital precession.The torque induced by the host star on the rotationally-flattened planet causes the planet's spin vector to precess around its orbital axis with a period equal to T α = 2π/(α cos ϵ), where ϵ is the obliquity and α is the spin-axis precession constant.In the absence of satellites orbiting the planet, α is given by (Neron de Surgy & Laskar 1997) where M ⋆ is the stellar mass, M p is the planet mass, R p is the planet radius, a is the semi-major axis, k 2 is the Love number, C is the planet's moment of inertia normalized by M p R p 2 , ω is the spin rate, and e is the eccentricity.
However, within its precession around the orbital axis, the spin vector is chasing a "moving target".The orbit plane itself undergoes precession due to interactions with other planets, the oblate host star, and any other gravitational sources that cause deviations from a 1/r potential.In the case of a single perturber, such as one companion planet, the orbital precession is uniform, meaning that the orbital inclination I and nodal precession frequency 1 g = Ω < 0 (where Ω is the longitude of ascending node) are constant.The orbital precession period is T g = 2π/|g|.
Within a uniformly precessing orbit frame, the equilibrium configurations of the spin vector are known as "Cassini states".Cassini state dynamics have been reviewed in many previous works (e.g.Colombo 1966;Peale 1969Peale , 1974;;Ward 1975;Ward & Hamilton 2004;Correia 2015;Saillenfest et al. 2019;Su & Lai 2020, 2022a,b).Here we recap the aspects that are relevant to our discussion.Cassini states are configurations in which the planetary spin axis, ŝ, and the unit orbit normal vector, n, precess at the same rate about the axis of the total system angular momentum vector, k.In the absence of dissipation, the three vectors are coplanar, whereas a Cassini state with dissipation has ŝ slightly shifted out of the plane defined by n and k (Su & Lai 2022a).For a given orbital inclination I, the obliquity ϵ of a Cassini state satisfies the relation g sin(ϵ − I) + α cos ϵ sin ϵ = 0. (2) Equation 2 has either two or four roots, depending on the value of the frequency ratio |g|/α and its comparison to a critical ratio, To avoid confusion, we note that many papers use g to denote apsidal precession frequencies and s nodal precession frequencies (e.g.Murray & Dermott 1999).However, literature on Cassini states often uses g to represent the nodal frequencies (e.g. Ward & Hamilton 2004), so we aim to be consistent with that.to be overwhelmed by the tidal damping torque (Millholland et al. 2020;Su & Lai 2022a).Figure 1 shows the obliquities of the four Cassini states as a function of |g|/α.Cassini states are strictly only defined in the case of uniform precession.However, in most configurations, the precession is non-uniform, as it is driven by multiple perturbers (e.g. more than one companion planet).In this case, the inclination/node evolution may be described as the superposition of various modes with multiple frequencies {g i } and amplitudes {I i } (Murray & Dermott 1999).The spin-vector may participate in what we define as "quasi-Cassini states", which are, roughly speaking, Cassini states with g in equation 2 replaced by one of the g i modes.In this paper, we will often use the phrase "Cassini state" even when we are actually referring to "quasi-Cassini states".We note that more complicated equilibria involving a mixture of frequency modes also exist (Saillenfest et al. 2019;Su & Lai 2022b).
There are multiple physical mechanisms that cause planets to enter Cassini states.One process that is relevant to close-in planets is tidal dissipation, which acts on the spin vectors and causes planets to dissipate into the equilibria (Millholland & Spalding 2020).Planets with |g|/α > (|g|/α) crit are dissipatively forced into state 2, whereas planets with |g|/α < (|g|/α) crit will dissipate into state 1 or 2, depending on the initial obliquity.If the initial obliquity is greater than the obliquity of state 2, there is a chance it will enter into state 2, depending on the dissipation timescale.Otherwise, the planet will end up in state 1.
A second process by which both close-in (tidally active) and more distant (tidally inactive) planets can enter Cassini states is through secular spin-orbit resonance encounters.This is the scenario that is thought to be relevant for the Solar System giant planets (e.g. Ward & Hamilton 2004).This process requires that the ratio of precession periods T α /T g = |g|/(α cos ϵ) crosses through unity from above.That is, either α, |g|, or ϵ must evolve such that T α /T g moves from values > 1 to < 1 as a response to some form of physical evolution of the system.If the evolution is slow enough, it can result in obliquity excitation into Cassini state 2.Many physical processes can cause the resonant evolution of T α /T g .One notable example is orbital migration, which modulates both |g| and α.
Before proceeding, we also note that Valente & Correia (2022) showed that planets participating in certain "traditional" spin-orbit resonances (involving commensurabilities between the rotation period and orbital period) can experience an unexpected tidal excitation of the obliquity, which can be maintained over the lifetime of the system.This mechanism can occur even without perturbing planets and represents yet another mode of obliquity excitation that may be relevant for the shortperiod planets considered in this study.
The parameters in Table 1 are the latest best-fit constraints from the papers indicated in the table caption.However, the best-fit parameters come from summary statistics applied to the marginalized posterior distributions and often do not correspond to resonant trajectories.For this reason, and in order to accurately represent the parameter uncertainties, we utilize samples directly from posterior distributions in the majority of this work.We use posterior distribution data from the following papers: Agol et al. (2021) for TRAPPIST-1,

SPIN EQUILIBRIA
Our first objective is to identify the set of spin equilibria of the planets in our sample.The resonant chains must have underwent some degree of migration in order to enter their observed configurations.Orbital migration causes both the |g| and α frequencies to evolve, so it is capable of generating a resonant evolution of T α /T g and exciting a high obliquity in Cassini state 2. We note that the concept of a "high obliquity" is somewhat arbitrary, so for the sake of definiteness, we define a high obliquity to be ≳ 20 • .In this section, we aim to simply identify the possible high obliquity states Table 1.Resonant chain systems and their basic parameters.We show the nominal periods, masses, and radii based on the references below.Note that we do not use these nominal values in our analysis but rather we sample from the posterior distributions as described in Section 4.

P [days]
Mp while remaining agnostic about how they might have originated.(We will investigate their origins later in Section 5 using migration simulations.)The exploration in this section will be split into two parts.First, we will use a frequency analysis to identify the possible equilibria.Second, we will use spin-orbit simulations to probe the stability of the equilibria in the presence of tidal dissipation.

Frequency analysis
In order to identify the equilibria of the spin-orbit resonances, we first need to calculate the spin-axis precession frequencies and orbital frequencies, which we will then use to calculate the associated Cassini state 2 values.Here we develop a method of calculating the spectrum of orbit nodal precession frequencies, {g i }, in a given system.It is not possible to use Laplace-Lagrange secular theory due to the multitude of MMRs in the res-onant chains.It is most straightforward to obtain the frequency spectrum through numerical methods.The nodal precession frequencies are primarily a function of the planet masses and semi-major axes, with a weaker dependence on the eccentricities and inclinations.Given the spread in these parameters within a given system's posterior distribution, each g i frequency has some uncertainty associated with it.To account for this, we will sample parameter vectors from the posterior distribution and average the results.Our analysis for each system proceeds as follows: 1. Select initial conditions: First, we pick a set of 10 random vectors of system parameters from the posterior distribution (described in Section 3), which constitute separate sets of initial conditions.The parameters are the masses, periods, eccentricities, and transit times or mean anomalies.We randomize all planets' inclinations according to i ∼ Rayleigh(0.05 • ).Non-zero inclinations are required for nodal precession to occur, but they are chosen to be small since the resonant chains are nearly coplanar (e.g.Agol et al. 2021).This assumption is also conservative since Cassini states are more stable to tidal dissipation with higher inclinations.2. Perform orbital integrations: For each of the 10 sets of initial conditions, we perform a short 10 5 year N -body integration using a Bulirsch-Stoer integrator that will be described in Section 4.2.Here we do not include tidal or spin dynamics, but we do include general relativistic precession (Mardling & Lin 2002).

Calculate nodal frequency spectra for each planet:
For each of the 10 N -body integrations, we use the resulting orbital elements to calculate the evolution of each planet's orbit normal vector, n, as a function of time.We then obtain a frequency spectrum of the xcomponent of n for each planet.(The result is not sensitive to which component of n is chosen.)We use a Lomb-Scargle periodogram, although we also tested a fast Fourier transform and found similar results.The peaks in the spectrum indicate the nodal frequencies.It is important to note that the full set of N p − 1 nodal frequencies is common to the system as a whole and does not vary from planet to planet.However, the strength of each mode does vary between planets, so it is necessary to calculate a spectrum for each planet in order to observe these differences.4. Compute average spectra for each planet: Having calculated a frequency spectrum for each planet and for each of the 10 sets of initial conditions, we now compute an "average" spectrum for each planet across the range of initial conditions.A modified average is necessary because the frequencies are shifted slightly between each of the 10 different iterations, with the shift being predominantly driven by the total mass in planets, M p,tot = Σ Np i=1 M p,i , where N p is the number of planets in the system.To account for this, we first scale the frequencies by M p,tot /M p,tot , where M p,tot is the average across the 10 sets of initial conditions.We then calculate the mean spectrum for each planet.5. Identify peak frequencies: For each planet's average spectrum, we locate the dominant frequencies using a peak finding algorithm.We define peaks as local maxima > 10% of the height of the tallest peak and separated by more than ∼ 0.1 dex in frequency space.We keep track of both the peak frequencies for each planet, as well as the collection of N p − 1 unique frequencies across the system as a whole.6. Calculate Cassini state obliquities: With the peak nodal frequencies identified for each planet, we can calculate the nominal equilibrium obliquities associated with each frequency.We assume the planet's spin rate is at equilibrium (at which dω/dt = 0).In the limit e ≈ 0, this depends on the obliquity as (Leconte et al. 2010) where n = 2π/P is the mean-motion.Combining equations 1, 2, and 4, we numerically solve for the Cassini state 2 obliquities that satisfy equation 2. We define α syn to be the spin-axis precession constant in the case of synchronous rotation (ω = n), and we calculate the average α syn across all the selected posterior parameter vectors.Here and elsewhere in the paper, we fix the planets' Love numbers to k 2 = 0.3 and the normalized moments of inertia to C = 0.3.Such values are reasonable for super-Earths/sub-Neptunes (Kramm et al. 2011), and although other values are possible, it is simplest to make a uniform choice with k 2 = C given the uncertainties in these parameters.We also assume that I = 0.1 • , which is arbitrary but commensurate with our assignment of inclinations.It is helpful to note that, in the limit and the Cassini state 2 obliquities can be computed as However, we use a numerical solution of the full equation 2 for all calculations.

Results
We now report the results of the frequency analysis.Table 2 displays the N p − 1 nodal frequencies, {g i }, in each resonant chain system.These were found by taking the set of unique frequencies among among all peak frequencies identified in Step 5 above.In some cases, the collection of peak frequencies across all planets yields more than N p − 1 unique frequencies.This occurs either when the average spectra do not perfectly account for the frequency shifts or when the system is chaotic.In these cases, we simply identify the N p − 1 most wellseparated frequencies.
As an example, Figure 3 displays the nodal frequency spectra for the planets in the TRAPPIST-1 system.The results show six (= N p − 1) clear unique peak frequencies, which form the {g i } set.The frequency modes have various strengths for each planet, such that only a subset of the six are significant in each case.For instance, only two of the six g i frequencies are strong for the outermost planet.The figure also includes the estimates of α syn , the spin-axis precession constant evaluated at ω = n.We consider the case with k 2 /C = 1 to be Figure 4 shows the nodal frequency spectra for all systems in our sample.There are many planets with overlap between the nodal frequencies and the α syn ranges.These include Kepler-80 d, e, b, and c, all of the Kepler-223 planets, TOI-1136 d and f, Kepler-60 b and c, GJ 876 d, K2-138 c through g, and TOI-178 c through g.For many of the planets, there are multiple frequencies that overlap with the α syn range, indicating multiple resonant equilibria.Planets with nodal frequencies less than the α syn range can also have significant Cassini state 2 obliquities, but they will be very close to 90 • and thus unlikely to be stable in the presence of tides.Lastly, planets with nodal frequencies greater than the α syn range will not have significant Cassini state 2 obliquities (since ϵ ∼ I in that case).We also note that the frequency spectra are more complicated for some systems than others.For TOI-178, the frequency shifts are imperfectly accounted for in the average spectra.For TOI-1136, the complicated spectra are driven by chaotic orbital behavior.
Table 3 displays the results quantitatively.For each planet, we show the nominal obliquities associated with each peak frequency in that planet's spectrum, as described in Step 6 above.The nominal obliquities use k 2 /C = 1 within the calculation of α syn , but we also show the range of obliquities resulting from the interval [0.3α syn , 3α syn ].If the obliquity is < 1 • for the whole range (which is true when α syn ≪ g i ), we do not include it in the table.In this sense, the table displays only the set of significant Cassini state 2 equilibria.Again, we observe that many planets have several high obliquity equilibria.Some of the equilibria have a large range of possible obliquity values, which is due to the order-ofmagnitude range in α syn we are considering.As shown in Figure 1, the Cassini state 2 obliquity can vary greatly within a small range of |g|/α.The equilibrium obliquities are also summarized visually in Figure 12 (late in the paper because we require more concepts to be introduced).

Numerical spin-orbit simulations
Although the previous analysis yielded each planet's Cassini state 2 equilibria, not all of these equilibria are stable with an arbitrary degree of tidal dissipation.Here we test the stability of the equilibria using numerical integrations that account for the planets' tidal, spin, and orbital evolution.Our N -body code was developed and described in Millholland & Laughlin (2019).It uses a Bulirsch-Stoer integrator and is constructed in the framework of Mardling & Lin (2002), which is itself based on Eggleton et al. (1998).All bodies are endowed with structure (i.e. they are not point mass objects), and Jacobi coordinates are used to calculate the orbital evolution.In addition to the instantaneous Newtonian gravitational accelerations, the code also accounts for the accelerations on the planets from the star's quadrupolar gravitational potential, and the accelerations due to tides raised on the planets from the host star.The tidal accelerations are calculated in the framework of equilibrium tide theory in the viscous approximation (e.g.Eggleton et al. 1998;Leconte et al. 2010) and then used to calculate the torques on the spin vectors.Further details of the code are provided in Millholland & Laughlin (2019).
We perform the stability assessments by running a suite of simulations for each system in which all of the planets are initially placed in one of their possible equi-Table 3. Equilibrium obliquities resulting from the frequency analysis.The middle column in each sub-table shows the value of αsyn (α when ω = n) with k2/C = 1.The right column shows the Cassini state 2 obliquities, solved as described in Step 6 of Section 4.1, and the range of obliquities associated with the interval [0.3αsyn, 3αsyn].(Recall that we set I = 0.1 • , such that ϵ = 0.1 • is the minimum possible value of ϵ.)The right column also shows the gi frequencies that the equilibrium obliquities are associated with, with the labeling corresponding to the ordering presented in Table 2.
000001 -librium states.The stability is then tested by slowly decreasing the planet's tidal quality factor Q such that the tidal torque increases in strength.Eventually, for a low enough Q (strong enough tidal dissipation), the Cassini state breaks because the tidal torque overwhelms the resonant torque (Millholland et al. 2020;Su & Lai 2022a).By recording the Q upon breaking, we identify the limits of stability of each equilibrium state.We perform 100 simulations constructed using the following set-up: 1. Orbital and physical parameters: We randomly choose one of the 10 sets of orbital and physical parameters from the posterior distribution (masses, radii, periods, eccentricities and transit times or mean anomalies).These are the same parameter vectors as used in Section 4.1.We randomize the inclinations and longitudes of ascending node ).We fix k 2 = 0.3 and C = 0.3 as discussed earlier.
2. Spin parameters: For each planet, we randomly choose the initial obliquity from the set of possible equilibria identified in Section 4.1.We initialize the phase angle of the spin vector such that it is in resonance (with ŝ in the plane formed by n and k, and with ŝ and n on opposite sides of k).

Tidal parameters:
We set the planet's initial tidal quality factor to Q 0 = 10 6 .We then decrease Q according to with the final value set to Q f = 10 3 and the timescale τ Q = 10 4 yr.We emphasize that this is not intended to correspond to an actual physical evolution of the planet's tidal properties over time but is rather a means of identifying the minimum Q for which the equilibrium can still be stable.
The functional form of Q(t) is arbitrary, but the results are insensitive to this as long as the evolution is slow relative to the spin-axis precession (τ Q > α −1 ).

Simulation parameters:
We set the timestep equal to 1% of the innermost planet's period and run the integration for 0.2 Myr.
For each simulation, we monitor the planets' orbital and spin parameters over time.The initial spin states are not always stable.However, when they are, the obliquities oscillate around the equilibrium for some time before breaking out when Q is sufficiently small.We calculate the approximate time at which this occurs by finding when the obliquity drops to below 5% of the equilibrium.We then find the tidal quality factor, Q break , corresponding to the breaking point based on equation 7.
The Q break values aren't always the same from simulation to simulation, but they follow some general trends.Figure 5 shows the results of Q break as a function of the equilibrium obliquities that the planets were initialized with.We display the data for all of the systems at once, with the colors indicating different systems.It is clear that the simulations show a variety of outcomes in terms of how long the obliquities stay stable in their initial equilibrium positions, even for the same planet and the same equilibrium obliquity.These outcomes depend on factors such as the orbital inclinations, which are randomized for each simulation.However, for a given equilibrium, the minimum value (which corresponds to the longest stability timescale of the equilibrium) is bounded by a lower envelope.This envelope indicates the value of Q below which the tidal torque overwhelms the spinorbit resonant torque, and the equilibrium can no longer be maintained.
The tidal breaking limit can be calculated analytically using the formalism provided by Su & Lai (2022a).The timescale t s of tidal realignment of the spin vector is given by In other words, t s is the approximate timescale for the obliquity of a planet initialized outside of a Cassini state to damp down due to the effects of tides.Stronger tidal forces yield a shorter timescale.For a planet in Cassini state 2, t s cannot be arbitrarily small, since tides cause a phase shift of the equilibrium position of the spin vector out of the plane defined by the orbit normal vector and the total system angular momentum vector (Fabrycky et al. 2007), and Cassini state 2 ceases to exist once this phase shift is 90 • .(For more details, see Su & Lai (2022a).)The phase shift is < 90 • and Cassini state 2 is stable as long as t s > t s,c , where Setting t s = t s,c and using equation 2 (with the assumption I ≫ ϵ) and equation 4, one can show that the critical tidal quality for tidal breaking is given by This result shows that Q break,t depends only on the obliquity and the inclination and not on any other system parameters.Higher obliquities correspond to higher values of Q break,t (since higher obliquity states are harder to maintain), whereas higher inclinations correspond to lower values of Q break,t (since the resonant torque is stronger).The theoretical tidal breaking value, Q break,t , is shown with a solid line in Figure 5.It is a good match to the lower envelope when I = 0.04 • , which is close to the median mutual inclination between planets.The additional curve corresponding to I = 0.4 • shows that Q break,t can be significantly lower if the inclinations are higher.Broadly speaking, the results of the numerical spinorbit simulations indicate that very high obliquity states (≳ 45 • ) can be maintained for some planets, but generally they require relatively high tidal quality factors Q ≳ 10 4 if the inclinations are low.Tidal quality factors of Q ∼ 10 4 are commensurate with expectations for sub-Neptune-sized planets (Morley et al. 2017;Puranam & Batygin 2018;Millholland et al. 2020), although rocky super-Earths likely have smaller tidal quality factors.For instance, Brasser et al. (2019Brasser et al. ( , 2022) ) used both long-term dynamical simulations and interior modeling to constrain the tidal parameters of the TRAPPIST-1 planets (except for planets g and h), finding Q/k 2 in the range of ∼ 100 to a few 1, 000.Despite similar tidal heat fluxes to Solar System bodies like Io and Enceladus, the tidal quality factors of the TRAPPIST-1 planets are larger, which could be explained if they have a partially molten mantle and thus a less efficient conversion of tidal energy to interior heating.Nevertheless, for planets with Q ≲ 10 4 , high obliquity states are likely only maintained for orbital configurations with inclinations higher than those considered here.
The simulations in this section only considered obliquity evolution when planets were placed directly into their equilibrium states, as a means of performing the stability assessment.In reality, planets are not dropped directly into their equilibria; they need to encounter them through dissipation or resonance crossing.A planet starting at an arbitrary obliquity will evolve due to the combined effects of tidal damping and secular orbital perturbations.These considerations motivate our next question: can planets end up in high obliquity states when undergoing typical orbital evolution? 5. MIGRATION SIMULATIONS So far we have discovered that high obliquity states are possible for many planets in resonant chains, and they can be stable as long as the planetary tidal quality factors are large enough.However, in terms of understanding whether these high obliquity states are likely, it is more meaningful to investigate how they could have Figure 5. Breaking Q vs. equilibrium obliquity.After initially placing the planets at one of their equilibrium obliquities and slowly decreasing Q, the quantity Q break is the tidal quality factor at which the obliquity falls below 5% of the equilibrium.The colors represent planets from different systems.The solid black curve corresponds to the theoretical tidal breaking Q (i.e.Q break,t ) using I = 0.04 • , which is a good match to the lower envelope.The solid gray curve corresponds to the theoretical Q break,t using I = 0.4 • , which is shifted to lower values.The horizontal dashed line indicates the minimum Q considered in the simulation.originated during the system's initial assembly.Here, we explore case studies of two systems in our sample: Kepler-223 and TOI-1136.These systems are fairly different from one another; the 6-planet TOI-1136 system has a more complicated orbital frequency spectrum than the 4-planet Kepler-223 system (Figure 4).We model the formation/migration of the resonant orbital architectures and the simultaneous spin evolution.The degree of migration responsible in forming resonant chains is still debated, and there are a multitude of possible pathways that can describe any given system.We consider our simulations to represent plausible but nonunique histories of the Kepler-223 and TOI-1136 systems.We explore scenarios of both long-range migration (e.g. Lee & Peale 2002;Mills et al. 2016;Izidoro et al. 2017) and short-range migration (e.g.MacDonald & Dawson 2018).For each system, we adopt a migration set-up used in previous literature.This is so that we can investigate the obliquity evolution with zero fine tuning on our part.The specifics of the migration schemes may affect the details of our results, such as the exact fraction of cases that end up in high obliquity equilibria, but they will not affect our main conclusions.

Kepler-223
The formation of the Kepler-223 system was previously explored by Mills et al. (2016), and we where λi is the mean longitude (Siegel & Fabrycky 2021).The three-body resonant angles of this particular simulation do not agree perfectly with the observed data (Mills et al. 2016), but based on our simulation outcomes, there does not appear to be a significant link between the final critical angles and the migration history.Note that planets b and c form a transient 7:5 resonance before settling into the 4:3 resonance.The three vertical lines before 0.5 Myr indicate when a planet pair forms an expected resonance and the rightmost vertical line marks the approximate end of the migration.
emulate their migration set-up here.
We adapt the modify orbits direct migration model in RE-BOUNDx (Tamayo et al. 2020) and also apply some features from the type I migration module (Kajtazi et al. 2023), including adaptive migration timescales and an inner disk edge that halts the migration.We explore a variety of migration distances including short-range (∼ 0.1 − 0.2 AU), medium-range (∼ 0.2 − 0.4 AU), and long-range (∼ 0.4 − 0.7 AU).Using the WHFast integrator (Wisdom & Holman 1991;Rein & Tamayo 2015) within REBOUND (Rein & Liu 2012), we integrate the short-range and medium-range migration simulations for 3 Myrs and long-range scenarios for up to 5 Myrs.The details of the set-up and variations are provided in Appendix A, and a representative example of the orbital evolution during short-range migration is shown in Figure 6.We perform 304 simulations total.
We use the spin dynamics module recently developed for REBOUNDx by Lu et al. (2023) to track the planets' spin evolution during and after the migration.This module is based on the equilibrium tides framework of Eggleton et al. (1998).Our simulations assume a constant tidal time lag, which for the i-th planet was calculated according to τ i = (2Qn i ) −1 (Lu et al. 2023), where n i is the mean motion and Q = 10 4 the initial tidal quality factor.This is a reasonable approximation based on the Q of Neptune (Zhang & Hamilton 2008) and the few extrasolar sub-Neptunes with Q constraints (Morley et al. 2017;Puranam & Batygin 2018).We set C = 0.3 and k 2 = 0.3 for all planets (as in Section 4.2).The planets are initialized with randomized obliquities (consistent with an early epoch of collisions of planetary embryos; Miguel & Brunini 2010), and we set the initial spin rate of each planet to be 20% − 67% of the breakup spin rate, ω break = GM p /R 3 p .For consistency, we also set the following parameters for the star: C ⋆ = 0.1, k 2⋆ = 0.1, Q ⋆ = 10 6 , P ⋆ = 27 days.The stellar obliquity is set to 0 • and the tidal time lag is calculated according to τ ⋆ = (2Q ⋆ n 1 ) −1 .
Looking broadly across our set of 304 simulations, we find that all planets in Kepler-223 are capable of exciting and maintaining high obliquities.This confirms our hypothesis that the migration process naturally excites secular spin-orbit resonances that capture the planets in high obliquity states.We first discuss the obliquity dynamics of a representative example simulation (Figure 7, which corresponds to the orbital evolution shown in Figure 6) before discussing the simulation results as a whole.As seen in Figure 7, the nodal precession periods, T gi = 2π/|g i |, experience rapid "kicks" when a planet joins the resonant chain and stabilize once the chain has fully formed.They also increase around the time when the torque caused by the inner disk edge fully and swiftly halts the migration.The spin-axis precession periods, T α = 2π/(α cos ϵ) where α is calculated directly using equation 1, experience slow changes due to the combined effects of migration and tidal damping of the spin rate.
In the case of planets b and c, T α never crosses with one of the T gi curves, so the obliquities tidally damp down to near zero.Before this though, the obliquities experience a transient increase while T α increases; we suspect that this occurs as the planet attempts to maintain the spin-axis precession frequency, but further work would be necessary to understand this behavior fully.As for planets d and e, the obliquities first experience minor "kicks" (≲ 20 • ) caused by secular spin-orbit resonance encounters when T α crosses one of the T gi curves from below.They later get trapped in Cassini state 2 with a high obliquity when T α crosses one of the T gi curves from above and the secular spin-orbit resonance is captured.
While Figure 7 is a single example, the obliquity dynamics are fairly similar throughout our simulation suite, with the main differences being the final obliquity outcomes.These results are summarized in Figure  8, from which we observe several interesting features.First, the final obliquities are sensitive to the initial positions of the planets.This is consistent with the fact that the secular spin-orbit resonances are crossed during the migration process and are thus sensitive to the timing and scale of the migration.The obliquity behavior is most diverse when the semi-major axis of planet b is initialized in the range ∼ 0.14 − 0.32 AU.
Another feature seen in Figure 8 is that the final obliquities of all planets display several equilibria.Planets d and e show, respectively, two and three equilibria where the final obliquity is larger than 10 • .For planets b and c, the equilibria are less obvious, but they each display at least one high obliquity equilibrium.In all cases the final obliquities oscillate around an equilibrium with amplitudes typically smaller than a few degrees.These results are summarized in Table 4.The oscillation amplitudes presented in Table 4 and Figure 8 are calculated as in Millholland et al. (2018), Additional details of the simulations are provided in Appendix A. To summarize, we find that the short-range migration simulations produce the lowest fraction of successful simulations (i.e. in which the simulated resonant chain is the same as the observed system) but the widest diversity of obliquity outcomes.The long-range migration simulations are most readily capable of producing the observed resonant chain.Overall, the biggest takeaway is that secular spin-orbit resonances are robustly encountered and captured during the formation process, even considering a variety of migration distances and timescales.

TOI-1136
Whereas the migration in our Kepler-223 simulations took place over a few Myrs, it is important to explore the effects of more rapid migration, which we do here for the 6-planet resonant chain TOI-1136.We model our TOI-1136 simulations on Dai et al. (2023).Again using the modify orbits direct routine within REBOUNDx (Tamayo et al. 2020) and the WHFast integrator (Wisdom & Holman 1991;Rein & Tamayo 2015), we simulate the disk migration of the planets into their present-day positions, using an artificial inner disk edge to halt the migration.We explore both short-range (∼ 0.05 AU) and long-range (∼ 0.5 AU) migration scenarios.The details of the migration set-up are provided in Appendix A, and a representative example simulation of the orbital evolution is shown in Figure 9.
After finding a set of migration parameters resulting in analogous systems to the present-day TOI-1136, we run an additional set of simulations where, simultaneous with the disk migration, the planetary spin axes are evolved using the Lu et al. (2023) spin dynamics module in REBOUNDx.All planets are initialized with randomized obliquities between 0 • − 20 • , and the initial spin frequencies are set to between 10% and 25% of the breakup spin rates.The initial tidal quality factor is set to Q = 10 3 for all planets, which gives us an opportunity to observe interesting behaviors during the short lifespans of our simulations.We perform seven distinct simulations and run them for 8τ a , around ∼ 2 Myr.
Across the various spin simulations, we observe similar obliquity dynamics.All the planets are temporarily excited to very high obliquities but none are able to achieve a long-term stable high-obliquity state, with the possible exception of planet g, whose final behavior is unclear.An example simulation is presented in Figure 10, where we plot the obliquities along with the spin-axis precession periods (T α ) and nodal precession periods (T gi ).Because of the complex nodal frequency spectra for the TOI-1136 system (see Figure 4), we consider a larger set of g i frequencies than just the N p − 1 peak frequencies discussed previously.We first divide our data into 40,000-year segments, and for each partition, we construct nodal frequency spectra as described in Step 3 of Section 4.1.We do this for each planet and then combine peak frequencies to create a set of g i frequencies, resulting in an evolution in the nodal frequencies over time.Interestingly, four frequencies remain by the end of our simulations, which are very similar in all simulation cases.
The general behavior seen in Figure 10 can be summarized as follows.Following the migration, T α increases as a result of tidal damping of the planetary rotation rate.T α eventually starts to track one of the T gi curves, which corresponds to the planetary spin axis settling into Cassini state 1 or a low-obliquity Cassini state 2. Planet g does not reach a stable state, but we suspect that it would with a longer simulation time, based on the behaviors of planets e and f.Before the T α curves reach their final states, the obliquities undergo "kicks" caused by secular spin-orbit resonance crossings when T α crosses one of the T gi curves from below (e.g.Hamilton & Ward 2004).The resonant kicks are particularly dramatic for planets e and f, whose obliquities increase significantly before the planets settle into the low-obliquity equilibria.Planets b, c, and d experience comparatively milder increases, in part because of the speed with which their T α curves reach their final states.We also note a few "kicks" in planet g's obliquity that do not appear to correspond with frequency crossings.These may be explained by the presence of mixed mode resonant encounters, consisting of the interaction of T α with multiple orbital frequencies (Su & Lai 2022b).
Figure 10 is a single example of our spin simulations, but in all cases we observe similar obliquity dynamics.Planets b and c end with low but non-zero final obliquities, below 20 • in all cases, averaging 8.2 • and 7.0 • respectively.Planet d has the lowest final obliquity, averaging less than 0.1 • , and both planets e and f do not average higher than 5 • in their final obliquities.We do not observe the completion of planet g's obliquity evolution in any simulation.

Phase space diagrams
The migration simulations of Kepler-223 and TOI-1136 indicated fairly disparate spin dynamical behavior; the Kepler-223 planets frequently experienced capture into secular spin-orbit resonances and maintained high obliquities, while the TOI-1136 planets did not.In order to better understand these outcomes, it is helpful to explore phase space portraits of the planetary spin states.These portraits illustrate the trajectories in a simplified analytical framework.When the nodal precession rate (g = Ω) is uniform and e ≈ 0, the Hamiltonian that governs the planet's spin evolution is given by (Colombo 1966;Peale 1969;Ward 1975;Su & Lai 2022a) Here ϕ and cos ϵ are canonically conjugate variables, where ϕ is the phase angle of the spin vector around the orbit normal vector.Trajectories of evolution in the phase space (ϕ, cos ϵ) fall on level curves of H.
The top section of Figure 11 shows the phase space for the planets in the Kepler-223 system.We include the separatrices associated with the g i frequencies for each Cassini state 2 obliquity (from Table 3).The level curves correspond to I = 0.1 • in equation 12, but the resonance widths increase as a function of the inclination.Here the separatrices are close together, but they do not overlap.This indicates that spin states with trajectories close to the resonant fixed points can remain stable, whereas trajectories near the bounds of the separatrices may be unstable due to the effects of overlapping resonances that drive chaos.These phase space portraits thus explain why the Kepler-223 planets (particularly planet e) are able to maintain high obliquity states in our migration simulations (recall Figure 8).
Figure 11 also shows the Hamiltonian phase space for the planets in the TOI-1136 system.Unlike for Kepler-223, the resonant regions associated with the different equilibria are overlapping significantly.Such overlap drives chaos and prevents planets from maintaining stable high obliquity states.This agrees with our overall findings from the migration simulations.

Chaotic spin dynamics
Our migration simulations revealed that the TOI-1136 planets failed to maintain high obliquities despite the existence of several high obliquity equilibria.We just argued that this is likely a result of chaos; Figure 11 showed that several TOI-1136 planets have resonance regions that are overlapping for various equilibria.Resonance overlap is well-known to generate chaos (Chirikov 1979), and chaotic spin dynamics can drive large-amplitude obliquity variations, such as on Mars (e.g.Laskar & Robutel 1993;Touma & Wisdom 1993;Correia et al. 2003).Is chaos an important consideration for the other resonant chain systems as well?
In order to assess the role of chaos for the rest of our sample, we create phase space portraits for all systems as in Figure 11.We then summarize the resonance overlap results in Table 6 and in Figure 12.Out of the 25 planets with high obliquity equilibria, 12 have no overlap of the resonance regions of their various equilibria, and 13 have overlap for either some or all of the equilibria.The phase space portraits here use I = 0.1 • , but there would be more or less overlap for higher or lower inclinations, respectively.Overall, these results indicate that secular spin-orbit resonance overlap is fairly common for planets  3, where the dot at the center of each separatrix indicates the Cassini State 2 corresponding to the associated gi frequency.
For the TOI-1136 system, planets e, f, and g do not have high obliquity states.
in resonant chains.We expect this to drive chaos and limit the prevalence of stable high obliquity states.The planets may also experience large amplitude obliquity variations as a result of chaos, but this will be limited by tidal damping.

Impacts on planetary atmospheres
Planets that manage to maintain stable high obliquity states will experience consequences for their climates and atmospheres.First, high obliquities could enhance habitability by promoting a more balanced distribution of stellar radiation on the planet (e.g.Spiegel et al. 2009 3. We omit those with ϵ < 1 • .The dot shows the center of the resonance, and the lines show the resonance width (which we calculate as the width of the widest part of the separatrix).We note that this plot does not depict the range in the obliquities associated with the uncertainty in α shown in Table 3. Table 6.Resonance overlap results.Based on phase space portraits like those in Figure 11, the first column shows the list of planets that have no overlap between the resonant regions of their various equilibria.The second and third columns show planets with overlap between the resonant regions for all equilibria and for some of the equilibria, respectively.Finally, the last column shows planets with no high obliquity equilibria.Resonance overlap for some or all equilibria is fairly common.Colose et al. 2019).Even more crucially, non-zero obliquities make it such that the equilibrium rotation rate is sub-synchronous rather than synchronous (equation 4), thus preventing tidal locking.The large-scale atmospheric structure is also affected by the obliquity.Non-zero axial tilts cause the tidal bulge raised on the planet by the host star to shift across the planet each orbit.This leads to tidal heating, which can inflate a planet's atmosphere dramatically.Millholland (2019) showed that the radii of sub-Neptunes with H/He envelopes comprising only a few percent by mass can be inflated to twice the size they would otherwise have in the absence of obliquity-driven tidal heating.Radius inflation may be occurring for some planets in resonant chains.For instance, our migration simulations showed that Kepler-223 d and e frequently maintain enhanced obliquities.These two planets are pretty puffy; with radii of 5.25 R ⊕ and 4.63 R ⊕ , they are the largest planets in the system.It is beyond the scope of this paper to explore the potential connection between the physical properties of planets in resonant chains and their spin states, but this would be valuable future work.

Results on TRAPPIST-1
Out of all the systems in our sample, TRAPPIST-1 is certainly the one with the highest interest to the community given that it hosts multiple temperate, roughly Earth-size planets.The rotation rates and tidal dynamics of the TRAPPIST-1 planets have been studied by multiple other authors (e.g.Makarov et al. 2018;Dobos et al. 2019;Vinson et al. 2019;Brasser et al. 2019Brasser et al. , 2022;;Shakespeare & Steffen 2023), but their axial tilts have received little attention.We thus summarize our results on TRAPPIST-1 here.TRAPPIST-1 b, c, and d all have multiple high-obliquity equilibria (see Table 3), indicating that it is theoretically possible for them to be trapped in one of these states.However, we deem it to be unlikely.Agol et al. (2021) showed that TRAPPIST-1 is extraordinarily coplanar, with a mutual inclination dispersion of around ∼ 0.04 • .Such low inclinations would require high Q values for the planets, Q ≳ 3 × 10 3 , in order for Cassini states to be stable with obliquities in excess of ∼ 45 • .However, the Q values of the TRAPPIST-1 planets are likely smaller than this (e.g.Brasser et al. 2019Brasser et al. , 2022)).
TRAPPIST-1 d is the only planet that may stand a chance of maintaining a non-negligible obliquity, since it possesses an equilibrium at ∼ 12 • , which could be stable with Q ∼ 10 3 .The resonant region associated with this equilibrium has no resonance overlap in phase space (Table 6), so it could avoid chaos.However, the planet would have had to find this equilibrium through an early resonant encounter.We have not studied the probability of this occurring in detail, but it is likely a low probability event.These results are important to consider in future studies of the atmospheres and habitability of the TRAPPIST-1 planets.

CONCLUSION
Resonant chains are extraordinary systems that warrant detailed investigations into their properties.The atmospheres, habitability, and even orbital and interior properties can all be affected by the dynamics of the planetary spin vectors.Given that resonant chains most likely formed via convergent orbital migration, their orbital precession frequencies would have been modified in the past.The orbital frequencies may have crossed through resonant commensurabilities with the natural frequencies of precession of the planetary spin axes, at which point the spin axes could be tilted to high obliquities via secular spin-orbit resonances (Section 2).We set out to test this hypothesis in this work by performing a comprehensive exploration of the spin dynamics of planets in observed resonant chains and identifying the prevalence of high obliquity equilibria.
We considered a sample of eight observed resonant chains, each with three or more planets in a sequence of mean-motion resonances (Section 3, Figure 2).We first analytically identified the spin equilibria with high obliquities.These are configurations in which the planet could maintain a stable high obliquity, even in the presence of tidal dissipation.We find such states to be very common: 25 out of the 42 planets in the sample have at least one high obliquity (≳ 20 • ) equilibrium, and seven out of the eight systems have at least one planet with a high obliquity equilibrium (Table 3, Section 4.1).A visual summary of these equilibrium obliquities is shown in Figure 12.
Although these states are common, not all of them can be stable for arbitrary tidal parameters.When tides are strong, they can overwhelm the resonant torque that otherwise maintains the high obliquity.We performed a stability assessment by placing the planets in their equilibrium states and evolving the systems using numerical spin-orbit simulations (Section 4.2).We found that the minimum stable Q (Figure 5) is well-described by an analytic result that just depends on the obliquity ϵ and inclination I (equation 10, adapted from Su & Lai 2022a): For a high obliquity ϵ ∼ 45 • and low inclination I = 0.05 • , this works out to about Q ≳ 3000.Most notably, this places stringent constraints on the TRAPPIST-1 system; planets b, c, and d all have theoretical high obliquity equilibria, but planet d is the only one that possesses an equilibrium (at ∼ 12 • ) that might be stable.
Our initial analysis only identified the existence of the spin equilibria but did not explore the likelihood of the planets attaining them.We considered two systems (Kepler-223 and TOI-1136) as case studies and simulated their formation via convergent migration and the resulting planetary spin evolution (Section 5) using a few hundred simulations.The Kepler-223 planets frequently encountered secular spin-orbit resonances during their migration; we found all four planets to be capable of maintaining stable high obliquity states, with planet e being most susceptible (Figures 7 and 8).On the other hand, the TOI-1136 planets did not readily maintain high obliquities despite encountering multiple secular spin-orbit resonances.We showed that these results are intuitive in the context of phase space portraits of the spin evolution, since the TOI-1136 planets show significant resonance overlap, whereas the Kepler-223 planets do not (Figure 11).Considering all systems in our sample, we find that resonance overlap is fairly common, showing up in about half of the planets with high obliquity equilibria.This resonance overlap drives chaos and can hinder the stability of high obliquities.We thus expect chaos to affect the spin evolution of some planets in resonant chains, although we suggest further study of the competing influences of chaos and tidal evolution.
In cases where planets do maintain high obliquities over long timescales, their atmospheres can be dramatically impacted, including through inflation driven by heating from obliquity tides.Tidal heating can puff up sub-Neptunes to twice the size they would otherwise have.This can have observable effects for systems like Kepler-223 that may readily maintain high obliquities.In this way, spin dynamics provide a crucial link between the planetary orbital evolution and physical properties, and they allow insight into the climate and atmospheric circulation.This may become more influential in the near future as JWST continues to advance the physical and chemical characterization of exoplanets.form a 8:5 resonance instead with b and c forming a 5:4 resonance.The frequent occurrence of higher order resonances might be the result of a relatively slow migration speed compared to medium-range or long-rage migration scenarios.

A.1.2. Medium-range migration
In the medium-range migration simulations, the initial a b is sampled according to a b ∼ Unif[0.2,0.4] AU.In general, the planets achieve the resonant chain within < 0.1 Myr and stop migrating ∼ 2.3 − 2.6 Myr after the start of simulation.Out of 82 medium-range migration simulations, the system achieves the currently observed configuration in 51 scenarios and becomes unstable in only one.In the remaining 30, b and c are again prone to forming higher-order resonances (7:5, 10:7, or 11:8) while d and e form a higher order resonance in only one case.

A.1.3. Long-range migration
We perform 85 long-range migration simulations where the initial a b is set according to a b ∼ Unif[0.4,0.7] AU.The planets form the resonant chains within less than ∼ 0.05 Myrs and stop migrating within ∼ 2.6 − 2.8 Myrs.In 67 simulations, the system achieves the observed configuration.There are no unstable cases.In the remaining 18 simulations, planets b and c form a 7:5 resonance.Planets d and e form a 7:5 resonance in only one scenario.The long-range migration scenarios thus display the most stable outcomes out of all variations, likely because of the smaller occurrence of high-order resonances.

A.2. Distinct sets of nodal frequencies
Our 304 simulations -including short-range, medium-range, and long-range migration scenarios -can be divided into three distinct sets based on their final nodal frequencies and eccentricities.Each of the three sets of nodal frequencies is associated with exactly one set of eccentricities, and vice versa.Each set is also associated with two well-defined libration centers of the three-body resonant angles.In addition, the obliquity behavior varies noticeably across the three sets.These results are summarized in Table 7 and Figure 13.
The nodal frequencies differ somewhat from those in Section 4.1.This could be the result of our simulations achieving slightly different eccentricities and inclinations from the ones used in Section 4.1, since the details of the migration affect the evolution of semi-major axis, eccentricity, and inclination, all of which contribute to the values of nodal frequencies.In addition, the nodal frequencies could be affected by our implementation of the inner disk edge.In fact, we can see from Figure 7 that nodal frequencies adjust when the migration is halted.
The three-body angles shown in Table 7 differ slightly from theoretical values found in previous studies (Siegel & Fabrycky 2021;Delisle 2017).The three-body angles ϕ bcd = 155.4• and ϕ cde = 119.7 • associated with set 1 come the closest to the observed values and, interestingly, all four eccentricities present in set 1 are in agreement with the observed data (Mills et al. 2016).This is also the three-body resonant angle configuration the system is most likely to achieve (in 48/172 ≈ 28% of successful simulations) which is in agreement with Delisle (2017).However, simulations with this final three-body angle configuration were less likely to maintain high obliquities on average.Out of 48 simulations, 23 ended up with at least one planet maintaining a high obliquity (> 10 • ); none of the simulations had more than two planets maintain high obliquities.
Table 7.The three distinct sets of nodal frequencies in the Kepler-223 system, along with their associated final eccentricities, libration centers of the three-body critical angles, and the frequency of such simulation outcomes.The orbital parameters are calculated as means of their values during last 10,000 years of the simulations.

A.3. TOI-1136
For TOI-1136, we retain the orbital and physical parameters from our spin-orbit simulations in Section 4.2, except the initial semi-major axes, which we start wide of the present-day positions.We initialize planet b between 5% and 25% outside its current orbit and set the other planets' semi-major axes such that the period ratios are between 2% and 25% wide of the current resonances.This allows us to test both short-range (moving planets ∼ 0.05 AU) and long-range (moving planets ∼ 0.05 AU) migrations.Eccentricities are initialized at 0 and inclinations are randomized within 0 • − 5 • .Like Dai et al. (2023), we use an inner disk edge at 0.05 AU with a width of 0.01 AU to halt the migration.We vary the disk surface density at 1 AU (Σ 1au ) within 5 − 10 3 g cm −2 randomly in logarithmic space, assuming Σ = Σ 1au a −1.5 .We vary the scale height, h ≡ H/R, between 0.01 and 0.1, kept constant throughout the disk.Converting these values to semi-major axis decay and eccentricity damping timescales (Pichierri et al. 2018), ; τ e ≃ 1 0.780 ) we obtain migration timescales between 0.005 Myr and 2.25 Myr.We generally evolve the system over 2τ a , although for longer migration timescales we resort to shorter simulations due to computational constraints.However, we find by visual inspection of period ratios and convergence of critical angles that migration finishes well before 0.5τ a in almost all cases.

A.3.1. Short-range migration
We find short-range migration to be the favored scenario to result in resonant chains resembling the observed TOI-1136 system.We define analogous simulations as those in which all planets ended in MMR with less than a 0.2% deviation from their present-day period ratios.Roughly 200 of our 1500 short-range migrations created analogous systems.Of these, shorter migration timescales are most successful, with inner disk densities of around 10 g cm −2 and scale heights near 0.03 giving a migration timescale of τ a ∼ 0.015 − 0.03 Myr.This relatively short timescale is explained by the inner disk edge, halting the fast migration of innermost planets and allowing outermost ones to fall into resonance (Izidoro et al. 2017).An example of one of our successful short-range migration simulations is shown in Figure 9.The middle panel shows that some planetary pairs initially undergo divergent migration, but once the innermost planets begin to reach the inner disk edge, this divergent migration becomes convergent.This is clearest for planets pairs d and e and f and g, which do not begin convergent migration until the innermost planet of the pair achieves its own MMR state.

A.3.2. Long-range migration
The long-range migration simulations are significantly less successful.Of around 250 long-range migration simulations (the exact number is unclear, as a small portion of simulations evolved into chaotic behavior), none resulted in resonant chain formation.This mirrors previous migration simulation results of Dai et al. (2023), which similarly failed to achieve MMR from longer-range migration.We also find that the weak 7:5 MMR of planets e and f is easily skipped or perturbed during longer-range migration.

Figure 1 .
Figure 1.Equilibrium positions of the planetary obliquity in the four Cassini states.The obliquities are plotted as a function of |g|/α.The inclination in equation 2 is set to I = 1 • .For |g|/α > (|g|/α)crit (equation 3), states 1 and 4 disappear, and the obliquity of state 2 tends towards ϵ = I (horizontal dotted line).States 1 and 2 are the only stable states in the presence of tidal dissipation, so stable high obliquities correspond to state 2.

Figure 2 .
Figure2.Architectures of the resonant chains in our sample.The planets are spaced according to their semimajor axes.The dot sizes are proportional to the planet sizes, and they are colored according to the planet masses.In between each near-resonance pair, we label the approximate period ratio.

Figure 3 .
Figure 3. Nodal frequency spectra for the TRAPPIST-1 planets.For each planet, we show the averaged spectrum of the x-component of n, as described in Step 4 in Section 4.1.The planets are ordered from smallest to largest orbital period.The dotted vertical line indicates the value of αsyn, and the shaded region indicates the interval [0.3αsyn, 3αsyn].The dots in each panel indicate the peak frequencies identified for that planet, as described in Step 5 in Section 4.1.The small thick tick marks on the bottom axis indicate the Np − 1 = 6 locations of the gi frequencies.

Figure 4 .
Figure 4. Nodal frequency spectra for all systems.For each planet, we show the averaged spectrum of the x-component of n, as described in Step 4 in Section 4.1.The planets within a given system are ordered from smallest to largest orbital period.The dotted vertical lines indicates the value of αsyn, and the shaded regions indicate the interval [0.3αsyn, 3αsyn].The dots in each panel indicate the peak frequencies identified for that planet, as described in Step 5 in Section 4.1.The small ticks on the bottom axis indicate the Np − 1 locations of the gi frequencies.(TRAPPIST-1 was already shown in Figure 3 but is included again for convenience.)

Figure 6 .
Figure 6.Example of a simulated migration history of the Kepler-223 system.From top to bottom, the panels display the orbital periods, period ratios, and three-body resonant angles.The resonant angles are calculated according to ϕ bcd = 3λ b − 6λc + 3λ d and ϕ cde = 2λc − 6λ d + 4λewhere λi is the mean longitude(Siegel & Fabrycky 2021).The three-body resonant angles of this particular simulation do not agree perfectly with the observed data(Mills et al. 2016), but based on our simulation outcomes, there does not appear to be a significant link between the final critical angles and the migration history.Note that planets b and c form a transient 7:5 resonance before settling into the 4:3 resonance.The three vertical lines before 0.5 Myr indicate when a planet pair forms an expected resonance and the rightmost vertical line marks the approximate end of the migration.

Figure 7 .
Figure 7. Dynamical evolution of the precession periods and obliquities in Kepler-223 corresponding to the migration history shown in Figure 6.Each panel corresponds to a different planet.Blue solid lines represent the spin-axis precession period Tα = 2π/(α cos ϵ), where α is calculated according to equation 1.The black dotted lines represent the three nodal precession periods Tg i = 2π/|gi| of the system, calculated using the frequency method described in Section 4.1.The colored solid lines in each panel represent the planetary obliquities.The first three vertical lines indicate the moments when a planet pair forms the expected MMR, and the vertical black lines on the bottom two panels indicate significant crossings of Tα and Tg i .

Figure 8 .
Figure 8. Final obliquities of the Kepler-223 planets across all simulations.The obliquities are plotted as a function of the initial semi-major axis of planet b.The dots represent the final equilibrium obliquities of each planet, and the errorbars represent the obliquities' oscillation amplitudes.We note that the oscillation amplitudes depend on our chosen inclinations.

Figure 9 .
Figure9.Example of a simulated migration history of the TOI-1136 system.From top to bottom, the panels show the time evolution of the orbital periods, the period ratios, and the two-body critical angles, calculated as ϕij = qλi − pλj + (p − q)ϖij.The inner disk edge halts the inward migration of the innermost planets and allows convergent migration of the outermost planets.Shown here is only the relevant portion of the simulation, not the entire timescale.

Figure 10 .
Figure 10.Dynamical evolution of the precession periods and obliquities in TOI-1136 corresponding to the migration history shown in Figure 9.Each panel represents a different planet, going from planet b at the top to planet g at the bottom.The dashed black lines indicate the precession periods corresponding to the nodal frequencies, Tg i = 2π/|gi|.The solid black lines represent the spinaxis precession periods, Tα = 2π/α cos ϵ.The colorful lines indicate the planetary obliquities, corresponding to the right y-axes.The vertical shaded regions have been added manually to mark obliquity spikes caused by resonance crossings of the orbital and spin-axis precession periods.Note minor obliquity spikes in planet g not related to these crossings.

Figure 11 .
Figure 11.Phase space portraits for the planets in the Kepler-223 and TOI-1136 systems.The thick dashed lines indicate the separatrices of H for I = 0.1 • , with the interior of the separatrix shaded in.We show the separatrices associated with each gi frequency indicated in Table3, where the dot at the center of each separatrix indicates the Cassini State 2 corresponding to the associated gi frequency.For the TOI-1136 system, planets e, f, and g do not have high obliquity states.

Figure 12 .
Figure12.Summary of equilibrium obliquities and resonance widths.For each planet in our sample, we show the Cassini state 2 obliquities from Table3.We omit those with ϵ < 1 • .The dot shows the center of the resonance, and the lines show the resonance width (which we calculate as the width of the widest part of the separatrix).We note that this plot does not depict the range in the obliquities associated with the uncertainty in α shown in Table3.

Figure 13 .
Figure 13.Final obliquities and eccentricities of the Kepler-223 planets, split up into the three distinct sets of nodal frequencies.From top to bottom, each set of two panels represents the obliquities and eccentricities associated with one collection of nodal frequencies.The colors correspond to different planets, with the same scheme as Figure 8.The likelihood of maintaining high obliquities is noticeably different for each set.

Table 2 .
Nodal precession frequencies.We show the set of Np − 1 frequencies, {gi}, for each system.These are calculated by finding the unique frequencies among the peak frequencies identified in Step 5 of Section 4.1.default value as mentioned earlier, but since k 2 and C are uncertain, we display the results for the interval [0.3α syn , 3α syn ].Planets for which the peak nodal frequencies are close to their α syn range can have Cassini state 2 obliquities that are enhanced and not too close to 90 the • .The planets that are most interesting in this regard are TRAPPIST-1 c, d, e, and f, which all have peak frequencies falling in their α syn range.

Table 4 .
Final obliquities of the Kepler-223 planets across all simulations.Only equilibria for obliquities > 10 • are listed.

Table 5 .
The relative occurrence of observed final combinations of high obliquity states in Kepler-223 migration simulations.