Kepler-discovered Multiple-planet Systems near Period Ratios Suggestive of Mean-motion Resonances Are Young

Before the launch of the Kepler Space Telescope, models of low-mass planet formation predicted that convergent type I migration would often produce systems of low-mass planets in low-order mean-motion resonances. Instead, Kepler discovered that systems of small planets frequently have period ratios larger than those associated with mean-motion resonances and rarely have period ratios smaller than those associated with mean-motion resonances. Both short-timescale processes related to the formation or early evolution of planetary systems and long-timescale secular processes have been proposed as explanations for these observations. Using a thin disk stellar population’s Galactic velocity dispersion as a relative age proxy, we find that Kepler-discovered multiple-planet systems with at least one planet pair near a period ratio suggestive of a second-order mean-motion resonance have a colder Galactic velocity dispersion and are therefore younger than both single-transiting and multiple-planet systems that lack planet pairs consistent with mean-motion resonances. We argue that a nontidal secular process with a characteristic timescale no less than a few hundred Myr is responsible for moving systems of low-mass planets away from second-order mean-motion resonances. Among systems with at least one planet pair near a period ratio suggestive of a first-order mean-motion resonance, only the population of systems likely affected by tidal dissipation inside their innermost planets has a small Galactic velocity dispersion and is therefore young. We predict that period ratios suggestive of mean-motion resonances are more common in young systems with 10 Myr ≲ τ ≲ 100 Myr and become less common as planetary systems age.


INTRODUCTION
Gravitationally mediated interactions between newly formed planets and their parent protoplanetary disks have long been predicted to cause planets to radially migrate within their parent protoplanetary disks (e.g., Goldreich & Tremaine 1980;Ward 1986Ward , 1997;;Lin & Papaloizou 1986).Further studies of these processes led to the prediction that the trajectories of migrating Corresponding author: Jacob H. Hamer jacobhhamer@gmail.complanets would converge and result in dynamically important planet-planet interactions (e.g., Bryden et al. 2000;Kley 2000;Masset & Snellgrove 2001).The Dopplerbased discovery of the resonant multiple giant planet system GJ 876 (Marcy et al. 2001) confirmed these predictions for convergent migration and showed that they can lead to resonant configurations.More detailed studies of the resonant capture process have since suggested that capture into low-order mean-motion resonances is a likely outcome if migration rates are slow and eccentricity damping efficient (e.g., Snellgrove et al. 2001;Lee & Peale 2002;Nelson & Papaloizou 2002;Thommes 2005;Kley et al. 2005;Quillen 2006;Pierens & Nelson 2008;Lee et al. 2009;Ogihara & Kobayashi 2013).Planet-planet scattering rarely leads to meanmotion resonances (e.g., Raymond et al. 2008).Though there are no planet-planet mean-motion resonances in the solar system1 , multiple giant exoplanet systems with period ratios suggestive of mean-motion resonances are more common than expected if orbital periods were uncorrelated (e.g., Wright et al. 2011).
The successful explanation of mean-motion resonances in multiple giant planet systems by the convergent migration scenario motivated its extension to systems of low-mass planets and the expectation that mean-motion resonances would be relatively common in systems of low-mass planets (e.g., Papaloizou & Szuszkiewicz 2005;Cresswell & Nelson 2006;Terquem & Papaloizou 2007).This turned out not to be the case.Data from the Kepler Space Telescope (Borucki et al. 2010) revealed that while multiple small planet systems are more likely to be found with period ratios suggestive of first-order meanmotions resonances than if period ratios were randomly distributed, period ratios suggestive of mean-motion resonances are not common (e.g., Lissauer et al. 2011;Fabrycky et al. 2014).Instead, planet pairs tend to favor period ratios just wide of resonances and tend to avoid period ratios just interior to resonances.
The properties of the Kepler-discovered multiple small planet system period distribution have been attributed to short-timescale processes related to the formation or early evolution of planetary systems.Multiple small planet systems that formed in situ would have no obvious reason to be found with period ratios suggestive of mean-motion resonances, though planet-planet interactions might still result in mean-motion resonances during in situ formation (e.g., Petrovich et al. 2013).Meanmotion resonances initially present in systems formed via convergent migration might be forced out of resonance by turbulence in their parent protoplanetary disks (e.g., Adams et al. 2008;Rein 2012;Batygin & Adams 2017) or by density waves generated by the planets themselves (e.g., Podlewska-Gaca et al. 2012;Baruteau & Papaloizou 2013;Cui et al. 2021).The migration and eccentricity damping timescales present in real protoplanetary disks may make resonant configurations transient (e.g., Goldreich & Schlichting 2014;Charalambous et al. 2022;Laune et al. 2022).Dynamical interactions with residual planetesimals can produce period ratios wide of mean-motion resonances as well (e.g., Thommes et al. 2008;Chatterjee & Ford 2015;Ghosh & Chatterjee 2023).
Long-timescale secular processes have also been proposed to explain the properties the Kepler-discovered multiple small planet system period distribution, especially secular interactions between near-resonant planets in the presence of weak dissipation.Because tidal dissipation can remove orbital energy from inner planets much more efficiently than from outer planets, a planet pair subject to differential tidal dissipation will naturally evolve to a larger period ratio (e.g., Papaloizou 2011;Lithwick & Wu 2012;Batygin & Morbidelli 2013b;Delisle et al. 2012Delisle et al. , 2014;;Lee et al. 2013).As expected in this tidal dissipation scenario, Delisle & Laskar (2014) found that as the period of the innermost planet in a Kepler-discovered system of small planets increases, the excess of planet pairs exterior to period ratios suggestive of mean-motion resonances decreases.In contrast, Lee et al. (2013) showed that tidal dissipation must be an order of magnitude more efficient than expected to explain the properties of the Kepler-discovered multiple small planet period distribution.Silburt & Rein (2015) further argued that in the tidal dissipation scenario the initial eccentricities required to explain the occurrence of planets exterior to 2:1 period ratios are unreasonably high.
A comparison of the relative ages of Kepler-discovered small planet systems close to/far from period ratios suggestive of mean-motion resonances would help differentiate between these possibilities.Since most Keplerdiscovered systems of small planets orbit mature stars with ages τ ≳ 1 Gyr, the lack of an age offset between systems with planets close to/far from period ratios suggestive of mean-motion resonances would imply that systems move out of mean-motion resonances early in their histories at ages τ ≲ 1 Gyr.On the other hand, the observation that systems with planets far from period ratios suggestive of mean-motion resonances are older than those systems with planets close to period ratios suggestive of mean-motion resonances would support the idea that long-timescale secular processes move systems out of mean-motion resonance.Additional secular processes beyond tidal dissipation would be necessary if systems with both small-separation and large-separation innermost planets displayed the same relationship between relative age and separation from period ratios suggestive of mean-motion resonances.
In this article, we use the Galactic velocity dispersion of a thin disk stellar population as a proxy for its age and find that Kepler-discovered systems of small planets with period ratios suggestive of second-order mean-motion resonances are younger than both singletransiting and multiple-planet systems that lack planets with period ratios suggestive of mean-motion resonances.We suggest that non-tidal secular processes operating over hundreds of Myr to Gyr are responsible for driving planet pairs away from period ratios suggestive of second-order mean-motion resonances.Among systems with at least one planet pair near a period ratio suggestive of a first-order mean-motion resonance, only the population of systems likely affected by tidal dissipation inside their innermost planets have a small Galactic velocity dispersion and are therefore young.We describe in Section 2 the assembly of our analysis sample.We detail in Section 3 the methodology we use to identify systems with period ratios suggestive of meanmotion resonances, evaluate the importance of tides for a system's evolution, and compare the relative ages of systems close to/far from period ratios suggestive of meanmotion resonances.We review the explanations for and implications of our analyses in Section 4. We conclude by summarizing our findings in Section 5.

DATA
We use in our analysis a catalog of Kepler planet candidates optimized for completeness as well as the detection and characterization of planet candidates in the presence of transit-timing variations (Lissauer et al. 2023). 2 Our sample includes all confirmed and candidate planets with a disposition 'P' in the Lissauer et al. (2023) catalog.Because virtually all Keplerdiscovered planet candidates in multiple-planet systems are true planets (e.g., Lissauer et al. 2012), we analyze all multiple-planet systems regardless of whether their constituent planets have been classified as planets or planet candidates.
Following Schlaufman & Halpern (2021), we calculate updated planet radii using the observed transit depths in our input catalog and the stellar radii derived in Brewer & Fischer (2018), Johnson et al. (2017), or Berger et al. (2020) (in order of decreasing priority).We account for both transit depth and stellar radius uncertainties in our updated planet radius uncertainties.We collect Gaia DR3 source id information for our sample using astroquery (Ginsburg et al. 2019) to query SIM-BAD (Wenger et al. 2000) for each Kepler Input Catalog 2 We have verified that our conclusions would be unchanged if we instead used the Thompson et al. (2018) Kepler Data Release 25 (DR25) Kepler Objects of Interest (KOI) list optimized for uniformity and statistical analyses and rejected all planet candidates classified as false positives.
(KIC) identifier in our sample.We use each Gaia DR3 source id to collect all other relevant information for each star in our catalog from the Gaia Archive 3 and require parallax over error > 10 to ensure the precision of our velocity dispersion inferences.
In addition to Gaia astrometry, we need radial velocities to calculate the space velocities of individual stars and thereby the velocity dispersion of our sample.In order of decreasing priority, we use the radial velocities from the California-Kepler-Survey (CKS; Petigura et al. 2017), the Apache Point Observatory Galactic Evolution Experiment (APOGEE) 4 , DR7 of the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) Medium-Resolution Spectroscopic Survey (MRS; Cui et al. 2012;Zhao et al. 2012), Gaia DR3 (Katz et al. 2023), and DR7 of the LAMOST Low-Resolution Spectroscopic Survey (LRS; Cui et al. 2012;Zhao et al. 2012).We combine multiple APOGEE radial velocities in the table apogeeStar for a single object using a vscatterweighted average.We combine multiple LAMOST radial velocities for a single object using a weighted average, weighting by radial velocity uncertainty for the MRS and by the g-band spectral signal-to-noise ratio for the LRS.If the difference between the maximum and minimum LAMOST radial velocities for a single object is larger than three times the uncertainty of the weighted average, we do not use the LAMOST radial velocity.Following Marchetti et al. (2022), to ensure reliable Gaia DR3 radial velocities we require rv nb transits > 10 and rv expected sig to noise > 5.
by the Keplerian orbits of stellar-mass companions with unimodal posteriors in the APOGEE-based catalog presented in Price-Whelan et al. (2017, 2020) The presence in the system of an object with these Keplerian orbital parameters-possibly a non-transiting, highly eccentric brown dwarf-is difficult to reconcile the properties of the confirmed planet Kepler-1717 b.In this case we decided that the transitbased system parameters are likely more realistic and therefore keep the system in our sample.
While the vast majority of known planet host stars are members of the Milky Way's thin disk stellar population, the presence in our comparison or analysis samples of planet host stars belonging to the thick disk or halo stellar populations could significantly inflate their velocity dispersions.We identify the star Kepler-292/KOI-1364/KIC 6962977 as a kinematic member of the thick disk after using galpy6 (Bovy 2015) to integrate its orbit in the McMillan (2017) Milky Way potential over several Gyr.We find that Kepler-292/KOI-1364/KIC 6962977 has z max ≳ 3 kpc and e = 0.47 and therefore remove it from our sample.Our final analysis sample of Keplerdiscovered transiting multiple-planet systems with the necessary data to calculate precise Galactic kinematics includes 576 systems with 1460 planets.

Identification of Systems with a
Plausibly-resonant Planet Pair It is non-trivial to prove that a pair of planets with a period ratio suggestive of a mean-motion resonance is actually resonant with librating orbital elements.Detailed dynamical analyses have shown that planet pairs with period ratios close to but not exactly equal to rational numbers may indeed be resonant with librating orbital elements.This is the case for GJ 876 (e.g., Peale 1976;Rivera & Lissauer 2001).Similarly, planet pairs with period ratios apparently equal to rational numbers may in reality have circulating orbital elements (e.g., Veras & Ford 2012).In short, detailed dynamical modeling is required to prove that a pair of planets is resonant with librating orbital elements.These dynamical models critically depend on planet masses and orbital eccentricities that are often unknown or poorly constrained for lowmass planets that represent the majority of known exoplanet systems.As a result, approximate methods are often necessary to evaluate whether or not a pair of lowmass planets can be thought of as plausibly resonant.
Because Kepler-discovered systems of small planets rarely have the precise constraints on planet masses and orbital elements necessary for rigorous dynamical modeling, heuristic methods have often been used to identify pairs of planets plausibly in mean-motion resonances.Virtually all heuristic methods in the literature assume that planet pairs with period ratios closest to rational numbers are also most likely to be in mean-motion resonances.We use two different heuristic methods to identify plausibly resonant pairs of planets, one that is more physically motivated but that depends on mostly unknown planet masses and orbital eccentricities and one that is less physically motivated but independent of planet masses and orbital eccentricity.As we will show, we reach the same conclusions using either methodology.
The first method we use to identify plausibly resonant planet pairs is the physically-motivated but planet massand orbital eccentricity-dependent heuristic method proposed by Steffen & Hwang (2015).The Steffen & Hwang (2015) approach considers planet pairs with period ratios separated from first-order mean-motion resonances by a small number of libration widths plausibly resonant, while planet pairs with period ratios separated from first-order mean-motion resonances by a large number of libration widths non-resonant.For each planet pair we calculate libration widths normalized to semimajor axis using Equation (8.76) of Murray & Dermott (1999) The ratio C r /n of the resonant part of the disturbing function C r to the orbital frequency n is defined as (m ′ /m c )αf d (α), where m c is the mass of the host star, m ′ is the mass of the outer planet in a planet pair, and αf d (α) are numerical factors from Table 8.5 of Murray & Dermott (1999).The quantity j 2 = −p where p is the denominator of the period ratio (e.g., 1 for a 2:1 meanmotion resonance).We use the Lissauer et al. (2011) mass-radius relation M p = (R p /R ⊕ ) 2.06 to infer masses for each planet in a pair based on their radii and calculate the value of the inner planet's eccentricity that minimizes libration width.This latter assumption ensures that we are as conservative as possible with our inferred libration widths.We plot in the top panel of Figure 1 the distribution of separations from the closest first-order mean-motion resonance in units of libration widths.We note that the libration width relation was derived under the assumption of a massive outer planet and a massless inner planet.While this assumption is not strictly valid, Steffen & Hwang (2015) reasoned that since libration width generally scales with total planet mass (e.g., Deck et al. 2013) this approach is still valid.7 The Steffen & Hwang (2015) heuristic method for the identification of plausibly resonant planet pairs described above is planet mass-and orbital eccentricitydependent.An alternative planet mass-and orbital eccentricity-independent heuristic is the ϵ parameter defined for first-and second-order resonances as where P out /P in is the ratio of the outer planet period to the inner planet period.The ϵ parameter is often used to describe the dynamics of systems near mean-motion resonances, and predictions for the evolution of ϵ away from resonances in initially resonant systems as a result of tidal dissipation or planetesimal scattering have been published (Delisle & Laskar 2014;Chatterjee & Ford 2015).The ϵ parameter does not account for the varying widths of resonances though, so a fixed separation in ϵ could be far from a relatively narrow 6:5 resonance but close to a relatively wide 2:1 resonance.In contrast, the ζ parameter described in Lissauer et al. (2011) and Fabrycky et al. (2014) defined for first-and second-order resonances as accounts for the varying widths of first-and secondorder resonances and is able to highlight asymmetries in the period ratio distribution of Kepler-discovered systems of small planets near first-order resonances.There is no physically-motivated way to specify the probability that a planet pair is truly in resonance as a function of ζ though.
The second method we use to identify plausibly resonant planet pairs is an extension of the ϵ parameterbased method described above.We define the parameter where p res = (j + i)/j and i ∈ (1, 2).This δ res parameter differs from the ϵ parameter used in Delisle & Laskar (2014) and Chatterjee & Ford (2015) in that it is normalized by the resonant period ratio.While this normalization has been neglected in theoretical analyses of departures from resonances (e.g., Lee et al. 2013), we include this normalization to account for the fact that widely-spaced resonances have larger "widths" in period ratio than more closely-spaced resonances.For example, a 3% fractional period-ratio separation from a 2:1 resonance would extend over 0.06 in period ratio while the same fractional period-ratio separation from a 6:5 resonance would extend over only 0.036 in period ratio.We plot in the bottom two panels of Figure 1 the distributions of separation from the closest first-or second-order mean-motion resonance in terms of δ res .
For the libration width-based heuristic we classify as "plausibly resonant" planet pairs within five libration widths of a single mean-motion resonance.Because the δ res -based heuristic lacks a physically motivated threshold for classification as plausibly resonant, we classify as plausibly resonant planet pairs with δ res < 0.015 to roughly match the size of the libration width-based plausibly mean-motion resonant sample.We fully acknowledge that the planet pairs pairs we classify as plausibly resonant may not have librating orbital elements and therefore may not truly be in mean-motion resonances.Since our definition of plausibly resonant refers only to mean-motion resonances, from here the phrase plausibly resonant implies plausibly mean-motion resonant.
Because it is well-established that resonance overlap leads to chaotic behavior (e.g., Wisdom 1980;Lecar et al. 2001;Quillen 2011;Deck et al. 2013;Batygin & Morbidelli 2013b;Morrison & Kratter 2016), if a pair of planets is within five libration widths or δ res < 0.015 of multiple resonances then we do not consider the pair to be plausibly resonant.As the widths of resonances increase with pair total mass, the existence of planet pairs apparently occupying multiple resonances may imply that the mass-radius relation we used for the libration width-based definition overestimated some planet pairs' total masses.In support of this possibility, we find that planet pairs occupying multiple resonances according to the libration width-based definition are more massive than plausibly resonant pairs or pairs outside of resonance entirely.These planet pairs apparently in multiple resonances likely instead have lower masses than predicted by the Lissauer et al. ( 2011) mass-radius relation, indicating the existence of a large number of planets with 2 R ⊕ ≲ R p ≲ 4 R ⊕ that are less dense than Uranus and Neptune in our own solar system (e.g., Schlaufman & Halpern 2021).We also remove from our sample of plausibly resonant systems two-planet systems that dynamical analyses by Veras & Ford (2012) have shown must be circulating instead of librating under the assumption that the two observed planets are the only planets in the system.
For each planet pair in our sample, we report in Table 1 separations in units of libration widths from each first-order resonance we considered.In Table 2 we report the separations from the same first-order and additional second-order resonances in terms of |δ res |.We list in Table 3 the number of plausibly resonant pairs in each first-and second-order resonance we considered.We plot in Figure 2 stacked histograms comparing the period ratio distributions of our plausibly resonant and implausibly resonant samples.

Identification of Systems Plausibly Affected by Tidal Dissipation
Tidal dissipation is likely important for many of the Kepler-discovered systems of small planets that comprise the majority of both our plausibly resonant and implausibly resonant samples described in the previous subsection.To identify systems where tidal dissipation likely plays a significant role in the evolution of systems' period ratios, we calculate circularization times for the innermost planet in each system in our sample.The timescale for orbit circularization due to tidal dissipation inside a planet can be written where Q p , k 2 , M p , M * , a, R p , and P are the planet tidal quality factor, planet Love number, planet mass, stellar mass, semimajor axis, planet radius, and orbital period (e.g., Lee et al. 2013).Both Q p and k 2 are sensitively dependent on a planet's interior structure.While useful constraints on Q p and k 2 are available for solar system planets, the internal structures of most exoplanets are poorly constrained leaving exoplanet Q p and k 2 highly uncertain.
A solar metallicity star's luminosity L * scales like where The instellation incident on an exoplanet F p in units of Earth insolations F ⊕ is therefore In our sample, the median innermost planet instellation is approximately 200 F ⊕ .More than 90% of the innermost planets in our sample experience F p ≳ 30 F ⊕ .Fulton & Petigura (2018) showed that planets with R p ≲ 1.7 R ⊕ and F p ≳ 30 F ⊕ are likely rocky, as any gaseous envelopes initially present would have been stripped by photoevaporation (e.g., Owen & Wu 2017), core-powered mass loss (e.g., Ginzburg et al. 2018), or some other process.In our sample, innermost planets with R p ≲ 1.7 R ⊕ are likely rocky while innermost planets with R p ≳ 2.0 R ⊕ likely have at least some hydrogen/helium (H/He) atmospheres.
To evaluate the circularization timescale τ e for innermost planets in our sample with R p < 4 R ⊕ , we assume the dissipative properties of Mars.We favor Mars over Venus and Earth because its dissipative properties are both well constrained and unaffected by the presence of oceans or a massive atmosphere (e.g., Tyler 2021;Farhat et al. 2022).Based on detailed observations of the orbits of Phobos and Deimos, Duxbury & Callahan (1982) found 30 ≲ Q Mars ≲ 130.We use the smaller value to ensure that we have identified all systems with innermost planets plausibly affected by tidal evolution.We Note-The number of libration widths from each first-order resonance considered for each planet pair in the sample ordered by Kepler Input Catalog (KIC) identifier.We also provide indicators whether or a not planet pair is adjacent or plausibly resonant according to our libration-width based definition.This table is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.
confirm below the these tidal parameter values do not qualitatively affect our results.Yoder (1995) suggested k 2 = 0.14 for Mars.While planets with R p > 1.7 R ⊕ likely have H/He atmospheres, there are no planets with similar structures in the solar system and their dissipative properties are consequently poorly constrained.We therefore use Mars's properties for all innermost planets with R p < 4 R ⊕ that likely have a small fraction of their masses in H/He atmospheres.
To evaluate the circularization timescale for innermost planets in our sample with R p ≥ 4 R ⊕ , we assume the dissipative properties of Neptune.Among the solar system's gas and ice giant planets, Neptune is the most dissipative by at least an order of magnitude.The dynamics of Neptunian satellites suggest 9000 ≲ Q Neptune ≲ 36000 assuming a theoretically derived Neptunian k 2 = 0.41 (Bursa 1992;Zhang & Hamilton 2008).As before, we adopt a value at the low end of the range.While the assumption of similar dissipative properties for exoplanets with R p ≥ 4 R ⊕ might overestimate the important of tidal dissipation and lead to underestimated circularization times, we use Neptune's dissipative properties to ensure that we have identified all systems with innermost planets plausibly affected by tidal evolution.
We plot in Figure 3 circularization time as a function of orbital period for innermost planets in the Keplerdiscovered systems of multiple planets in our sample.Lee et al. (2013) showed that the eccentricities of planet pairs near first-order resonances decay due to tidal dissipation according to a shallow power law in time.
Those authors found that more than 50 circularization timescales were necessary to evolve a system initially in a mean-motion resonance to a period ratio 3% away from a rational number.Because most Kepler-discovered systems of multiple planets in our sample orbit mature solar-type stars, we assume that tidal dissipation can only meaningfully alter the near-resonant period ratios of a system if 50τ e ≲ 1 Gyr or equivalently if τ e < 20 Myr.

Relative Age Inferences
As argued above, the relative ages of Keplerdiscovered systems of small planets close to/far from period ratios suggestive of mean-motion resonances have the potential to both (1) decide whether departures from resonances occur relatively early or relatively late in the evolutionary histories of exoplanet systems and (2) evaluate the role of tidal evolution in moving systems out of resonances.If all multiple-planet systems are found to be the same age regardless of the period ratios of their constituent planets, then that observation would suggest that the features of the period ratio distribution of Kepler-discovered systems of small planets is set in the first Gyr of planetary system evolution.If multipleplanet systems with a plausibly resonant planet pair are younger than multiple-planet systems without a plausibly resonant planet pair, then planetary systems must evolve out of mean-motion resonances as they age.If an age offset is found and restricted to systems with small innermost planet circularization times, then that would strongly support tidal evolution as the explanation for Note-The absolute value of distance in δres from each first-and second-order resonance considered for each planet pair in the sample ordered by KIC identifier.We also provide indicators whether or a not planet pair is adjacent or plausibly first-/secondorder resonant according to our δres-based based definition.This table is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.
the features on the period ratio distribution of Keplerdiscovered systems of small planets.If an age offset is found in systems regardless of innermost planet circularization time, then that would suggest the importance of an additional non-tidal secular process.
Taking the approach pioneered in Schlaufman & Winn (2013) and Hamer & Schlaufman (2019, 2020, 2022), we use the Galactic velocity dispersions of these samples to compare their relative ages.Because the velocity dispersion of a thin disk stellar population grows with the average age of that stellar population (e.g., Binney et al. 2000), populations with significantly different velocity dispersions have significantly different ages.The relationship between age and velocity dispersion in the thin disk of the Milky Way is a function of the specific region of the Galaxy.As a result, calibrations like that of Binney et al. (2000) based on solar neighborhood samples cannot be directly applied in the Kepler field.While it is possible to calibrate the relationship between velocity dispersion and age in the Kepler field, virtually all stellar age inference techniques ultimately rely on stellar models.They are therefore model dependent.On the other hand, the Hamer & Schlaufman (2019, 2020, 2022) approach enables precise, model-independent relative age comparisons.Conclusions based on model-independent relative ages avoid the systematic uncertainties associated with other age inference methodologies.We therefore prefer these uncalibrated but model-independent relative age inferences for the analyses that follow.For a particular sample of stars, we first convert astrometric and radial velocity data into Galactic U V W space velocities.For each star, we next use pyia to generate 100 random sets of equatorial coordinates, parallaxes, and proper motions in accord with Gaia astrometric solution covariance matrix for that star.We then use astropy to generate 100 random radial velocities in accord with each star's observed radial velocity and its uncertainty (Astropy Collaboration et al. 2013, 2018;Price-Whelan & Brammer 2021).We next calculate 100 sets of U V W velocities for each star and then infer 100 Galactic velocity dispersions for a particular sample according to where U , V , and W are the components of the Galactocentric velocity and bars denote median values.We finally take the median Galactic velocity dispersion across these 100 realizations as a particular sample's Galactic velocity dispersion.The typical U , V , and W velocity uncertainties we infer in this way are less than about 0.2 km s −1 , much smaller than the typical 1-σ ranges of our inferred velocity dispersion distributions.Individual stellar U V W velocity uncertainties are therefore much too small to impact our analyses.

Plausibly First-order Resonant Systems
We first compare the Galactic velocity dispersions and therefore relative ages of Kepler-discovered systems of small planets with/without plausibly first-order res-onant planet pairs.This sample of systems without plausibly first-order resonant planet pairs does include systems with plausibly second-order resonant planet pairs, though we have confirmed that our conclusions are unchanged if we also remove systems with plausibly second-order resonant planet pairs.We use 5000 bootstrap samples to characterize the underlying Galactic velocity dispersion distributions of each subsample from the single realization we observe and plot in Figure 4 the resulting velocity dispersion distributions assuming both the libration width-based and δ res -based definitions for plausibly first-order resonant systems.We find that systems with plausibly first-order resonant planet pairs have a Galactic velocity dispersion consistent with the Galactic velocity dispersion of systems that lack such a pair.The implication is that the age distribution of systems with plausibly first-order resonant planet pairs overlaps with the age distribution of systems that lack such a pair.
We then divide our sample into two subsamples: one in which tidal dissipation in a system's innermost planet is likely important (i.e., τ e < 20 Myr) and one in which such dissipation is unlikely to be important (i.e., τ e ≥ 20 Myr).We apply the same bootstrap resampling strategy described in the preceding paragraph and plot in Figure 5 the resulting velocity dispersion distributions.For systems with an innermost planet likely affected by tides, we find that systems with a plausibly first-order resonant planet pair have a colder Galactic velocity dispersion than systems that lack a plausibly resonant planet .pair.This is so for both the libration width-based and δ res -based definitions for plausibly first-order resonant systems.For the libration width-based definition, the probability that an offset this large could be produced by chance is about 1 in 72 or equivalently about 2.2 σ.For the δ res -based definition, the probability that an offset this large could be produced by chance is about 1 in 1500 or equivalently about 3.1 σ.These observations suggest that that tidal dissipation plays an important role in moving away from first-order resonances multiple-planet systems left in first-order resonances after the dissipation of their parent protoplanetary disks.
It has been proposed in certain models accounting for secular interactions in compact multiple-planet systems (e.g., Hansen & Murray 2015) or obliquity tides (e.g., Millholland & Laughlin 2019;Millholland & Batygin 2019;Millholland 2019) that the efficiency of tidal dissipation inside an innermost planet may exceed our assumed value.To investigate this possibility, we divided our sample based on its median innermost planet circularization time τ e = 400 Myr.Dividing the sample at its median value allows for a qualitative examination of the role of tides independent of a particular model while also equalizing the sizes of the tidal/non-tidal samples and maximizing the statistical power of their comparison.We report the outcome of this exercise in Figure 6.Regardless of the definition used to identify plausibly resonant planet pairs, systems with circularization times τ e < τ e and a plausibly first-order resonant planet pair have colder Galactic velocity dispersion than sys- tems that lack such a pair.The probabilities that offsets in the tidally affected samples this large could occur by chance are about 1 in 18 (equivalent to about 1.6 σ) for the libration-width definition and 1 in 93 (equivalent to about 2.3 σ) for the δ res -based definition.On the other hand, systems with circularization times τ e > τ e do not display an offset in Galactic velocity dispersion between plausibly first-order resonant/implausibly firstorder resonant systems.The similarity between our results using two different definitions for the importance of tidal evolution shows that our results are not sensitively dependent on the choice of tidal parameters.

Plausibly Second-order Resonant Systems
We next compare the Galactic velocity dispersions and therefore relative ages of Kepler-discovered systems of small planets with/without plausibly second-order resonant planet pairs using the same procedure described in the preceding subsection.This sample of systems without plausibly second-order resonant planet pairs does in-clude systems with plausibly first-order resonant planet pairs, though we have confirmed that our conclusions are unchanged if we also remove systems with plausibly first-order resonant planet pairs.The libration-width based definition is only applicable for plausibly firstorder resonant systems, so we only use the δ res -based definition for plausibly second-order resonant systems in this subsection.We plot in Figure 7 the resulting velocity dispersion distributions under the δ res -based definition for plausibly second-order resonant systems.We find that systems with plausibly second-order resonant planet pairs have a significantly colder Galactic velocity dispersion than systems that lack such a pair.The implication is that systems with plausibly second-order resonant planet pairs are younger than systems that lack such a pair.The probability that an offset this large could occur by chance is less than 1 in 100000, or equivalently about 4.6 σ.
To investigate the possible importance of tidal evolution in moving systems away from second-order reso-  The Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems with a plausibly first-order resonant planet pair is consistent with the Galactic velocity dispersion of systems that lack such a pair regardless of the heuristic used to identify systems with a plausibly first-order resonant planet pair.
nance, as before we split our sample of systems into two subsamples: one in which tidal dissipation inside the innermost planet is likely important (i.e., τ e < 20 Myr or τ e < τ e = 400 Myr) and one in which tides are unlikely to be important (i.e., τ e ≥ 20 Myr or τ e ≥ τ e = 400 Myr).We apply the same bootstrap resampling strategy described in the preceding subsection and plot in Figure 8 the resulting velocity dispersion distributions.We find that systems with a plausibly second-order resonant planet pair have a significantly colder Galactic velocity dispersion than systems that lack such a pair regardless of the possible impact of tidal dissipation.When splitting our sample at τ e = 20 Myr, the probabilities of offsets as large as we observe by change are about 1 in 210 (equivalent to about 2.6 σ) and about 1 in 76000 (equivalent to about 4.2 σ) for the systems affected/unaffected by tides.When splitting our sample at τ e = 400 Myr, the probabilities of offsets as large as we observe by change are about 1 in 280 (equivalent to about 2.7 σ) and about 1 in 76000 (equivalent to about 4.2 σ) for the systems affected/unaffected by tides.
For systems with a plausibly second-order resonant planet pair, the lack of a dependence of velocity dispersion offset on the exact criteria used to identify systems affected by tidal evolution suggests that tidal dissipation is not the best candidate to explain our observation that systems with a plausibly second-order resonant planet pair are younger than systems that lack such a pair.As we will argue in the next section, some other non-tidal secular processes is likely responsible for drawing systems out of second-order resonances.

DISCUSSION
The Galactic velocity dispersion of the sample of Kepler-discovered multiple-planet systems with at least one plausibly second-order resonant planet pair is colder We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top panels include systems in which the innermost planet is likely affected by tidal evolution (i.e., τe < 20 Myr), while bottom panels include systems in which the innermost planet is unlikely to be affected by tidal dissipation (i.e., τe ≥ 20 Myr).Left: plausibly resonant/implausibly resonant according to the libration width-based definition.Right: plausibly resonant/implausibly resonant according to the δres-based definition.The Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems likely affected by tides and with a plausibly first-order resonant planet pair is colder than the Galactic velocity dispersion of systems that lack such a pair, regardless of the definition of plausibly first-order resonant systems.The probabilities that offsets this large could occur by chance is about 1 in 72 (equivalent to about 2.2 σ) for the libration width-based definition and about 1 in 1500 (equivalent to about 3.1 σ) for the δres-based definition.The large and statistically significant offsets present in the likely tidally affected samples support the idea that tidal evolution is responsible for moving planetary systems formed in first-order resonance away from those resonances over Gyr timescales.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top panels include systems in which the innermost planet is likely affected by tidal evolution (i.e., τe < τ e = 400 Myr), while bottom panels include systems in which the innermost planet is unlikely to be affected by tidal dissipation (i.e., τe ≥ τ e = 400 Myr).Left: plausibly resonant/implausibly resonant according to the libration width-based definition.Right: plausibly resonant/implausibly resonant according to the δres-based definition.Regardless of the definition used to identify plausibly first-order resonant systems, systems with innermost planet circularization times less than the median value τ e = 400 Myr have a colder Galactic velocity dispersion and are therefore younger than implausibly resonant systems.The probabilities that offsets in the tidally affected samples this large could occur by chance are about 1 in 18 (equivalent to about 1.6 σ) for the libration-width definition and 1 in 93 (equivalent to about 2.3 σ) for the δres-based definition.Our conclusion that plausibly first-order resonant systems with an innermost planet likely to be affected by tidal dissipation are young does not sensitively depend on the criteria used to identify systems likely affected by tides.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.The Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems with a plausibly second-order resonant planet pair is significantly colder than the Galactic velocity dispersion of systems that lack such a pair.The probability that an offset this large could occur by chance is less than 1 in 100000, or equivalently about 4.6 σ.This observation implies that systems initially in second-order resonance are driven out of resonance over Gyr timescales.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top panels include systems in which the innermost planet is likely affected by tidal evolution, while bottom panels include systems in which the innermost planet is unlikely to be affected by tidal dissipation.Left: Tidal and non-tidal samples separated at τe = 20 Myr.The probabilities that offsets this large could occur by chance is about 1 in 210 (equivalent to about 2.6 σ) for the systems likely affected by tides and about 1 in 76000 (equivalent to about 4.2 σ) for systems unlikely to be unaffected by tides.Right: Tidal and non-tidal samples separated by the median circularization time τe = 400 Myr.The probabilities that offsets this large could occur by chance is about 1 in 280 (equivalent to about 2.7 σ) for the systems likely affected by tides and about 1 in 76000 (equivalent to about 4.2 σ) for systems unlikely to be unaffected by tides.Regardless of the importance of tidal dissipation inside innermost planets, the Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems with a plausibly second-order resonant planet pair is colder than the Galactic velocity dispersion of systems that lack such a pair.The implication is that tidal dissipation is not the process responsible for moving systems initially in second-order mean-motion resonance out of resonance.
than the Galactic velocity dispersion of the sample of Kepler-discovered multiple-planet systems without such a pair.The implication is that systems with a plausibly second-order resonant pair are systematically young.The same is true for Kepler-discovered multiple-planet systems with at least one plausibly first-order resonant planet pair, but only when the innermost planet in the system is close enough to its host that tidal dissipation inside the planet likely affects the secular evolution of the system.

Evaluation of Potentially Observable Consequences of the Tidal Evolution Implied by Our Observations
Because tides appear to be responsible for the Galactic velocity dispersion and therefore age offsets we observe between first-order resonant and non-resonant systems, we now assess whether or not the tidal heating of innermost planets in first-order resonant systems with forced eccentricities would be observable.Following Jackson et al. (2008) the tidal heating power experienced by a planet can be written where G is the gravitational constant, M * its host star's mass, R p its radius, Q ′ p its modified tidal quality factor, a its semimajor axis, and e its orbital eccentricity.For the sample of Kepler-discovered multiple-planet systems we identified as plausibly first-order resonant, nonzero innermost planet orbital eccentricities and Equation (11) can be used to determine the tidal heating experienced by innermost planets.We use the same Q ′ p values we assumed in out circularization timescale calculations and assume e = 0.02 corresponding to the median libration width-minimizing eccentricity in our sample.We find a typical tidal heating power H ∼ 10 15 W and a maximum H ≳ 10 17 W.
While important for the internal structures of rocky planets, the observable effects of even H ≳ 10 17 W of tidal heating are unlikely to be observable.Valencia et al. (2007) showed that 10 17 W/M ⊕ of tidal heating can cause the partial melting of a planet's lower lower mantle, and we identify four planets in our sample experiencing tidal heating exceeding that threshold: Kepler-1371 c/KOI-2859.02(H ≈ 2.6 × 10 17 W), Kepler-1669 c/KOI-1475.01 (H ≈ 3.2 × 10 17 W), Kepler-9 d/KOI-377.03(H ≈ 3.7 × 10 17 W), and Kepler-342 e/KOI-1955.03(H ≈ 5.2 × 10 17 W).Valencia et al. (2007) also showed that H ≈ 6.8 × 10 17 W of tidal heating in GJ 876 d would only increase its radius by about 100 km.As a result, in the four systems listed above the observational consequences of the tidal heating are much smaller than the observational consequences of uncertainties in mass-radius relations due to planet composition.They are therefore unobservable as larger-thanexpected planet radii.
These tidal heating powers are also unlikely to be observable in the form of higher-than-expected planet equilibrium temperatures.For a possibly measurable effect, the tidal heating experienced by a planet would need to be comparable to its instellation flux multiplied by its cross-sectional area and one minus its Bond albedo.Only planets with small orbital separations experience significant tidal heating, and these plants also experience intense instellations that are much more important for setting their equilibrium temperatures than tidal heating.On the other hand, planets in high planetary obliquity states affected by obliquity tides may experience tidal heating with sufficient intensities to increase their equilibrium temperatures (e.g., Millholland & Laughlin 2019).An observation that revealed an enhanced planetary equilibrium temperature would be hard to confidently attribute to obliquity tides though, as apparently increased equilibrium temperatures could be caused by an overestimated planetary Bond albedo or by the impact of an atmospheric greenhouse effect.A possibly more easily observable consequence of intense tidal heating from obliquity tides would by inflated radii of planets with significant hydrogen/helium atmospheres wide of mean-motion resonances (e.g., Millholland & Laughlin 2019;Millholland 2019).

Consistency Between our Results and Published
Explanations for the Period Ratio Distribution Observed by Kepler Our observation of a velocity dispersion and therefore age offset between our samples of plausibly secondorder resonant and implausibly resonant systems unlikely to be affected by tidal evolution suggests that some other process(es) beyond tidal dissipation may be necessary to explain our observation.While our observation has no implications for the absolute ages of plausibly second-order resonant/implausibly resonant systems, the Galactic velocity dispersion offset we observe implies a relative age difference between plausibly second-order resonant/implausibly resonant systems of at least 100 Myr.To investigate the possibility that plausibly resonant systems are a young population with age ≲ 1 Gyr and implausibly resonant systems are an old population with age ≳ 1 Gyr, we examine the distribution of lithium equivalent width and abundance as well as Ca II H and K log R ′ HK and S HK activity indices in both populations using data from Berger et al. (2018) and Brewer & Fischer (2018).We plot these dis-tributions in Figure 9. Neither population displays obvious signs of youth and plausibly resonant systems do not have obviously larger observed lithium equivalent widths, inferred lithium abundances, or log R ′ HK /S HK activity indices.The implication is that neither population has a typical age ≲ 1 Gyr.
Because the host stars of both plausibly second-order resonant and implausibly resonant systems have lithium properties and activity indices suggesting ages in excess of 100 Myr, we assert that it is unlikely that the age difference between the two populations we identify is imparted by a process driving second-order resonant planet pairs out of resonance within the first 100 Myr of a system's evolution.Indeed, if initially resonant systems had already been driven away from resonance by the time of their parent protoplanetary disks' dissipations, then the age difference between plausibly resonant/implausibly resonant systems would be small and significant Galactic velocity dispersion offsets between populations would not exist.We therefore argue that in situ formation of planets (e.g., Petrovich et al. 2013), planetary interactions with protoplanetary disk density waves raised by additional planets (e.g., Podlewska-Gaca et al. 2012;Baruteau & Papaloizou 2013;Cui et al. 2021), planetplanet interactions in the presence of protoplanetary disk-mediated eccentricity damping (e.g., Charalambous et al. 2022;Laune et al. 2022), and system disruptions due to progressive protoplanetary disk dissipation (e.g., Liu et al. 2022) are all unlikely to explain the multipleplanet system period ratio distribution observed by Kepler.
Because most Kepler-discovered multiple-planet systems orbit solar-type stars that spin down on the order of 100 Myr, we also assert that evolving host star quadrupolar potentials are unlikely to be responsible for disrupting resonant systems and generating the age offset we inferred.Stellar oblateness due to rotation can cause small separation planetary systems to depart from idealized Keplerian orbits.In these cases, the quadrupole term of the stellar potential that scales like J 2 R 2 * is usually dominant (e.g., Schultz et al. 2021).The stellar quadrupole moment J 2 can be approximated as where k 2 is related to the apsidal motion constant k 2 /2 and Ω * , M * , and R * are the stellar angular velocity, mass, and radius (e.g., Schultz et al. 2021).We use the stellar models of Ekström et al. (2012) that account for the rotational evolution of stars to evaluate the time evolution of J 2 R 2 * .We take k 2 from Batygin & Adams (2013) assuming apsidal motion constants k 2 /2 (e.g., Csizmadia et al. 2019) as they do for fully convective (k 2 /2 = 0.14) and fully radiative stars (k 2 /2 = 0.01).Because a solar mass and composition main sequence star has a radiative core and a convective envelope, we assume k 2 = 0.14 for the limiting case where evolution of the quadrupole term of the stellar potential will be most important.We find that J 2 R 2 * evolves by an order of magnitude during the first 100 Myr, by a factor of several during next few hundred Myr, and by very little after 1 Gyr for the duration of the main sequence.We therefore argue that evolving stellar quadrupolar potentials are unlikely to explain the multiple-planet system period ratio distribution observed by Kepler.
It has been shown that initially resonant multipleplanet systems can become unstable shortly after the dissipations of their parent protoplanetary disks with instability times inversely proportional to the number of planets in a system and their characteristic mass (e.g., Matsumoto et al. 2012;Izidoro et al. 2017;Pichierri et al. 2018;Pichierri & Morbidelli 2020;Izidoro et al. 2021;Goldberg & Batygin 2022;Goldberg et al. 2022;Rice & Steffen 2023).Instabilities can also arise in initially stable systems as a result of either stellar or planetary mass loss.Tightly packed resonant systems close to the threshold for dynamical instability can become unstable if host stars loses just 1% of their masses or if systems' total planetary masses change by more than 10% (e.g., Matsumoto & Ogihara 2020).Larger mass losses can produce instabilities even in systems initially far from the threshold for dynamical instability.Wood et al. (2002) derived a relationship between observed X-ray fluxes and stellar mass loss rates and uses that relationship to predict the mass loss rates of solar-type stars as function of stellar age.They made the provocative suggestion that the Sun may have lost up to 3% of its mass by 1 Gyr.Even if this extreme estimate is accurate, the vast majority of that mass loss takes place in the first 100 Myr of a system's evolution.That timescale is too short to explain our observations.Late time dynamical instabilities in multiple-planet systems driven by stellar mass loss are therefore unlikely to be responsible for disrupting resonant systems and generating the age offset we inferred.
Late time dynamical instabilities in multiple-planet systems driven by planetary mass loss may play some role in the disruption of initially resonant systems and the production of the age offsets we inferred though.To evaluate this scenario for the typical system in our sample, we use the models of Lopez & Fortney (2013).The typical planet in our sample has M p ≈ 5 M ⊕ and experiences an instellation F p ∼ 100 F ⊕ .Assuming the Lopez & Fortney (2013) 2018) lithium abundance as a function of GBP − GRP.At constant GBP − GRP, large equivalent widths or lithium abundances should be observed in the spectra of stars younger than about 650 Myr.We see no significant equivalent width or abundance offsets between plausibly resonant/implausibly resonant systems, suggesting that these plausibly resonant systems are not young in an absolute sense.Bottom left: from Brewer & Fischer (2018) log R ′ HK as a function of GBP − GRP.Bottom right: from Brewer & Fischer (2018) SHK as a function of GBP − GRP.log R ′ HK and SHK are measures of the strength of emission in the core of the calcium H & K lines.At constant GBP − GRP, stronger emission should be observed in the spectra of younger stars.We see no significant log R ′ HK or SHK offsets between plausibly resonant/implausibly resonant systems, suggesting that these plausibly resonant systems are not young in an absolute sense.equal to 0.1, the typical planet in our sample that initially possessed a massive hydrogen/helium (H/He) atmosphere could have lost more than 80% of its mass in H/He in just 50 Myr.A planet that originally had just 10% of its mass in H/He would lose its entire atmosphere in 100 Myr.Both timescales are too short explain our result.On the other hand, a planet with 10% of its original mass in H/He experiencing F p ≈ 50 F ⊕ would lose 50% of this mass in 100 Myr and the rest in 1 Gyr.At face value then, dynamical instabilities due to the loss of primordial H/He atmospheres in initially resonant systems experiencing intermediate instellations could play some role in the explanation of our result.Dynamical instabilities caused by planetary mass loss would result in collisions that leave behind systems with large mutual inclinations and high eccentricities inconsistent with the small mutual inclinations and low eccentricities inferred for Kepler-discovered system of multiple planets (e.g., Wu & Lithwick 2013;Fabrycky et al. 2014;Van Eylen & Albrecht 2015;Mills et al. 2019;He et al. 2020).
Dynamical instabilities do represent a possible solution to the so called "Kepler dichotomy" though (e.g., Johansen et al. 2012;Izidoro et al. 2017Izidoro et al. , 2021)).In this scenario, single-planet systems may have initially been multiple-planet, possibly resonant systems that experienced a dynamical instability that excited eccentricities and inclinations.If this scenario is accurate, then single-planet systems should be older than multiple-planet systems.To evaluate this proposed explanation for the Kepler dichotomy, we use the Thompson et al. (2018) Kepler DR25 KOI list and compare the Galactic velocity dispersions of the populations of single-transiting systems, plausibly first-order/secondorder resonant multiple-planet systems, implausibly resonant multiple-planet systems, and all multiple-planet systems.We show the results of these calculations in Figure 10.The population of single-transiting systems has a Galactic velocity dispersion significantly warmer than the population of plausibly second-order resonant multiple-planet systems.The Galactic velocity dispersions of single-transiting, first-order resonant, implausibly resonant multiple-planet, and all multiple-planet systems are all consistent.Our results contrast with LAMOST-based inferences that found Galactic velocity dispersion decreases with planet multiplicity (Chen et al. 2021a,b).We find no systematic differences between single-transiting and multipletransiting systems in lithium equivalent width, lithium abundance, log R ′ HK , or S HK .In accord with the idea that single-transiting systems are the outcomes of dynamical instabilities, these analyses suggests that plausibly second-order resonant systems are younger than single-transiting systems.Plausibly second-order resonant systems aside, our analyses provide no evidence that single-transiting systems are older than multipletransiting systems in general.

Interactions with Drifting Planetesimals as a Possible Explanation
Interactions with a residual disk of planetesimals are a mechanism that can explain the observations that planet pairs tend to favor period ratios just wide of resonances and tend to avoid period ratios just interior to resonances.Chatterjee & Ford (2015) and Ghosh & Chatterjee (2023) have simulated the interaction of resonant planet pairs with co-located planetesimal disks remaining after the dissipation of gaseous protoplanetary disks.Those analyses have shown that planets interacting with a total mass of planetesimals at least 1% of the total planetary mass can produce departures from resonances in less than 1 Myr similar in magnitude to those observed by Kepler.While that timescale is too short to explain our results, the arrival of a net planetesimal mass 1% of the total planetary mass over a 1 Gyr interval could in principle reproduce both the Kepler period-ratio distribution and our observation that plausibly second-order resonant systems have a colder Galactic velocity dispersion than implausibly resonant systems.
A delayed interaction with an exterior disk of residual planetesimals is a key ingredient for explaining the architecture of the solar system in the Nice model (e.g., Gomes et al. 2005;Morbidelli et al. 2005;Tsiganis et al. 2005).It has also been suggested that it would take a few hundred million years for the terrestrial planets to clear the inner solar system of small bodies and a few billion years for the gas & ice giants to clear the outer solar system of small bodies (e.g., Goldreich et al. 2004).Simulations of interactions between Jupiter, Saturn, and a residual planetesimal disk have shown that the timescale for the initiation of planetesimal-driven migration is sensitive to the location of the inner edge of the residual planetesimal disk, and can take more than 1 Gyr (e.g., Gomes et al. 2005;Thommes et al. 2008).
While Chatterjee & Ford (2015) and Ghosh & Chatterjee (2023) considered instantaneous interactions with co-located planetesimals, extended interactions with migrating planetesimals could also take place.Planetesimals can migrate long distances over billions of years due to a process called the Yarkovsky effect in which a rotating planetesimal heated on its day side by its host star will experience a non-gravitational acceleration as its hotter day side radiates more than its cooler night side (prograde rotators migrate outwards and retrograde rotators migrate inwards).The net force experienced by a planetesimal from the Yarkovsky effect has complex dependencies on its thermal parameters, albedo, emissivity, size, rotation rate, and obliquity.Yarkovsky effect-driven accelerations can cause significant semimajor axis evolution for planetesimals with sizes from meters to tens of kilometers.
To evaluate the potential of the Yarkovsky effect to deliver planetesimal mass into the vicinities of Keplerobserved multiple-planet system over Gyr timescales, we first roughly estimate the total mass in planetesimals necessary to deliver via the Yarkovsky effect planetesimals amounting to about 1% of the typical total planetary mass in our sample of plausibly second-order resonant systems.Assuming the mass-radius relation from Lissauer et al. (2011), the typical total planet mass in our sample of plausibly second-order resonant systems is M p,tot ≈ 10 M ⊕ .Assuming planetesimal properties like those in the solar system, the Yarkovsky effect most strongly affects planetesimals with sizes R in the range 1 km ≲ R ≲ 10 km.To roughly calculate the total disk mass necessary to provide M tot ≈ 0.1 M ⊕ of migrating 1 or 10 km planetesimals, we assume an equal number of prograde/retrograde rotators and first assume the planetesimal size distribution produced by detailed simulations of the streaming instability dN/dR ∝ R −2.8 (Simon et al. 2016  We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top: the populations of plausibly first-order (blue) and second-order (orange) resonant multiple-planet systems.Top-middle: the population of implausibly resonant multiple-planet systems.Bottom-middle: the population of all multiple-planet systems.Bottom: the populations of confirmed (orange) and candidate (blue) single-transiting systems.The plausibly second-order resonant system population has a significantly colder Galactic velocity dispersion than the single-transiting system population.The population of single-transiting systems has a Galactic velocity dispersion consistent with the populations of implausibly resonant multipleplanet systems and all multiple-planet systems.
inate the migrating planetesimal mass budget, then we find that a total disk mass M tot ≈ (32, 20) M ⊕ would be necessary to supply sufficient planetesimal mass to affect period ratios.We next assume the planetesimal size distribution appropriate for steady-state collisional evolution dN/dR ∝ R −3.5 (Wyatt 2008).If objects with R = (1, 10) km dominate the migrating planetesimal mass budget, then we find that a total disk mass M tot ≈ (3, 9) M ⊕ would be necessary to supply suffi-cient planetesimal mass to affect period ratios.Because we assumed that only planetesimals with a single size migrate, these are upper limits on the total disk mass.
We use the simplified model of the Yarkovsky effect presented in Veras et al. (2019) and implemented in REBOUNDx 3.1.0(Tamayo et al. 2020) to infer an upper limit on the amount of migration experienced by planetesimals as a function of R. We first assume that the semimajor axes of planetesimals migrating inwards shrink until they reach the exterior 2:1 mean-motion resonance with the outermost planet in an initially secondorder resonant system.Upon reaching that location, an inwardly migrating planetesimal will experience eccentricity excitation leading to orbit crossing and eventual collision with a planet in the system.In our sample of plausibly second-order resonant systems, the typical semimajor axis of the outermost observed planet is a = 0.15 AU.The corresponding exterior 2:1 meanmotion resonance would be located at a 1 = 0.24 AU.We assume a planetesimal density ρ plan = 2.7 kg m −3 and a Bond albedo A = 0.017.The median stellar mass in our sample of plausibly second-order resonant systems M * = 0.92 M ⊙ corresponding to a luminosity L * ≈ 0.72 L ⊙ .We evolve the orbits of planetesimals with 1 km and 10 km for 0.1 Myr at an array of semimajor axes values to estimate the instantaneous da/dt at many a.At time t = 0 we place a planetesimal at the location of the exterior 2:1 mean-motion resonance for our median system and integrate its orbit backward in time for 1 Gyr.We find that the migration as a result of the Yarkovsky effect can move (1,10) km planetesimals to a 1 = 0.23 AU from a 0 = (0.45, 0.26) AU in 500 Myr and from a 0 = (0.62, 0.29) in 1 Gyr.Our inferred total planetesimal disk masses and migration rates imply surface densities Σ plan in the range 10 2 g cm −3 ≲ Σ plan ≲ 10 3 g cm −3 consistent with the solid surface density inferred for the minimum-mass solar nebula 10 1 g cm −3 ≲ Σ MMSN,solid ≲ 10 2 g cm −3 (e.g., Weidenschilling 1977;Asplund et al. 2009).
While the total amount of solid mass and the solid surface density required for Yarkovsky effect-driven planetesimal migration are plausible, maintaining enough mass in (1,10) km planetesimals for 1 Gyr could be challenging.Planetesimal disk mass inferences based on observable dust properties suggest that planetesimal disks can contain M tot ∼ 100 M ⊕ shortly after protoplanetary disk dissipation (e.g., Wyatt 2008).Those masses are expected to drop below M tot ∼ 1 M ⊕ after about 10 Myr due to collisional grinding.At the same time, apparent episodes of late accretion in the solar system provides evidence for planetesimals disks lasting for a few hundred Myr.Indeed, a residual disk of planetesimals with 10 M ⊕ ≲ M tot ≲ 100 M ⊕ in the outer solar system is thought to be necessary to explain the orbital properties of the gas and ice giants (e.g., Hahn & Malhotra 1999).Likewise, a residual disk of planetesimals in the inner solar system with M tot ∼ 1 M ⊕ is thought to be necessary to damp the eccentricities and inclinations of the terrestrial planets after the giant impact phase of planet formation a few tens of Myr after the formation of calcium-aluminum rich inclusions at the dawn of the solar system (e.g., Schlichting et al. 2012).While more detailed simulations will be necessary to confirm or reject the scenario outlined above, these solar systembased inferences indicate that our scenario is at least a plausible explanation for our observations.

CONCLUSION
We find that while plausibly second-order meanmotion resonant multiple-planet systems discovered by the Kepler Space Telescope are not young in an absolute sense, they do have a colder Galactic velocity dispersion and are therefore younger than both Keplerdiscovered implausibly resonant multiple-planet systems and single-transiting systems.The same is true for firstorder mean-motion resonant systems, but only when tidal dissipation inside a system's innermost planet is likely important for the system's secular evolution.Our observations are inconsistent with any proposed physical mechanisms that drive initially resonant systems away from resonance in a gaseous protoplanetary disk, during its dissipation, or that take place in the first 100 Myr of a system's evolution.
We confirm that the age offset between plausibly second-order resonant and implausibly resonant systems persists at high statistical significance even when tidal dissipation inside the innermost planet is likely too weak to act as an angular momentum sink.The implication is either that tides are more efficient than expected (perhaps due to obliquity tides) or that an additional nontidal secular process that acts over hundreds of Myr to one Gyr is responsible for moving initially second-order resonant systems out of resonance while leaving them with small eccentricities and mutual inclinations.We propose that interactions with km-sized planetesimals migrating inward due to the Yarkovsky effect from a residual planetesimal belt at a ∼ 1 AU with total mass 1 M ⊕ ≲ M tot ≲ 10 M ⊕ are a plausible explanation consistent with our observations and worthy of more careful future study.Based on our observation that plausibly second-order resonant systems are younger than both implausibly resonant multiple-planet and singletransiting systems, we predict that plausibly secondorder resonant systems should be frequently found orbiting young stars with ages between about 10 Myr and 100 Myr.The same should be true for first-order resonant systems with an innermost planet likely affected by tidal dissipation.Both plausibly first-and second-order resonant systems should then become less common as systems age.
2009415.We thank Jason Steffen for helpful comments and for providing their catalog of Kepler multipleplanet systems in advance of publication.We also thank the anonymous referee for an insightful suggestion that significantly improved this article.This research has made use of the NASA Exoplanet Archive (Akeson et al. 2013), which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000).This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (10.26093/cds/vizier).The original description of the VizieR service was published in Ochsenbein et al. (2000).Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.  2013,2018), galpy (Bovy 2015), matplotlib (Hunter 2007), pandas (McKinney et al. 2010), pyia (Price-Whelan & Brammer 2021), REBOUNDx (Tamayo et al. 2020)

Figure 1 .
Figure 1.Distributions of separation from closest mean-motion resonance.In each panel the blue histogram gives the number of adjacent planet pairs in each histogram bin, while the orange histogram stacked on top of the blue histogram gives the number of non-adjacent planet pairs in each histogram bin.Top: distribution of separation from first-order mean-motion resonances in units of libration widths.We consider planet pairs within five resonance libration widths "plausibly resonant" in our subsequent analyses.Middle: distribution of separation from first-order mean-motion resonances in δres.Bottom: distribution of separation from second-order mean-motion resonances in δres.We consider planet pairs with δres < 0.015 plausibly resonant in our subsequent analyses.Since our definition of plausibly resonant refers only to mean-motion resonances, from here the phrase plausibly resonant implies plausibly mean-motion resonant.

Figure 2 .
Figure 2. Distribution of planet-pair period ratios.The orange histogram gives the number of implausibly resonant pairs in each period ratio bin while the blue histogram stacked on top of the orange histogram gives the number of plausibly resonant pairs in each period ratio bin.Top: plausibly resonant/implausibly resonant according to the libration width-based definition.Plausibly resonant pairs far from period ratios suggestive of resonances are high total-mass planet pairs with larger libration widths.Bottom: plausibly resonant/implausibly resonant according to the δres-based definition.

Figure 3 .
Figure3.Circularization time as a function of orbital period for the innermost planet in each multiple-planet system in our sample.We plot as blue points planets for which we assume the dissipative properties of Mars Qp = 30 & k2 = 0.14(Duxbury & Callahan 1982;Yoder 1995) and as orange points planets for which we assume the dissipative properties of Neptune Qp = 9000 & k2 = 0.41(Bursa 1992;Zhang & Hamilton 2008).A point's size encodes the mass of the planet it represents calculated using theLissauer et al. (2011) mass-radius relation.According toLee et al. (2013), tidal evolution can measurably alter a system's period ratio distribution when the circularization time of its innermost planet τe ≲ 20 Myr.

Figure 4 .
Figure 4. Distributions of Galactic velocity dispersion based on bootstrap resampling for systems with (blue) or without (orange) a plausibly first-order resonant planet pair.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top: plausibly resonant/implausibly resonant according to the libration width-based definition.Bottom: plausibly resonant/implausibly resonant according to the δres-based definition.The Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems with a plausibly first-order resonant planet pair is consistent with the Galactic velocity dispersion of systems that lack such a pair regardless of the heuristic used to identify systems with a plausibly first-order resonant planet pair.

Figure 5 .
Figure5.Distributions of Galactic velocity dispersion based on bootstrap resampling for systems with (blue) or without (orange) a plausibly first-order resonant planet pair.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top panels include systems in which the innermost planet is likely affected by tidal evolution (i.e., τe < 20 Myr), while bottom panels include systems in which the innermost planet is unlikely to be affected by tidal dissipation (i.e., τe ≥ 20 Myr).Left: plausibly resonant/implausibly resonant according to the libration width-based definition.Right: plausibly resonant/implausibly resonant according to the δres-based definition.The Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems likely affected by tides and with a plausibly first-order resonant planet pair is colder than the Galactic velocity dispersion of systems that lack such a pair, regardless of the definition of plausibly first-order resonant systems.The probabilities that offsets this large could occur by chance is about 1 in 72 (equivalent to about 2.2 σ) for the libration width-based definition and about 1 in 1500 (equivalent to about 3.1 σ) for the δres-based definition.The large and statistically significant offsets present in the likely tidally affected samples support the idea that tidal evolution is responsible for moving planetary systems formed in first-order resonance away from those resonances over Gyr timescales.

eeFigure 6 .
Figure6.Distributions of Galactic velocity dispersion based on bootstrap resampling for systems with (blue) or without (orange) a plausibly first-order resonant planet pair.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top panels include systems in which the innermost planet is likely affected by tidal evolution (i.e., τe < τ e = 400 Myr), while bottom panels include systems in which the innermost planet is unlikely to be affected by tidal dissipation (i.e., τe ≥ τ e = 400 Myr).Left: plausibly resonant/implausibly resonant according to the libration width-based definition.Right: plausibly resonant/implausibly resonant according to the δres-based definition.Regardless of the definition used to identify plausibly first-order resonant systems, systems with innermost planet circularization times less than the median value τ e = 400 Myr have a colder Galactic velocity dispersion and are therefore younger than implausibly resonant systems.The probabilities that offsets in the tidally affected samples this large could occur by chance are about 1 in 18 (equivalent to about 1.6 σ) for the libration-width definition and 1 in 93 (equivalent to about 2.3 σ) for the δres-based definition.Our conclusion that plausibly first-order resonant systems with an innermost planet likely to be affected by tidal dissipation are young does not sensitively depend on the criteria used to identify systems likely affected by tides.

Figure 7 .
Figure7.Distribution of Galactic velocity dispersion based on bootstrap resampling for systems with (blue) or without (orange) a plausibly second-order resonant planet pair according to the δres-based definition.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.The Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems with a plausibly second-order resonant planet pair is significantly colder than the Galactic velocity dispersion of systems that lack such a pair.The probability that an offset this large could occur by chance is less than 1 in 100000, or equivalently about 4.6 σ.This observation implies that systems initially in second-order resonance are driven out of resonance over Gyr timescales.

eFigure 8 .
Figure8.Distributions of Galactic velocity dispersion based on bootstrap resampling for systems with (blue) or without (orange) a plausibly second-order resonant planet pair according to the δres-based definition.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top panels include systems in which the innermost planet is likely affected by tidal evolution, while bottom panels include systems in which the innermost planet is unlikely to be affected by tidal dissipation.Left: Tidal and non-tidal samples separated at τe = 20 Myr.The probabilities that offsets this large could occur by chance is about 1 in 210 (equivalent to about 2.6 σ) for the systems likely affected by tides and about 1 in 76000 (equivalent to about 4.2 σ) for systems unlikely to be unaffected by tides.Right: Tidal and non-tidal samples separated by the median circularization time τe = 400 Myr.The probabilities that offsets this large could occur by chance is about 1 in 280 (equivalent to about 2.7 σ) for the systems likely affected by tides and about 1 in 76000 (equivalent to about 4.2 σ) for systems unlikely to be unaffected by tides.Regardless of the importance of tidal dissipation inside innermost planets, the Galactic velocity dispersion of the sample of Kepler-discovered multiple small-planet systems with a plausibly second-order resonant planet pair is colder than the Galactic velocity dispersion of systems that lack such a pair.The implication is that tidal dissipation is not the process responsible for moving systems initially in second-order mean-motion resonance out of resonance.

Figure 9 .
Figure9.Indirect age indicators for our sample of Kepler multiple-planet system host stars.We plot as blue points systems with at least one plausibly resonant planet pair and as orange points systems that lack a plausibly resonant planet pair.Top left: fromBerger et al. (2018) the equivalent width in m Å of the 670.8 nm lithium feature as a function of GBP − GRP.Top right: fromBerger et al. (2018) lithium abundance as a function of GBP − GRP.At constant GBP − GRP, large equivalent widths or lithium abundances should be observed in the spectra of stars younger than about 650 Myr.We see no significant equivalent width or abundance offsets between plausibly resonant/implausibly resonant systems, suggesting that these plausibly resonant systems are not young in an absolute sense.Bottom left: fromBrewer & Fischer (2018) log R ′ HK as a function of GBP − GRP.Bottom right: fromBrewer & Fischer (2018) SHK as a function of GBP − GRP.log R ′ HK and SHK are measures of the strength of emission in the core of the calcium H & K lines.At constant GBP − GRP, stronger emission should be observed in the spectra of younger stars.We see no significant log R ′ HK or SHK offsets between plausibly resonant/implausibly resonant systems, suggesting that these plausibly resonant systems are not young in an absolute sense.

Figure 10 .
Figure10.Distributions of Galactic velocity dispersions for the indicated populations.We indicate median values as solid vertical lines, 1-σ ranges between vertical dashed lines, and 2-σ ranges between vertical dot-dashed lines.Top: the populations of plausibly first-order (blue) and second-order (orange) resonant multiple-planet systems.Top-middle: the population of implausibly resonant multiple-planet systems.Bottom-middle: the population of all multiple-planet systems.Bottom: the populations of confirmed (orange) and candidate (blue) single-transiting systems.The plausibly second-order resonant system population has a significantly colder Galactic velocity dispersion than the single-transiting system population.The population of single-transiting systems has a Galactic velocity dispersion consistent with the populations of implausibly resonant multipleplanet systems and all multiple-planet systems.
S. Department of Energy Office of Science.The SDSS-III web site is http://www.sdss3.org/.SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofísica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah.The SDSS website is www.sdss.org.SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics -Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.This research made use of Astropy 8 , a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018).This research has made use of NASA's Astrophysics Data System.Facilities: ADS, Exoplanet Archive, Gaia, Kepler, LAMOST, Sloan Software: Astropy (Astropy Collaboration et al.

Table 1 .
Separations from Resonances in Libration Widths

Table 3 .
Resonance-occupation DistributionNote-The number of plausibly resonant pairs in the first-and second-order mean-motion resonances we considered in this analysis.Resonances are ordered by decreasing period ratio Pout/Pin, so planets with initial period ratios larger than 3:1 experiencing convergent disk-driven migration should encounter these resonances from top to bottom.