Period Ratio Sculpting Near Second-Order Mean-Motion Resonances

Second-order mean-motion resonances lead to an interesting phenomenon in the sculpting of the period ratio distribution due to their shape and width in period-ratio/eccentricity space. As the osculating periods librate in resonance, the time-averaged period ratio approaches the exact commensurability. The width of second-order resonances increases with increasing eccentricity, and thus more eccentric systems have a stronger peak at commensurability when averaged over sufficient time. The libration period is short enough that this time-averaging behavior is expected to appear on the timescale of the Kepler mission. Using N-body integrations of simulated planet pairs near the 5:3 and 3:1 mean-motion resonances, we investigate the eccentricity distribution consistent with the planet pairs observed by Kepler. This analysis, an approach independent from previous studies, shows no statistically significant peak at the 3:1 resonance and a small peak at the 5:3 resonance, placing an upper limit on the Rayleigh scale parameter, $\sigma$, of the eccentricity of the observed Kepler planets at $\sigma=0.245$ (3:1) and $\sigma=0.095$ (5:3) at 95% confidence, consistent with previous results from other methods.


INTRODUCTION
The interactions between two planets orbiting a star become especially strong when the planets' periods are integer multiples of one another. These period commensurabilities are called mean-motion resonances (MMRs), and their order is determined by the difference in the integer multiples; e.g. a first-order resonance has an m : (m − 1) period ratio. Different order MMRs have different characteristics, and in this paper we will examine how the symmetrical shape of second-order MMRs affects the aggregate period ratio distribution for a population of planets.
A general overview of MMRs is given in Murray & Dermott (1999) and a recent detailed analytical treatment can be found in Hadden (2019). An analysis of the period ratio distribution of Kepler multi-transiting sys-tems, including features related to MMRs, was in Fabrycky et al. (2014).
When planets are in resonance, their orbital periods vary. The typical rate at which this variation occurs is faster than the orbital precession by a factor of ( M mp ) 1/2 (Murray & Dermott 1999). Therefore, to keep the resonant angles in a librating state, the orbital period ratio averages to very close to the ratio of integers that defines the resonance.
Orbital eccentricities play a role in the maintenance of a resonance, because it is the eccentricity that supplies torque between two planets (Peale 1976). Therefore, the more eccentric the planets of a system are, the more likely they will be in resonance, and the more likely their period ratios will average out to the nominal value of the resonance. The high-order resonances are more sensitive to eccentricity; the strength of the coupling is proportional to the eccentricity to the power of the order of the resonance (Murray & Dermott 1999).
Eccentricities are challenging to measure for transiting exoplanets. The transit timing variation (TTV) method is the most precise for measuring eccentricities, but it favors planets near first-order MMRs and presents a degeneracy between mass and eccentricity (Lithwick et al. 2012).
Examining the effect of eccentricity on transit durations is a more general but less precise technique. Fabrycky et al. (2014) (building on Moorhead et al. (2011)) compared transit durations between planets in the same system to derive a typical Rayleigh scale. Ford et al. (2008) proposed a similar technique with dependence on the stellar density; later studies were able to apply this technique as improved stellar parameters became available (Van Eylen & Albrecht 2015;Xie et al. 2016;Van Eylen et al. 2019;Mills et al. 2019). Shabram et al. (2016) used transit and occultation detections to measure projected eccentricity via transit duration ratios and phase offsets, independent of the stellar density, but this probes generally the Hot Jupiter population.
Alternatively, forward modeling can use observed populations to constrain intrinsic population distributions, e.g. He et al. (2019He et al. ( , 2020. The general picture that has emerged from these studies are that multi-planet systems of small (i.e. nongiant) exoplanets tend to have low eccentricities. Both multi-planet giant systems, such as those probed by radial velocity measurements, and single-planet systems tend to have a distinctly excited eccentricity distribution compared to multi-planet sub-giant systems.
We develop a new method to measure the eccentricity distribution of multi-planet systems observed by Kepler using second-order MMRs. This method has different associated systematics and does not depend on precisely measured stellar properties or specific transit observations. Although this method relies on MMRs, given that migration likely doesn't easily trap planets in second-order MMRs (Delisle et al. 2014;Migaszewski 2015;Xu & Lai 2017), we expect it traces post-scattering systems.
The structure of this paper is as follows: Section 2 discusses the methodology used to generate and analyze our simulated systems; the results are covered in Section 3; an analytical approximation is discussed in Section 4; Section 5 discusses the importance of the results and further considerations, including consideration of the TTVs of the observed planets in K02261; and Section 6 summarizes our key findings.

METHODS
For our investigation, our goal was to create a synthetic population of planet pairs that was similar to those found in the Kepler sample. We intentionally em-ulate the properties of the sample as observed; i.e., we seek to probe the observed distribution and not the intrinsic distribution. To isolate the effects of eccentricity, we use a wide variety of planet pairs like those observed by Kepler to marginalize over other factors, particularly the mass which has a strong degeneracy with eccentricity in MMRs (see Section 5.4). We then seek to quantitatively compare the resulting second-order MMR peaks at 3:1 and 5:3 for various eccentricity distributions, particularly in comparison to that of the Kepler sample.

System Generation
In order to numerically investigate planetary behavior near the second-order MMRs, we needed a population of planet pairs. Our intent was to create a population that was similar to that of Kepler, although without being rigorously Kepler-like.
To develop a population of planet pairs near the 5:3 and 3:1 MMRs, we first randomly generated a large number of planet pairs evenly in period ratio near the commensurability location. The mutual inclination of the two planets was chosen randomly from a Rayleigh distribution with scale 0.032 radians (Fabrycky et al. 2014). The inner planet period, each planet mass, and the stellar mass were randomly assigned using distributions based on the population of multi-planet systems observed by Kepler. Arguments of pericenter were chosen uniformly between −π and π.
This set of generated planet pairs is intended to be used to probe pairs at various eccentricities. Therefore, to ensure a wide coverage of eccentricity, several methods of assigning eccentricities were used, including uniformly between 0 and 1, uniformly between 0 and 0.5, and linearly decreasing from 0 to 0.5. For each planet pair, multiple systems were created varying only in the mean longitude of the planets.
This initial population of planet pairs was then integrated for 10 6 orbits of the outer planet to remove any "unstable" pairs. In this context, pairs were considered unstable when they departed from the period ratio range of interest and not as a strict determination of instability. The bounds at which pairs were removed were <2.5 or >3.5 for the 3:1 resonance and <1.55 or >1.88 for the 5:3 resonance. The integrations were computed using the IAS15 integrator in the REBOUND package (Rein & Liu 2012;Rein & Spiegel 2015).
At the end of the first stability run, 50 surviving systems were randomly chosen for each resonance to be integrated for ten times longer. Of these 50 systems, an additional 4 went "unstable" (departed the specified period ratio bounds) for the 3:1 systems and an additional 8 for the 5:3 systems, with no apparent preferential lo-cation in period ratio. These low percentages indicated that the initial stability run was sufficient to ensure a robust sample that would not be dominated by planet pairs on the verge of instability.
From the planet pairs that remained in the sample after the stability cut, one pair was randomly chosen from each initial set of mean-longitude systems to be used in the analysis, and thus each pair used in the analysis has unique parameters. Finally, of these pairs, those that had an average period ratio in close proximity to exact commensurability were selected for the analysis. This last cut was done between 2.94 and 3.06 for the 3:1 resonance (3.00 ± 0.06) and between 1.62 and 1.713 for the 5:3 resonance ( 5 3 ± 7 150 ). The limits are chosen to be wide enough to capture the expected resonant behavior but narrow enough not to include additional higher-order resonances. Specifically, the limits around the 5:3 were chosen to exclude the 8:5 and 7:4 meanmotion resonances.
Detailed distributions of the properties of the final generated sets of 913 (near 3:1) and 670 (near 5:3) planet pairs are shown in Appendix B.

Analysis
The final results of simulations described in Section 2.1 provide an instantaneous snapshot of the orbital properties of the planet pairs. These properties change over time, however, as the planets perturb each other, particularly planets that are caught in the mean-motion resonance and librate around a fixed point configuration. Given the time baseline of the Kepler mission, we can actually observe the average orbital properties. Thus the first step of the analysis was to integrate each planet pair for 3.5 years and to time-average the osculating orbital properties over this period. This is done by saving the orbital properties at 2000 steps during the 3.5 year integration and then taking the mean.
The approximate TTV period of planets near a m : n mean-motion resonance is given by Deck & Agol (2016): and pairs actually in the resonance will experience libration around the center of the resonance on a timescale shorter than that (Nesvorný & Vokrouhlický 2016). Calculating the TTV period for the final planet pairs in our set of generated systems and comparing to the time-averaging period of 3.5 years, we find that the majority of these planet pairs, 54.2% for the 3:1 and 81.8% for the 5:3, will have completed at least one libration cycle. We can also calculate an analytical approximation of the libration periods with the Andoyer module of the celmech package. This predicts a libration period of less than 3.5 years for 31.2% for the 3:1 pairs and 55.7% for the 5:3 pairs. Note however that this approximation generally assumes low eccentricities, which is not the case for many of these pairs. Thus we expect that for both generated populations, the average period ratio would fall much closer to the nominal resonance period ratio than the instantaneous period ratio.
This peaking can be seen in Figure 1. There is also some stability sculpting from the initial uniform period ratio distribution (see Figure 21). This sculpting is expected, as being near resonance can increase instability, but being in resonance can provide protection against instability, thus leading to an excess of planet pairs near each resonance. However, at any random point, the peak of the excess planet pairs is much more spread out and indistinct than the time-averaged peak. This narrow time-averaged peak is a distinct signature, and is even expected when there is no excess planet pairs (see Figure 9). In our generated systems, we find that no more than approximately one-third of the planet pairs located in the time-averaged peak can be attributed to a stability-induced excess. Thus the effect of sculpting from instability is both expected in real populations and not necessary for the appearance of the time-averaged resonant peak.

Model
To quantify the size and shape of this peak, we use a model comprised of a linear background and a Gaussian peak centered on the nominal resonance value. The normalized equation for this model is and it has three parameters to be fitted: c 1 , c 2 , and c 3 , where c 1 is the slope of the linear background distribution, c 2 is the amplitude of the Gaussian peak, and c 3 is the width of the Gaussian peak (standard deviation). The N indicates a normal distribution for the given center and width. The period ratio is denoted by r, and r 0 is the nominal period ratio of the mean-motion resonance. The limits on the period ratio range, r min and r max , are taken to be 2.94 and 3.06 for the 3:1 resonance (3.00 ± 0.06) and 1.62 and 1.713 for the 5:3 resonance ( 5 3 ± 7 150 ).

Fitting
To fit for these parameters to our simulated period ratio distribution, we implemented a Markov chain Monte Carlo (MCMC) sampler using the emcee package  Figure 1. The probability density of period ratios near the (a) 3:1 and (b) 5:3 mean-motion resonances. The nominal resonance value is marked with a dotted line. The instantaneous final and averaged final distributions are shown for the specified set of generated systems as described in Section 2.1. There is an excess of planet pairs due to the resonance seen in the instantaneous distribution (a result of stability sculpting), but the signal of interest is the narrow peak that is clearly revealed by time-averaging. (Foreman-Mackey et al. 2013). We treat each period ratio in our sample as a Poisson process, calculating their probabilities as individual independent events, which requires no binning of the data. Given that we expect to find a well-defined peak for this set of systems, the priors are uninformative, using a flat log prior for the peak height and width and a prior for the slope that is flat in the angle of the line with the axis, to avoid biasing the slope fit towards higher slopes.
We initialized the MCMC analysis with 32 walkers based on an initial least squares fit to the data. The chains were considered converged when the number of steps was more than 100 times the maximum autocorrelation time calculated using emcee and the autocorrelation time changed less than one percent over 100 steps (<10,000 steps in total). We used a burn-in of twice the maximum autocorrelation time and thinned the sample by half the minimum autocorrelation time.
After fitting for the full set of systems, we then use the resulting posterior of the peak width as a prior for the remaining fittings, including fits to the observations. As some of these fits do not have well-defined peaks, keeping this strict prior allows us to compare whether similar peaks are present in other sets of systems. Without this prior, for a set of systems with no peak, the peak height and width are highly correlated and the peak width is unconstrained. Additionally, this narrow width is a signature of the peak that we are interested in, and using this strict prior ensures we are investigating the peak arising from the time-averaging effect and not, for example, from stability-enhanced MMR protection.
We fit for the Kepler data set around each commensurability to find values for the three parameters. The KOIs included in this data set are listed in Table 4. We take the KOIs from the NASA Exoplanet Archive 1 for all candidate and confirmed KOIs. Including candidate KOIs may mean some pairs are not genuine planets. Looking at planet pairs is an intrinsically multi-planet phenomena, which leads to some detection biases from masking of transits in the light curve (Zink et al. 2019). For planets with periods <200 days, this is only a 5.5% efficiency loss, so we neglect this potential bias. Additionally, eccentricity affects the probability of transit (Burke 2008). However, these effects are offsetting, and any net effect would be to make it more likely for eccentric planets to be detected.
To quantify how many additional planets are present due to the time-averaged resonance peak, we find the 1 https://exoplanetarchive.ipac.caltech.edu fractional area of the Gaussian curve that is in excess from the background. Figure 2 shows the eccentricity distributions of the generated system sets. While there is a clear amount of stability sculpting in the eccentricity distributions, there are many highly-eccentric planet pairs that remain in our stable samples. In order to compare the Kepler distributions to various eccentricity distributions, we take sub-samples of our generated system sets that are consistent with a given eccentricity distribution. This sampling is done by calculating a probability for each eccentricity using the desired eccentricity probability distribution function and then using those probabilities as weights for drawing planet pairs without replacement. The sub-sample size was fixed at the same number of planet pairs as in the Kepler samples, 30 for the 3:1 resonance and 73 for the 5:3 resonance to allow for the same amount of systematic error from the sample size. The sub-samples were checked against a sample of the same size generated from the desired distribution using a 2-sample Kolmogorov-Smirnov test and a 2sample Anderson-Darling test and accepted only if both p-values were ≥ 0.1.

Sampling and Fitting Eccentricity Scales
The eccentricity distributions used are visualized in Figure 3. Given that we are considering planet pairs, the eccentricity distributions are those for combined eccentricities, e 1 +e 2 . We considered the individual eccentricities to be independent.
The combined Rayleigh distributions were calculated analytically. The distribution of the sum of two independent variables is expressed by the convolution of the two distributions, P (x 1 + x 2 ) = P (x 1 ) * P (x 2 ). Thus, the sum of two independent Rayleigh random variables with the same scale parameter σ has a probability density function of the form ) .
(3) The combined beta distribution was calculated numerically, incorporating the uncertainties in the distribution parameters from Kipping (2013). A set of 10 6 eccentricities were drawn, each time with parameters drawn from a normal distribution with the mean and sigma given by the error bars from Kipping (2013), α = 0.867 ± 0.044 and β = 3.03 ± 0.17 (this is slightly modified from β = 3.03 +0.17 −0.16 as listed by Kipping (2013) for ease of computation, with negligible effects). The probability density function was then calculated using a Gaussian kernel-density estimate. and end (solid green) of the stability run. Their similarity shows they are in statistical steady-state. The gray background shows the initial combined eccentricity for all planet pairs, including those that went unstable and were excluded from the final selection. There is some stability sculpting, particularly at higher eccentricities, but many highly eccentric planet pairs remain. probability density combined beta, K13 combined Rayleigh, HL14 combined Rayleigh, 3:1 this work, upper limit combined Rayleigh, 5:3 this work, upper limit Figure 3. Several probability density functions with example parameters for the combined eccentricity of the planet pairs (e1 + e2). The combined beta distribution was derived numerically from samples taken in accordance with Kipping (2013) parameters and plotted using a Gaussian kernel-density estimate. Combined Rayleigh distributions for various scales were calculated analytically from Equation 3 and plotted for several samples scales: the best-fit value from Hadden & Lithwick (2014), σ = 0.018; the upper limit for 3:1 from this work, σ = 0.245; and the upper limit for 5:3 from this work, σ = 0.095.
We repeated this ten times, using a new random subsample, for each eccentricity distribution and combined the resulting chains to get the final posterior for a given eccentricity distribution. For the beta distribution, we increased the number of calculations to fifty, but found negligible improvement in the fit and so only ten were used for the remaining sub-samples. Finally, we used the MCMC analysis on the sub-samples of the generated systems and calculated the fitted fractional peak area at each resonance. The peak area was found by calculating the area due solely to the continuum distribution and subtracting that from one (the normalized area). The results of this fitting and comparison with the Kepler data are discussed next.

RESULTS
Our first result is the detection and quantification of the presence of an excess in planet pairs near the 3:1 and 5:3 resonance in our simulated systems. The fitted parameters are shown in Table 1, with errors listed from the 68.2% interval, along with the derived parameter of peak area. The data and median models, along with a sample of posterior draws, are shown in Figure 4.
Next, we classify the presence of a peak near the 3:1 and 5:3 resonances in the Kepler data. The 3:1 for Kepler has no indication of a peak, and in fact the median model is a slight dip at the resonance location; however, the model fit is consistent with a flat background and zero peak. The 5:3 for Kepler has a small peak at the resonance location. The fitted parameters are shown in Table 1, with errors listed from the 68.2% interval, along with the derived parameter of peak area. The data and median models are shown in Figure 5. While the models are not as closely fit to the data as with the generated systems, this can be attributed to the much smaller sample size. For illustration, we drew ten samples with size equal to that of the Kepler data from the median model for each resonance, also shown in Figure 5. For both resonances, the Kepler distributions fall within the range of these samples.
We use our sub-sampling as described in Section 2 to generate a posterior for the peak area of a set of planet pairs with eccentricity consistent with a beta distribution as described in Kipping (2013). These posteriors, along with the posteriors of the peak area for the Kepler samples, are shown in Figure 6. The Kepler samples are not consistent with the beta distribution to a level of 1.8σ (92.0% confidence) for the 3:1 and 3.2σ (99.9% confidence) for the 5:3. The fitted parameters for the beta distribution sub-samples are shown in Table 1, with errors listed from the 68.2% interval, along with the derived parameter of peak area.
We do a similar comparison for sub-samples of sets of planet pairs with eccentricity consistent with a Rayleigh distribution of various scales, from σ = 0.005 to 1.000 in steps of 0.005. Some scales were poorly sampled by our generated systems and were not analyzed or included in the results. Figure 7 shows the results for two of the eccentricity scales for each resonance, including the cumulative distribution for each of the ten random subsamples as well as the model resulting from the median of the combined MCMC posteriors. At low eccentricities, the distribution has very little peak, as expected. This shows up as almost entirely continuum in our model, which describes the distribution well. At the higher eccentricity, the peak starts to become clear. Our model captures this changing behavior of the distribution, allowing us to quantify the fractional peak area for comparison with the model fitted from the Kepler data.
The combined results for all the fitted scales are summarized in Figure 8, using a representation of the posteriors of the peak area (error bars are the 68.2% interval) at various scales and the Kepler data for comparison.
We find that, at 95% confidence, the upper limit for the Rayleigh scale consistent with the Kepler data is (75 ± 6) × 10 −5 0.48 ± 0.08 Table 1. Summary of the fitted parameters from our MCMC posteriors. c1 is the slope of the linear background distribution, c2 is the amplitude of the Gaussian peak, and c3 is the width of the Gaussian peak (standard deviation). The prior on c3 used on the all generated systems set (*) was different from the other fittings; see text for details. The fractional peak area, A pk , is derived from the fitted parameters to quantify the additional excess or lack at the resonance location. We include here only the Rayleigh distribution results corresponding with 1-, 2-, 3-, 4-, and 5-σ discrepancy with Kepler.
0.245 for the 3:1 and 0.095 for the 5:3 resonance. The fitted parameters for a select number of the Rayleigh distribution sub-samples are shown in Table 1, with errors listed from the 68.2% interval, along with the derived parameter of peak area. These selected scales are 1-, 2-, 3-, 4-, and 5-σ discrepant with Kepler, and are also shown in Figure 8 as dotted vertical lines.

ANALYTICS
To better understand the physical underpinnings of our results from numerical simulations, we investigate an analytical approximation for the effects of a secondorder mean-motion resonance on the observed period ratio distribution. We use the framework of the integrable model developed by Hadden (2019), hereafter H19. This model is applicable for MMRs of any order interior to 2:1, and so we focus our application here on the 5:3 resonance. While this analytical model is not suitable for the 3:1 resonance, we expect the actual dynamics to be similar, as seen as well in our simulations.
There are two relevant pieces of information that we obtain from the model: the resonance width (in period ratio space) and the libration timescale. We use the integrable one-degree-of-freedom Hamiltonian H(J, θ; J * ) as given by Equation 19 of H19; see that paper for full details on the derivation and application.
In brief, we first calculate the conserved quantity J * for a given system (from H19 Equations 9, 20, and the relation J * = J − kP ) and then find the unstable fixed point at θ = 0 by finding the value of the conjugate momentum J that maximizes H and therefore determines the energy of the separatrix (E sep ). The maximum width is calculated by solving H(J, π/2; J * ) = E sep for J. The two J values can be converted back to orbital elements using provided relations (H19 Equation 21), giving a minimum and maximum period ratio for the resonance. To find the period of libration, we use an approximation that linearizes the equations of motion in the vicinity of the equilibrium, using the Andoyer module of celmech. This approximation is only accurate for small-amplitude librations, and a more detailed analytical model would develop this further.
To investigate the period ratio distribution, we use a simple model with the resonance width and the period of libration. First, we calculate the resonance width for a given system configuration of masses, periods, and eccentricities.
If the period ratio of a planet pair is within the resonance width, the libration over time is represented by a sinusoid with the calculated libration period, initial phase as given by kθ 0 (see H19 for details), and centered at 5 3 . The amplitude of the sinusoid is drawn randomly    from a linear probability distribution between minimum at zero amplitude and maximum at an amplitude of half the resonance width. For period ratios outside the resonance width, the period ratio is assumed to be constant. This is a simplification, as the period ratios would likely change over time, but they would not follow the resonant model and so these changes are negligible to our result.
Finally, we take a time period to average over and "observe" the period ratio based on our libration model over that period n times where n = int tavg P2 .
We use this analytical model in two complementary ways. First, we apply the model to a Kepler-like population of planet pairs for comparison with our N-body results. Second, we apply the model to a single fiducial planet pair and investigate the effect of changing individual parameters of the system.
The systems for the analytical population model are generated according to the same process as described in Section 2.1. These pairs are uniform in their initial period ratio between 5 3 ± 7 150 , and each planet's eccentricity is randomly assigned using a Rayleigh distribution. We use three scales, a low-eccentricity σ = 0.05, mediumeccentricity σ = 0.15, and high-eccentricity σ = 0.50. Figure 9 shows the resulting probability distributions and cumulative distribution functions for the three populations. As expected, a very narrow peak forms at the exact resonance location, the height of which increases with increased eccentricity. Because the initial period ratio is exactly uniform, the addition of planet pairs to this peak must result in a lower density of planet pairs outside the peak. These troughs are seen in the distributions. However, because the resonant and eccentric pairs end up at a very narrow location but can start from a wide range of initial period ratios that varies based on the properties of each system, the troughs are very spread out compared to the peak.
To examine the effect of various parameters on the period ratio distribution, we establish a fiducial system and then vary parameters one at a time. The fiducial system properties are given in Table 2. The properties were chosen to be consistent with Kepler-like systems but without any special consideration. Then we investigate the period ratio libration for a grid of 500 period ratios (by changing the period of the outer planet) between our limits of interest, 5 3 ± 7 150 . The resulting cumulative distribution function in period ratio is shown in Figure 10, both for the average period ratio and a period ratio chosen randomly within the time period. As expected, a very clear and narrow peak forms when the period ratios are averaged over time. We also ran an N-body simulation with the same parameters for 3.5 0.5 -10 Table 2. Variable parameters used for analytical investigation. When varied, each parameter is varied between the minimum and maximum values in 20 linear steps. P2 is determined for each system based on the period ratio.
comparison. While the distributions are not exact, the agreement is good for this simplified model. It is also interesting to note in the N-body distribution what seem to be similar (though much weaker) signals at the 13:8, 18:11, and 17:10 resonances. In this analytical model, we do not include a continuum such as is present in our model in Equation 2. Instead, this analytical model examines only how the peak super-imposed on the continuum is affected by the changing system parameters. This, and the fact that all of the systems in this analytical model have the same properties aside from the period ratio, is the difference in the shape of the cumulative distribution functions in Figures 4 and 10. The height and width of the peak is our interest, and thus what we investigate via this analytical model.
The results of varying the nine parameters are shown in Figure 11. Only the averaged period ratios are shown here for clarity, with the color indicating the parameter's value as indicated by the color bar to the right of each plot. For each parameter, the fiducial value is plotted as a white dot on the color bar for reference.
First we see that varying the masses of the planets (11a and 11b) has a noticeable but small effect. The effect does not depend on which planet's mass changes, as the resonance is sensitive to the planets' combined mass. The amplitude of the peak does not change much, but the width does change slightly and the center of the peak shifts. This is due to a change in the libration period-if the libration period increases past t avg , the averaged value will no longer lie at the center as a cycle has not completed. Changing the mass of the star (11c) has a generally similar effect, though there is a stronger effect when reaching much larger mass ratios (i.e., a smaller stellar mass). As expected, we see a very strong dependence of the peak on eccentricity (11d and 11e). As eccentricity is increased, the amplitude of the peak at resonance becomes much larger. In this case, the effect is notably stronger when e 2 is increased. The orientation of the eccentricities (11g and 11h) also has an effect, although less than from the magnitude. In our fiducial system, the eccentricities are 0.5 rad (∼ 30 • ) out of alignment. As either of the changes towards ±π we see an increase in the amplitude of the peak, suggesting a stronger signal for anti-aligned eccentric planet pairs.
Changing P 1 (11f) has no apparent effect on the amplitude or width of the peak, as is expected. However, the location of the center of the peak changes significantly. This is because the libration period depends on the planets' periods, and so for longer period planet pairs, the libration cycle has not yet completed during this t avg . We can similarly see this effect from changing t avg (11i). As t avg decreases, fewer of the systems have completed their libration and so the center shifts away from the resonant location. In this case, the center of the peak shifts leftward for incomplete libration cycles. This depends on the initial phase of the libration. For different initial phases, the peak can also shift rightward, and so the direction of the peak shift should not be seen as significant.
We note a few remarks about these analytical results. First, we use a very simplified model for the libration. Several N-body comparisons were done to ensure this model was reasonable to use as an approximation, but it does not give exact behavior. Second, we use the simplified Hamiltonian from H19 that is based on a leadingorder approximation in eccentricity. This introduces additional inaccuracies, especially for our high-eccentricity systems. Third, there are no stability cuts or considerations, and it is probable that not all of these systems are viable. Fourth, this analytical model assumes a planar configuration and no inclination effects.
Given these points, we do not provide this analytical investigation for a quantitative comparison. However, we expect the qualitative trends to be accurate and find the overall behavior to be in line with our expectations, confirming the resonant behavior as seen in our numerically simulated results.

DISCUSSION
It is particularly striking how narrow the resonance peak is. This arises from the time-averaging of the orbital parameters and is a clear feature of the secondorder resonance. The narrowness of the peak is, for example, what allows us to ignore the effects of any associated troughs in our analysis.
Our sets of generated systems are intended to create sets of planet pairs with a large array of eccentricities limited only by stability. We are agnostic to formation mechanisms, i.e., whether such systems are physically realistic. This allows us to marginalize over other properties to isolate the effects of eccentricity on the period ratio distribution. Other potential effects are discussed below.

Comparison with Previous Results
The first results pertaining to the eccentricity of exoplanet populations were obtained via radial velocity measurements. These systems were typically composed of giant planets and generally excited in eccentricity. Wright et al. (2009) found a mean eccentricity of 0.22 for multi-planet systems, their data containing many more giant systems than in the Kepler observations. Similarly, Kipping (2013) quantifies the eccentricity distribution of radial velocity planets with a beta distribution that gives a mean eccentricity of ∼0.22.
As more transiting planet data became available via Kepler, and different techniques were developed for constraining the eccentricity of these planets, a different picture began to emerge. These results indicated much lower eccentricities are typical for small planets. Kane et al. (2012) found that for giant planets the Kepler observations were in agreement with previous radial velocity results, but that planets smaller than Neptune became "rapidly and significantly more circular." Our constraints on the planets near the 3:1 and 5:3 are consistent with the results obtained previously for small, multi-planet systems observed by Kepler. Our 5:3 constraint, which would predict a mean eccentricity of no more than 0.12 at 95% confidence (for a Rayleigh distribution, E = σ π 2 ), is not consistent with the eccentricities of giant planets such as those analyzed by Wright et al. (2009). We find a 1.8σ (3:1) and 3.2σ  centricity decreasing as n increases. This multiplicity dependence has also been found by Bach-Møller & Jørgensen (2021). Additional evidence for a difference between single-planet and multi-planet eccentricities was found in Xie et al. (2016), Mills et al. (2019), andVan Eylen et al. (2019). This dichotomy may be due to unobserved outer companions (Poon & Nelson 2020). As our method depends on measuring a period ratio, which must have two planets, we probe the observed n ≥ 2 Kepler sample. He et al. (2019) and He et al. (2020) use forward modeling to constrain properties of the intrinsic population of Kepler planets. While in this work we are constraining the eccentricities of the observed population, rather than the intrinsic population, we find the comparison of interest. He et al. (2019) find low eccentricities, with characteristic Rayleigh scales σ 0.01 − 0.02. He et al. (2020) uses an angular momentum deficit approach for the modeling and finds a median eccentricity of 0.138 for two-planet systems, decreasing down 0.029 for fiveplanet systems and to 0.009 for ten-planet systems (see their Table 5). These results are consistent with ours, despite the differing populations being described.
Interestingly, He et al. (2020) find that a log-normal distribution is a better fit for the eccentricity distribution, rather than a Rayleigh distribution which is more typically assumed (e.g. He et al. 2019, Fabrycky et al. 2014). Tremaine (2015) and Shabram et al. (2016) also recommended a non-Rayleigh distribution, while Van Eylen et al. (2019) consider several models and do not advocate one over another. We leave the investigation of additional eccentricity distributions to future work.
Our results, then, are consistent with other studies and further evidence that the eccentricity distribution of small (i.e. not-giant) planets in multi-planet systems have a different underlying eccentricity distribution than single-planet or giant-planet systems.

Tides
Tides might act to remove some of the planet pairs from resonance or decrease their eccentricity over time (Novak et al. 2003;Delisle et al. 2014). Here we consider the impact of tides on these relatively close-in planets.
We calculate the circularization timescale for our planet pairs using Equation 9 of Rasio et al. (1996). The resulting timescales, scaled by Q, are shown in Figure 12. For terrestrial planets, typical Q values range from 10 to 500 (Goldreich & Soter 1966). We would thus expect approximately half of our sample would have undergone one tidal circularization timescale in 100 Myr.
However, this does not affect our results, which are focused solely on the effect of eccentricity in shaping the period ratio distribution and not in the feasibility of such an eccentricity distribution. This tidal timescale is a consideration for other limits on the stable eccentricity of a similar population.
Furthermore, we can split our sample sets into two equal parts based on the scaled circularization timescale to see if systems with faster circularization timescales have different properties in this context that might be affecting our results. Instead, we find that the resulting averaged period ratio and eccentricity distributions are not statistically different. The p-values from a 2sample Kolmogorov-Smirnov test for the split samples of the 3:1 set are 0.70 for the averaged period ratio distribution and 0.17 for the averaged combined eccentricity distribution. For the 5:3 set, the corresponding p-values are 0.13 and 0.34, respectively. Thus, we expect tides to be an interesting but negligible consideration for our results.

Eccentricity at Formation
In our approach to creating our sets of planet pairs, we did not use a physical motivation for the eccentricities of the planets. Indeed, we purposely attempted to cover a large eccentricity space to allow for sub-sampling by eccentricity. We did, however, use a stability analysis to cut unstable planet pairs from our sample (see Section 2.1). This stability cut produced some shaping of the eccentricity distribution (i.e. Figure 2), in particular preferentially removing high-eccentricity pairs.
The stability-shaped eccentricity distribution is still much too eccentric to match the Kepler data, however. Because stability alone cannot explain the low observed eccentricities, this implies some other process must prevent high-eccentricity pairs, likely related to formation mechanisms.
The eccentricity distribution predicted by formation models varies depending on many parameters, particularly related to the properties of the protoplanetary disk (Fortier et al. 2013). Effects that serve to excite eccentricity include reduced solid surface density (Van Eylen et al. 2019), shallower solid surface density profiles (Moriarty & Ballard 2016), the presence of additional companions (Huang et al. 2017;Anderson & Lai 2017), and others. Conversely, there are effects that can dampen eccentricity, such as the presence of planetesimals to provide dynamical friction during formation or gas in the disk (Mulders et al. 2019). Recent findings by Yee et al. (2021) indicate that the observed eccentricities (of those measured by Hadden & Lithwick (2017)) are lower than predicted stability limits or from giant impact formation theory. Fitting the exact formation   Figure 10. The analytically-predicted cumulative distribution function in period ratio for the fiducial system. The time-averaged period ratio is shown in solid magenta and a random period ratio from within the time period is shown in dashed magenta. The same distributions are shown for a comparable N-body simulation in cyan (solid is averaged and dashed is random). The black dashed line shows the initial flat period ratio distribution. The resonance location at 5 3 is indicated with a vertical dotted line. There is reasonable agreement between the numerical and analytical distributions, and in both cases, the peak at resonance only appears when the period ratios are averaged.
conditions consistent with the Kepler planets is the beyond the scope of this paper.

Mass Dependence
Mass and eccentricity both have strong effects on the behavior of planets in mean-motion resonance. This can be seen, for example, in Figure 11. By using a distribution of masses and mass ratios in our set of simulated planet pairs that approximates that seen in the Kepler data, we intend to marginalize over the mass effect to isolate the effect from eccentricity. To illustrate how the mass and eccentricity effects interact, we split our set of planet pairs into 2 groups in mass ratio and 2 groups in eccentricity. The split was chosen at the median values in each case; e 1 + e 2 = 0.58 and log m1+m2 M = −4.42 for the 3:1, e 1 + e 2 = 0.28 and log m1+m2 M = −4.44 for the 5:3. Figures 13 and 14 show the distributions for these split sets of planet pairs. The eccentricity distribution is largely unaffected by the mass ratio cut (13a/14a), and the same is true for the mass ratio distribution and the eccentricity cut (13c/14c). However, strong differences are seen in the distributions of the period ratios.
When looking across all eccentricities, the higher mass ratio half of the simulated systems has a more notable peak at the exact commensurability location, but the peak is still clear for the lower mass ratio planet pairs (13b/14b). Conversely, when looking across all mass ratios, the higher eccentricity half of the data has almost the entirety of the effect at the commensurate period ratio, and the peak almost vanishes for the lower eccentricity half of the data (13d/14d).
Thus, while mass and eccentricity both have strong effects on the resonant behavior, the eccentricity effect is much stronger. Given this difference in strength and the marginalization over the apparent mass ratio distribution of Kepler planets, our results robustly isolate the eccentricity effect near the 3:1 and 5:3 resonances.

Resonance Probability
Using our fitted model for the Kepler period ratio distribution, we can calculate a resonance probability for a given period ratio by comparing the height of the full model at that period ratio to the height the model predicts for the continuum at that period ratio.
When applied to the Kepler 3:1 sample, the resonance probability for all the planet pairs is zero. This is due to the fact that our fitted model predicts no peak at the resonance.
For the Kepler 5:3 sample, there are 10 KOI pairs with non-zero resonance probability. These are listed in Table 3 with their calculated resonant probability. The resonance probabilities are plotted in Figure 15. In Section 5.6 we investigate one of these KOI pairs, K02261, in more detail.

Investigating the Case of K02261 (Kepler-1164) via Timing Analysis
There are two KOIs in this system: K02261.01 and K02261.02. The former is confirmed (Kepler-1164 b, Morton et al. 2016); the latter is a candidate. The K02261 system is also known as Kepler-1164.
To investigate this system, we use the transit numbers, times, and variations from Rowe et al. (2015). There are 315 transits of K02261.01 and 157 transits of K02261.02 covering a time period of approximately 1455 days. The pair's average period ratio over this time is 1.666(4).
Transit data cannot give an instantaneous value of the osculating period ratio as measured exactly from an Nbody simulation. Instead, we approximate the period ratio at any given point in time using the last two observed transits of each planet prior to that point. The average period for each planet is approximated as the difference in transit times divided by the difference in transit numbers, and these periods are used to calculate an approximate period ratio for that point in time.

KOIs
Kepler  Table 3. KOIs near the 5:3 with non-zero resonance probability. The associated Kepler name is listed where applicable.
The variation of period ratio over time for K02261 according to this method is shown in Figure 16. Errors are incorporated numerically by repeating the calculation 1000 times, drawing transit times randomly from a Gaussian described by the transit time and TTV uncertainty for each transit from Rowe et al. (2015). We use 600 evenly spaced times from the observing period to calculate the period ratio.
Due to the scatter and uncertainty of the observational data, it is possible but not evident that the period ratio of the planet pair may be oscillating around the exact commensurability over time. This is the case even when binning the data or attempting to phase fold on various periods. We use N-body models of the planet pair to investigate further.
First, we need to have a defined set of system parameters in order to produce an N-body model. Given the uncertainty of the data, we do not aim to constrain a fit but merely find sample systems that are consistent with the observational TTV data.
From the least squares optimization function in SciPy (Virtanen et al. 2020), we obtain a fit and covariance matrix using the observed transit times. We use these to generate 10 5 samples from a multivariate normal distribution, enforcing a requirement that the samples must be physical: that is, not allowing negative masses or eccentricities greater than 1. We use the analytical stability criterion from Hadden & Lithwick (2018) to ensure these samples are stable, with Z < Z crit , where Z is the complex relative eccentricity (see their Equations 15 and 19).
From these stable, physical samples that are consistent with the observational data, we can select examples to illustrate the possible behavior of K02261. We select a subset of 10 random samples that show libration in the resonant arguments of the 5:3 resonance and a sub-set of 10 random samples that do not. From these two example models, we repeat the period ratio over time calculation that was done for the observed transit times in Rowe et al. (2015). This model uses the same transit numbers as Rowe et al. (2015), and so any observational gaps in the data will be re-created and associated systematic uncertainty in this method of approximating the period ratio from transits is maintained. There is no uncertainty associated with the transit times themselves, however, being calculated exactly from the N-body simulation. For comparison, we also directly take the osculating period ratios from the N-body simulation at the same points in time.
From these twenty sample systems, we can compare the osculating period ratios estimated from the transit times and taken exactly from the simulation. We compare these two period ratio measurements in Figure 17 and find that the estimation from the transit times is highly correlated with the exact period ratios over time, with a Pearson correlation coefficient of 0.834 for the entire set of twenty. Looking at the same relation in the sample systems categorized as resonant, the Pearson correlation coefficient is 0.925, and for those categorized as non-resonant, the Pearson correlation coefficient is 0.777. These high correlations indicate that this estimation approach is a useful measurement of the period ratio changing over time, without requiring an assumption about the resonant state.
Next, we use these samples to examine how the period ratio changes over time. During the Kepler observation period, there is little structure to the period ratio change. This is likely due to the small expected masses of the planets, both of which are consistent with zero from our stable, physical samples resulting from the least squares fit. However, if we extend the time period another ten Kepler observation periods, structure be-  Figure 11. The analytically-predicted, time-averaged cumulative distribution functions in period ratio for various systems. The parameter being varied is indicated by the color bar to the right of each plot, where the fiducial value is plotted in white for reference. The resonance location at 5 3 is indicated with a vertical dotted line. All the parameters affect the height and width of the peak, the most notable being the eccentricity of either planet, shown in panels (d) and (e). The circularization timescale of the inner planet of our generated populations, scaled by Q. Based on a typical Q value of ∼100 for terrestrial planets, we expect that perhaps half of our sample may have undergone one tidal circularization timescale in 100 Myr, indicating that maintaining these planets at such high eccentricities may be difficult.
gins to emerge, as shown in Figure 18. The period ratio varies periodically, although not always sinusoidally, and while there are many types of variation within these 20 samples, some difference between the resonant and nonresonant samples is visible. In particular, the mean of the period ratio over the entire time is shown with a horizontal line for each sample, and the resonant samples clearly cluster closer to the 5:3 commensurability while the non-resonant samples tend to vary around different period ratios. However, when we compare these models to the observational data, it is clear the errors from measuring the transit times dominate over any of these effects-contrast the range in the y-axis in Figures 16 and  18. From this, we infer that the observational data from Kepler for specific systems near the 5:3 resonance do not contradict our conclusions in this paper, but that more precise transit times observations are needed to place any substantive constraints on these individual systems.

CONCLUSION
In this work, we demonstrated that highly-eccentric planets produce a very narrow peak near the secondorder resonance locations of the 5:3 and 3:1 MMRs when averaged over sufficient time (i.e., longer than the typical libration period), primarily via numerical simulations but also supported by an analytical framework from Hadden (2019). While we focused on second-order resonances in this work, as they are the widest and strongest resonances of this shape and therefore the likeliest to have a detectable signal, we note that this type of narrow, time-averaged peak shaping is expected to be present for higher-order resonances as well.
We fit a model to the 5:3 and 3:1 peaks in the Kepler period ratio distribution, finding a small peak at the 5:3 and no evidence of a peak at the 3:1.
The lack of these peaks in the planet pairs observed by Kepler therefore places an upper limit on the eccentricity distribution of these planets. We place that upper limit at σ = 0.245 (3:1) and σ = 0.095 (5:3) at 95% confidence for a Rayleigh distribution with scale σ for each planet.
Using our model of the 5:3 Kepler data, we can calculate a resonance probability for the planet pairs located near the resonance. For the K02261, a potentially resonant planet pair identified by this calculation, we investigate the observed transit times from Rowe et al. (2015) and conclude that it is possible to estimate the change in period ratio over time based on the transit times. However, in this case the observational errors in the TTVs are too large to place any constraints on the individual system.
Given that not all of the highly eccentric planets in our numerical simulations are ruled out by stability, we infer that some other mechanism, likely from formation processes, drives the observed low eccentricities of the planet pairs observed by Kepler.
Our method is an independent approach to constraining the eccentricities of the population observed by Kepler, and we find results consistent with other methods: the small planets typical of Kepler observations tend to have low eccentricities, lower than those typically observed of giant planets.

ACKNOWLEDGMENTS
This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. We thank the anonymous referee for comments that greatly improved this work. We thank Darin Ragozzine for a helpful discussion on quantifying our model via equation 2 and Sam Hadden for assistance in using the celmech package and his analytical model. ), shown in the top panels as the combined eccentricity (a) and period ratio (b) distributions for these two halves of the population. In the bottom panels, we split the population into two halves by the combined eccentricity and show the log mass ratio (c) and period ratio (d) distributions for these two halves. While both mass and eccentricity affect the period ratio distribution, the time-averaged peak at resonance is still present even in the low-mass sample, but it disappears in the low-eccentricity sample, indicating the eccentricity is more strongly contributing to the effect. . The probability of the Kepler planet pairs being part of an excess peak associated with 5:3 mean-motion resonance and not part of the background continuum of period ratios. All pairs with a non-zero median probability are plotted; see Table 3 for planet names and values.

Software:
REBOUND 2 , Forecaster 3 , emcee 4 , celmech 5 , SciPy 6 APPENDIX A. KOIS The KOIs used in the analysis of the Kepler population in this paper are listed in Table 4.

B. POPULATION COMPARISONS
In this section, we illustrate various properties of the synthetic population of planet pairs that we generated as described in Section 2.1. When applicable, the corresponding properties of the Kepler sample are shown as well to highlight the general similarity in parameter coverage. Data for the Kepler sample, including errors, were taken from the NASA Exoplanet Archive for all candidate and confirmed KOIs on October 28, 2020. This similarity is important to ensure that we are marginalizing over other contributing factors to isolate the eccentricity effect on the period ratio distribution near second-order MMRs. Figure 19 shows the cumulative distribution of the logarithm of the planet-to-star mass ratios for the generated planet pairs and Kepler planet pairs. The planet masses for the Kepler planets were forecast from their radii using the Forecaster package. Observational uncertainty and the probabilistic nature of the Forecaster mass-radius relation were incorporated numerically by drawing 10 radius and stellar mass samples for each planet. Figure 20 shows the distribution of the inner planet's period for the generated planet pairs and Kepler multi-planet systems. This period generally sets the dynamical timescale of the systems, including tidal circularization and the period ratio libration periods. Figure 21 shows the distribution of the period ratio for the generated planet pairs before and after the stability run. More systems went unstable near to but not right at the exact commensurabilities, leading to some peaking even in the instantaneous (osculating) period ratios. This effect is intensified when time-averaged (see Figure 1). Table 4. KOIs in the vicinity of the 3:1 and 5:3 resonances. KOI data were taken from the NASA Exoplanet Archive for all candidate and confirmed KOIs on October 28, 2020. Observational errors are not listed for the periods (and therefore period ratios) as they are generally of an order of one part per million and may not be representative of the accuracy of the listed periods (i.e. Lissauer et al. 2020 Figure 22 shows the distribution of the mutual inclination for the generated planet pairs before and after the stability run. There is little difference from stability sculpting, and we do not anticipate mutual inclination to play a significant role in this resonance behavior. The mutual inclination does affect the ability to detect both planets via transit; however, as we are modeling the observed population and not the intrinsic population, this effect is not examined in this work. Figure 23 shows the distribution of the magnitude of the difference in the longitudes of the pericenter (| 1 − 2 |, such that 0 is aligned and π is anti-aligned) for the generated planet pairs, and Figure 24 the same for the arguments of the pericenter (|ω 1 − ω 2 |). There is little indication of stability sculpting in these distributions, except perhaps a slight preference for aligned arguments of the pericenter in stable pairs near the 5:3 resonance. The effect of the alignment on the dynamical behavior on non-secular timescales is expected to be small.
The distribution of the combined eccentricities for the generated planet pairs can be found in text in Figure 2.  Figure 16. The period ratio over time for K02261 estimated from transit time data from Rowe et al. (2015). Black points are calculated from the given transit times, and errors are the standard deviation of 1000 transit time realizations drawn using the uncertainty on the TTVs. The uncertainty from observations is too large to reveal any structure in the behavior of the period ratio over time.  Figure 17. Comparison of the measured period ratio over time using exact osculating periods from an N-body simulation and the estimated osculating period ratio using the previous two transits at a given point in time to estimate the planet periods. The Pearson correlation coefficient is 0.834, indicating a highly linear correlation. Exact correlation is shown by the black line for reference. The comparison is done for twenty realizations of the K02261 that are consistent with the observed data, ten which show signs of libration in the resonant argument and ten that do not. Contours of the top 1% and 10% of the point density are shown to aid in interpreting the densest regions. The high correlation indicates that this estimation approach is a useful measurement of the period ratio changing over time directly from the transit times.  final initial all initial Figure 24. Distribution of the instantaneous magnitude of the difference in the arguments of the pericenter (such that 0 is aligned and π is anti-aligned) of the final selected planet pairs at the beginning (dashed black) and end (solid green) of the stability run. The gray background shows the initial magnitude of the difference in the arguments of the pericenter for all planet pairs, including those that went unstable and were excluded from the final selection. Distributions are shown for the (a) 3:1 and (b) 5:3 resonances.