Kepler Multitransiting System Physical Properties and Impact Parameter Variations

We fit a dynamical model to Kepler systems that contain four or more transiting planets using the analytic method AnalyticLC and obtain physical and orbital parameters for 101 planets in 23 systems, of which 95 are of mass significance better than 3σ, and 46 are without previously reported mass constraints or upper limits. In addition, we compile a list of 71 Kepler objects of interest that display significant transit impact parameter variations (TbVs), complementing our previously published work on two- and three-transiting-planet systems. Together, these works include the detection of significant TbV signals of 130 planets, which is, to our knowledge, the largest catalog of this type to date. The results indicate that the typical detectable TbV rate in the Kepler population is of order 10−2 yr−1 and that rapid TbV rates (≳0.05 yr−1) are observed only in systems that contain a transiting planet with an orbital period less than ∼20 days. The observed TbV rates are only weakly correlated with orbital period within Kepler’s ≲100-day-period planets. If this extends to longer periods, it implies a limit on the utility of the transit technique for long-period planets. The TbVs we find may not be detectable in direct impact parameter measurements, but rather are inferred from the full dynamics of the system, encoded in all types of transit variations. Finally, we find evidence that the mutual inclination distribution is qualitatively consistent with the previously suggested angular momentum deficit model using an independent approach.


Introduction
The Kepler space telescope (Borucki et al. 2010), utilizing the transit method, has brought a major increase in the number of known exoplanets.Following the path of Kepler, the number of transit detections is rapidly increasing with the current Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2010) and the future PLAnetary Transits and Oscillations of stars (PLATO; Rauer et al. 2014) missions.
The vast amount of photometric data calls for interpretation methods that enable us to extract the most out of the data accurately and efficiently.One of the most direct ways to extract planetary physical and dynamical information from photometric (and other types of) observational data is to run multiple N-body integrations.Such integrators were implemented at high accuracy (e.g., Wisdom & Holman 1991;Chambers 1999;Deck et al. 2014).The drawback of using N-body integrations is the relatively high computational cost, which increases with the time span of the data, a consideration that becomes more and more significant as the available data increase and is now frequently spanning multiple decades.In order to enable a more efficient interpretation of photometric data, different authors approached the interpretation process in an analytic manner, usually using transit timing variations (TTVs; Agol et al. 2005).When TTVs are detectable, the orbital configuration and the planetary masses may be estimated.This made TTVs a tool of great utility and thus a topic of thorough investigation, a few examples being Nesvorný & Morbidelli (2008), Nesvorný (2009), and Nesvorný & Beaugé (2010), a series of papers that analyzed the TTV problem and developed the method ttvim; Lithwick et al. (2012), who analyzed the fundamental TTV mode to first order in eccentricity; Deck & Agol (2015), who highlighted the synodic "chopping" TTV; Hadden & Lithwick (2016), who extended the method of Lithwick et al. (2012) to second order in eccentricity; and Agol & Deck (2016), who published the publicly available code TTVFaster that analytically calculates the TTV to first order in eccentricity.All of the above, however, analyze the TTVs-which are only a by-product of the main measured quantity: the stellar flux.There are two main drawbacks to performing a TTV analysis rather than a full light-curve fitting.First, for small planets with shallow transits, an individual transit (unlike the full "folded" light curve) may not be significant in the data, making timing severely limited or impossible, thus preventing TTV analysis.In such cases, a full, global, light-curve model can extract information that cannot be extracted by attempting to fit transit times.Second, restricting the analysis to TTV ignores other types of transit variations, such as transit duration variations (TDVs).These are of high interest in probing forces out of the plane but are unfortunately even more difficult to measure than TTVs on long-cadence data such as Kepler's; the utility of such information, in the rare cases where it is observable, can be seen in the detection of the mutual inclination between planets b and c in the Kepler-108 system (Mills & Fabrycky 2017).This mutual inclination was later explained to have potentially arisen from a resonant interaction with a binary member of the host star, called ivection resonance (Xu & Fabrycky 2019); thus, in this case, the detection of TDVs eventually led to the understanding of a formerly unrecognized dynamical phenomenon.
In a previous study (Judkovsky et al. 2022a, hereafter Paper I) we introduced AnalyticLC, a method and code implementation for modeling light curves.This analytic method balances accuracy and efficiency and can also model radial velocity (RV) and astrometry to enable joint fitting of these data types.We used this method to analyze a set of twoand three-transiting-planet systems (Judkovsky et al. 2022b, hereafter Paper II).This work led to 140 new planetary masses and orbital properties, along with a list of Kepler objects of interest (KOIs) that display transit impact parameter variations (TbVs).Statistics of such transit variations are important for characterizing the distribution of planetary system properties, offering opportunities to derive inferences on population properties (Fabrycky et al. 2014;Xie et al. 2016;Millholland et al. 2021).
In this work, we extend the study to systems that contain four or more transiting planets.Such systems are interesting owing to the dynamical richness embedded in many planetary interactions involving multiple resonances.Multiplanet interactions can also lead to phenomena that are more complex than direct pairwise interactions.An example of such a scenario was presented in the KOI-500 system by Ford et al. (2011), who mentioned that the similarity between the distances from resonance might result in a dynamical interplay that could be manifested in the observations.In Paper I, we formulated this kind of interaction and denoted it super-mean-motion-resonance (SMMR).This dynamical phenomenon arises from the simultaneous interaction of three planets in a near-resonant chain.In multitransiting planet systems, such effects are more likely to be significant.Here we report three such cases in KOI-707, KOI-1589, and KOI-2038 (see Section 4).
This paper is organized as follows.In Section 2 we review the methods used in this work.In Section 3 we review the main findings from this work in two aspects: planetary masses and statistics of measured TbVs.In Section 4 we refer to each of the analyzed systems, including comparison to former literature results and discussion of the main dynamical features of each system.In Section 5 we conclude and discuss future prospects.Numerical results are provided in tables in the Appendix.

Methods
The methods used in this work are the same ones used in Paper II, with a slight adjustment: the number of walkers in the DE-MCzs process (Differential Evolution Monte Carlo with a sampling of past states using snooker updates; ter Braak & Vrugt 2008) is larger because the number of planets, as well as the corresponding number of model parameters, is larger.We repeat the main stages below.

Data Reduction
The detrended light curve is obtained by applying a cosine filter to the Presearch Data Conditioning (PDC) maximum a posteriori (MAP) Kepler data (Stumpe et al. 2014).The cosine filter is based on the one used by Mazeh & Faigler (2010) and modified such that the shortest timescale is four times the transit duration time, as provided by NASA Exoplanets Science Institute3 (NExScI).The aim of this choice is to avoid the filter affecting the transit shape itself.Outliers beyond 5σ were iteratively rejected (outliers were removed separately for Kepler long-and short-cadence data).
The raw data detrending and the model fitting affect each other, so we performed an iterative process: having detrended the light curve, we fitted a model, normalized the raw data by the best-fitting model, and then detrended again.This iterative process was regarded as converged when the χ 2 score did not improve from the previous iteration by more than the equivalent of 1σ, which is defined as the 68th percentile of the χ 2 (N) distribution, where N is the number of model degrees of freedom.This process has already been proven robust in Paper II, and its robustness was reaffirmed here.
We used all available Kepler data (long and short cadence).Binning was applied for the long-cadence data points, where the number of binning points was set by requiring that the error caused by the binning is an order of magnitude smaller than the typical data uncertainty, using the formulae given by Kipping (2010).
The model function to evaluate the instantaneous flux was AnalyticLC, which was described in Paper I and was already used to study Kepler's two-and three-transiting-planet systems in Paper II.The model uses an analytic 3D dynamical model coupled with the flux calculation formulae of Mandel & Agol (2002).For flux modeling, we took the quadratic limbdarkening parameters from NExScI.

Model Parameters
The model includes eight parameters per transiting planet.Two relate to the physical properties of the planet: μ (planet-tostar mass ratio) and R p /R * (planet-to-star radius ratio).The other six are (equivalent to) the orbital elements: P (orbital period); T mid0 (reference time of midtransit); Δe x , Δe y (eccentricity vector component differences from the neighboring interior planet; for the innermost planet these are just the eccentricity vector components); I I I I cos , sin x y = W = W (inclination vector components, where I is the inclination with respect to the z-axis and Ω is the longitude of ascending node).The axis system is defined such that the x-axis points from the center of the star to the observer, y lies on the sky plane such that the innermost planet crosses the sky plane on it, and the zaxis lies on the sky plane and perpendicular to the x-and yaxes.Using this axis system, e ecos is the eccentricity component on the line of sight and is positive if the periapse is between the star and the observer, where e is the eccentricity magnitude and ϖ is the longitude of periapse.A diagram illustrating the orbital elements and the axis system is shown in Figure 1 of Paper II.
For the innermost planet, we fix Ω = π/2, i.e., we set I x = 0.An additional model parameter is the ratio of the innermost planet semimajor axis to stellar radius, a/R * (as in Ofir & Dreizler 2013), leaving it still with eight parameters.Using Kepler's law, the semimajor axis of the other planets in the system is then scaled from the innermost one.

Nonlinear Fitting Method
We use our code implementation of DE-MCzs (ter Braak & Vrugt 2008) to generate a posterior distribution and obtain a best-fitting model.This method suits multidimensional problems with correlated parameters, as in our case of fitting a dynamical model to a light curve.
We limited the eccentricity component magnitude to smaller than 0.6 and the inclination component magnitude to larger than 50°since the model is usually invalid beyond these relatively extreme configurations; in practice, the solutions usually converged to much smaller values of a few percent/a few degrees.
The choice of the number of walkers was guided by the value of N N log , where N is the number of fitted parameters.There is no general prescription for the optimal number of walkers; we based this choice on experimenting with different types of nonlinear problems and on its satisfactory performance in our previous work on two-and three-transiting-planet systems (Paper II).Convergence was declared upon reaching a Gelman-Rubin criterion (Gelman & Rubin 1992;Brooks & Gelman 1998) of less than 1.2 for all parameters over 10,000 generations.After the removal of the burn-in, we were left with a well-mixed sample, from which statistical inference of the parameters was made.The parameter values described in the tables in the Appendix are the sample median and the difference between the median and the 15.865th and 84.135th percentiles, corresponding to 1σ in a normal distribution.
The walkers' initial states in orbital periods, reference times of midtransit and planet-to-star radius ratio, and planet-to-star mass ratio are distributed based on values and errors from NExScI.If no planetary mass is provided, we translate the archive planetary radius to an initial guess on mass by using an empirical formula given by Weiss et al. (2017, Figure 8).
We ran the nonlinear fitting procedure five times per system in order to check whether it converges to the same minimum consistently, using a different realization of the initial walkers' distribution at each run.In our former study of two-and threetransiting-planet systems, we empirically found that this number suffices for most cases for obtaining repeatable results.We also note that each run contains 112 walkers (for fourplanet systems), 144 walkers (for five-planet systems), and 184 for six-planet systems.Multiplied by tens of thousands of generations at each run, the generated samples covered a large fraction of the parameters space.In Section 4 we discuss the results system by system and note cases with ambiguous solutions that require further investigation.
The initial guess of a (1) /R * was obtained from the literature stellar mass and planetary orbital period; whenever possible (for the decisive majority of planets), we used stellar data from Berger et al. (2020).For KOI-593 and KOI-1422, where this source had no data, we used NExScI.
The fit parameters of AnalyticLC are dimensionless; we translate them to absolute masses and radii using stellar literature data from Berger et al. (2020).

Verification versus N-body Integration
In order to verify the validity of the best-fitting analytic model generated by AnalyticLC, which is calculated using truncated series expansion of the full 3D gravitational interaction, we compare it to a model generated by an N-body integration using Mercury6 (Chambers 1999).Following Ofir et al. (2018), the data uncertainties normalize the flux differences between the analytic model and the N-body model, and the sum of the squares of these normalized differences provides a χ 2 -like estimate of the mismatch between the analytic and the N-body model at the best-fitting parameter set: where F ALC (t) and F Nbody (t) are the flux values obtained from AnalyticLC and the N-body integration at each data point, respectively, and dF(t) is the data uncertainty of each data point.Nbody 2 c measures the ability to statistically discern between the two models given the data available.
We used the empirical cumulative distribution function (CDF) obtained from the posterior distribution to translate the value to CDF by interpolation and then to the equivalent number of standard deviations assuming a normal distribution.We refer to this as σ Nbody , which quantifies the systematic error of the model relative to the statistical error of the data.
Figure 1 shows an example of the comparison between AnalyticLC and the results of a full N-body integration for the four-transiting-planet Kepler-79 system (KOI-152).The best-fit solution of AnalyticLC for this system agrees well with an N-body integration (σ Nbody = 0.0321), implying that the model error with respect to full N-body integration is two orders of magnitude smaller than the statistical error arising from the data uncertainty.By performing the computationally heavy multidimensional search using AnalyticLC, while validating only the best-fit point using an N-body integrator, we ensure both the efficiency of our search and the accuracy of our result.
In Paper II we have shown a map of the validity of the models generated by AnalyticLC.Here we extend this map to include both the two-and three-transiting-planet systems from that work and the systems from this current work.This map is shown in Figure 2. Systems with eccentricities and inclinations of up to 0.1-0.2 and orbital separations of 10 mutual Hill radii can be adequately modeled by AnalyticLC in the vast majority of cases to the precision of the Kepler data.There are also systems of smaller separations and/or larger eccentricities in which the model is adequately correct.Systems that include close approaches closer than ∼7R H -8R H are usually not consistent with the N-body-integrator-based result.This is an expected outcome of our analytic approach, which inherently assumes slow variations of all orbital elements except the mean longitudes-an assumption that breaks when close encounters occur.
It is noteworthy that our analysis is a type of model fitting, and therefore it has some inherent limitations.For example, it will not model any real part of the system that was not included in the model, e.g., nontransiting planets.It is clear from Figure 2 that AnalyticLC found that some systems have relatively high eccentricities-much higher than is typical for multiplanet systems.We believe that at least some of these systems actually contain additional, nontransiting planets that affect the dynamics of the system, and the observed eccentricities are simply a reflection of AnalyticLCʼs attempt to use the degrees of freedom open to it to match the observed data.
Regarding the current sample of systems, from the 23 systems that passed all our tests (see below) and attain σ Nbody < 1.5, for 16 systems 0.1 Nbody 2 s < , implying a very good match (such that when summing in quadrature the systematic and statistical errors, the systematic error contribution is an order of magnitude smaller than the statistical one).

Comparison to Strictly Periodic Circular Model
Our model has eight parameters per planet, while a minimal model of a circular strictly periodic orbit would involve only five for the innermost planet (P, T mid0 , R p /R * , a (1) /R * , b, where b is the impact parameter) and only four for the other planets (a/R * for the other planets can be deduced from the period ratios).In order to justify the larger number of parameters in the full model, we calculated the Akaike information criterion (AIC; Akaike 1974) for both the dynamical and strictly periodic models, and we validate only solutions where the AIC of the dynamical model is better.A more detailed discussion of the selection of AIC as our information criterion is given in Paper II.

Solution Consistency with Stellar Parameters
One of our model parameters is a (1) /R * , the ratio between the semimajor axis of the innermost planet and the stellar radius.This parameter can be obtained alternatively from literature values for the stellar mass and radius (Fulton & Petigura 2018;Berger et al. 2020) and the measured orbital period of the innermost planet.We checked the consistency of our converged solution with the literature-obtained value with respect to the literature error and our posterior distribution error on a (1) /R * .In the systems for which we provide a dynamical solution, there are no cases in which this difference is larger than 2.7σ, where σ is the root sum square of the error estimates on a (1) /R * from our posterior distribution and the literature.In fact, the differences between our estimated a (1) /R * and literature values distribute nicely around zero with a standard deviation close to unity, resembling a normal distribution.We therefore consider our results to be in agreement with the literature's stellar parameters.

Consistency among Solutions
As described above, we ran DE-MCzs five times for each planetary system, each time with a different random seed.We discarded solutions for which our model is not compatible with a full N-body integration (see Section 2.4) and solutions for which the number of model parameters is statistically unjustified (see Section 2.5).In addition, we discarded solutions that imply unreasonably high planetary densitymore than 2σ above 12 g cm −3 , the approximate density at the base of Earth's outer core (Sorokhtin et al. 2011, chap. 2).This process left us with a subset of runs for each KOI.If all of them converged to the same maximum likelihood region in parameter space, the solutions are consistent with each other, and we report the obtained solution posterior distribution statistics.If they converged to different regions in parameter space, we report all of them and regard one of them as the "adopted" solution, recognizing that this solution is not unique.Selection of the adopted solution is based on differences in fit quality (χ 2 ).If there are a few solutions of similar fit quality, we publish all of them and select one of them as the adopted solution based on physical reasoning, e.g., favoring solutions with plausible planetary densities, with small eccentricities for short-orbit planets, and with smaller mutual inclinations, as we treat here systems with four or more transiting planets, which are statistically likely to possess small inclinations (Ragozzine & Holman 2010).For systems with more than one solution, we refer to the various solutions in Section 4.
For 11 out of the 23 systems with valid dynamical solutions, we report one solution, either because it was more likely than the other solutions by more than ∼3σ (e.g., , or because only one solution passed all our validity criteria (e.g., KOI-720).In some cases, a few runs converged to the same local minimum, yielding consistent solutions; in such cases, we report it once.For 12 systems, we report more than one solution (five with two solutions, seven with more).
In order to test the solutions' consistency, for each pair of runs we check the overlapping area between the probability density function (pdf) of the posterior marginal distributions of each parameter.This overlapping area, denoted as η, is a measure of the similarity of two distributions (Pastore & Calcagnì 2019), defined for two pdf's of the parameter x, f 1 (x) and f 2 (x), as where the integration is restricted to the range of validity of the variable x.A value of η = 1 means that f 1 , f 2 are the same distributions.To calibrate the values of this quantity, two Gaussians with the same standard deviation shifted from one another by one standard deviation yield η ; 0.62; two Gaussians shifted by two standard deviations yield η ; 0.31.Any value much smaller than that would mean that the distributions are significantly different; in cases where at least one of the parameters η  0.31 we publish both solutions involved (as long as they have a similar fit quality and pass all the criteria detailed above).This is the case, for example, for KOI-1336.

Sample of Fitted Planetary Systems
Our former work (Paper II) analyzed the data of two-and three-transiting-planet systems.Here we selected systems of four or more transiting planets of Kepler.
Our initial sample consisted of 56 four-planet systems, 17 five-planet systems, and 1 six-planet system (KOI-157).This sums up to 74 systems for which we began the computation process.For 67 of these systems, at least one run converged to a final solution.The validity criteria detailed above were applied to these solutions, namely, (i) N-body matching (illustrated in Figure 1), (ii) AIC improvement over a strictly periodic model (Section 2.5), and (iii) plausible planetary density (Section 2.7).This process filtered out some of these systems.In the end, we were left with a sample of 23 systems (15 of four planets, 7 of five planets, and 1 of six planets), consisting of 101 planets with a dynamical solution, for which we reported masses and orbital elements.
We combine the error obtained from our posterior distribution of the planet-to-star mass ratio and the literature error on stellar mass in quadrature to obtain the total error on the absolute planetary mass.Similarly, we combine the error obtained from our posterior distribution of the planet-to-star radius ratio and the literature error on stellar radius in quadrature to obtain the total error in absolute planetary radius.In most cases, the majority of the error in mass stems from the fit and not the stellar parameters.For the radius, it is the other way around: the radius ratio is well constrained by the fit, and most of the uncertainty arises from the uncertainty in stellar absolute radius.

General
Of the 23 systems for which we found a valid solution, 15 contain four transiting planets, 7 contain five transiting planets, and 1 contains six transiting planets (KOI-157, Kepler-11), summing to 101 planets in total.In Figure 3 we show a map of orbital periods, planetary radii and planetary densities obtained in this work, and resonance locations.This map gives a brief overview of the planetary systems for which physical properties were obtained in this work.Many systems are near at least one resonance location, as the proximity to resonances generates the high-amplitude TTVs that enable mass estimation.The uniform spacing and the similar planetary radii are visually apparent in many systems (the so-called "peas in a pod"; Weiss et al. 2018).In some of these systems, this similarity also seems to occur in planetary masses and densities; this similarity has been recently studied quantitatively and has been shown to exist, though to a somewhat smaller extent than the radius similarity (Otegi et al. 2022).

New Constraints on Planetary Masses
In the left panel of Figure 4 we show the mass-radius spread of the 101 planets in these 23 systems, with curves of constant density and the one-dimensional marginal distribution of this planet population.The radius gap at R p ∼ 1.8R ⊕ (Fulton et al. 2017) is visually apparent; the typical planetary mass is 5-10 m ⊕ .The sample is dominated by planets with densities of 1-3 g cm −3 .While these densities are definitely not consistent with a mostly rocky composition, they are also consistent with the solar system's ice giants only at the lower end of this range.This suggests that the typical composition of these exoplanets is unlike anything An overview of the N-body match with AnalyticLC for all our runs, including both this work and the former work of the two-and threetransiting-planet systems, including runs that both did and did not pass our acceptance criteria.This map demonstrates the accuracy limits of Analy-ticLC.The horizontal axis is the largest magnitude of all free eccentricities and inclinations in the system at the best-fitting solution derived by AnalyticLC.The vertical axis is the minimal closest approach in the system in units of mutual Hill radii.The color of the points represents the value of σ Nbody , with points at which σ Nbody > 1.5 (which do not pass our acceptance criterion) also marked with black edges.
found in the solar system, with the total mass more evenly divided between the rocky core and volatile envelope.Some planets have densities lower than 1 g cm −3 -these are planets massive enough to keep their large gas atmosphere, and they are not extremely hot (none of them has an orbital period shorter than 10 days).The uppermost point in this panel belongs to KOI-191.01; the derived mass and radius imply a density of ∼0.05 g cm −3 ; as we discuss in Section 4, the inference on this planet is questionable and deserves further study.The next low-density planet belongs to KOI-152.01, with a planetary mass of ∼9.4 m ⊕ and a radius of ∼7.2 R ⊕ ; this value is more reliable and is consistent with past literature to ∼1.5σ (Jontof-Hutter et al. 2014;Hadden & Lithwick 2017;Jontof-Hutter et al. 2021).The curve of 5.5 g cm −3 was chosen because it is the approximate Earth's bulk density; the curve of 12 g cm −3 was chosen because it is the estimated approximate density at the base of Earth's outer core (Sorokhtin et al. 2011, chapter 2), and as such it is used to delineate the upper limit of accepted solutions (Section 2.7).
In the right panel of the same figure, we show a comparison of the masses obtained in this work with values of HL17 and JH21, which are studies that have a large number of common KOIs with our work.Our results are in good agreement with the masses obtained by JH21 (some of which were given as upper limits; these are indicated by arrows pointing left).
HL17 applied a default mass prior and a high-mass prior and got two results for each planet; the masses we obtained are typically larger than their low-prior mass and are close to their masses obtained from high-mass priors.This is a reasonable outcome, as our prior is uniform in mass, the same as their high-mass prior.
For KOI-834.03, HL17 report a mass of 239.1 m ⊕ , which is out of the bounds of this plot; the value is not plausible, as this is a 1.95 R ⊕ planet candidate.For six planets, the results we obtained are different from the results obtained by HL17 by more than 2σ.All of these cases are indicated in Figure 4 specifying their KOI numbers: KOI-232.04,KOI-232.05,KOI-834.01,KOI-157.03,KOI-707.04, and KOI-1589.03.For KOI-707.04 our estimate is only 1 m ⊕ below the upper limit given by HL17; for KOI-232.05we might overestimate the planetary density, and the value given by HL17 seems more physically plausible; however, it is significant to less than only 3σ.In the other cases we believe that the solutions we provide are more physically plausible than the ones given in HL17.We refer in more detail to all these specific cases in Section 4.
Overall, there is good agreement between the masses obtained in this work and previously reported planetary masses, thus giving confidence in the fitting process and the reliability of the newly reported masses.The masses, radii, and orbital elements of our adopted solution are tabulated in the Appendix, along with machine-readable files.
In the left panel of Figure 5 we show the distribution of fractional error in planetary mass, σ m /m, for all the planets for which we obtained a valid dynamical solution in this work (including the 140 planetary masses obtained in Paper II).The typical fractional error in planetary mass derived from this study is better than the typical fractional error from past studies by roughly a factor of two.We attribute this improvement to the global fitting approach, which integrates all of the transits together and takes into account all types of transit variations, rather than using TTV fit only (the technique used in most of the referred studies).In addition, in the right panel of the same figure, we show that most of the new planetary masses we obtained are of planets with small TTV magnitudes, quantified by the TTV standard deviation σ TTV , calculated over the TTV values of all transit events in the best-fitting model.The current study evidently includes more TTVs of lower amplitude, and we believe that the explanation for this is that systems with weak dynamical interaction (and consequently small-amplitude TTVs) are usually of small planets, whose individual transit times are anyway poorly constrained.The combination of shallow transits with weak interactions makes it difficult, and sometimes impossible, to extract individual transit times and detect a TTV pattern; in such systems, the global flux fitting Each blue error bar is related to a single planet.The red lines are normalized histograms of the masses and radii obtained by summing up the pdf's of all points.As the sample is small and radii are well determined, the radius-weighted histogram appeared more as a collection of discrete values, so it was smoothed using a Gaussian kernel of width 0.3 R ⊕ to better show population-wide trends.The gray contours are constant-density curves.(b) Comparison of masses obtained in this work with literature values: JH21 values (blue) and upper limits (red); HL17 default mass prior, three of which show upper limits only (yellow); and high-mass prior (purple).We omit planets for which literature values are within less than 2σ of zero.For six objects labeled with their KOI numbers, the results obtained in this work disagree with the results of HL17 by more than 2σ; these are discussed in more detail in the text.The black line shows the identity function.
can exploit the data better to obtain an estimate of planetary mass.

Out-of-plane Forces
Forces out of the plane are manifested in transit light curves as TDVs, or, alternatively, TbVs.In Paper II we have shown how TbVs can be interpreted as arising from the interaction among the transiting planets, or from an interaction with a nontransiting companion.
Only in a handful of cases was mutual inclination constrained from transit variations; however, such cases are of high interest and exhibit dynamical richness.A good example is the Kepler-108 system, in which the detection of TDVs led to a photodynamical analysis that showed that the planets in this system are highly inclined (Mills & Fabrycky 2017).Later on, the appearance of such inclination was explained as arising from a binary member via an interaction denoted as ivection resonance (Xu & Fabrycky 2019).This term refers to a situation in which a resonance between the binary orbital period and the planetary nodal motion acts to pump the planetary orbital inclination.The example of Kepler-108 shows how the detection of a mutually inclined system led to the discovery of a new dynamical phenomenon.
Another avenue in researching mutual inclination in planetary systems involves population analysis of the detected transit variations, such as the usage of the TDV catalog of Shahaf et al. (2021) by Millholland et al. (2021).The latter authors addressed the distribution of mutual inclinations in planetary systems by analyzing the number of expected TDVs in the case of a dichotomic model that assumes a bimodal distribution in mutual inclination and a model that assumes that mutual inclinations can be described by a distribution arising from a model based on angular momentum deficit (AMD; Laskar & Petit 2017;He et al. 2020).The results of Millholland et al. (2021) have shown that in terms of the number of detected TDVs the AMD-based model shows better consistency with the Kepler population than the dichotomic model; this is a significant conclusion, as it offers a nondichotomic solution to the so-called "Kepler dichotomy."This exemplifies the importance of population analysis of transit variations.It is therefore a natural motivation for us to try to expand the knowledge we have on transit variations arising from forces out of the plane.
In Paper I, we produced a catalog of TbVs of planets in twoand three-transiting-planet systems; here we produce a similar catalog for systems having four or more transiting planets (Table 4).The current catalog contains 71 planets undergoing TbVs to better than 2σ, out of which 52 are better than 3σ; together with the former catalog from Paper II, we present 130 planets undergoing TbVs to better than 2σ (77 to better than 3σ), an extension to the catalog of TDVs published by Shahaf et al. (2021) that contains 31 KOIs that display statistically significant TDVs (the only KOI of multiplicity larger than 3 included in their catalog is KOI-841.02,which is also included in our catalog as displaying TbVs).
We note that for most analyzed systems the signs of the impact parameters are not definitive, despite the typical small errors on b.This is because |b| may be well constrained by the shape of the transit, even if the dynamics do not distinguish between the same-hemisphere configuration (positive b) and the opposite-hemisphere one (in the context of impact parameters and orbital alignment; see also Fabrycky et al. 2014).We tested this using the following procedure.We examined all of our valid solutions and switched the signs of the impact parameters in all possible configurations.In the language of our fitted parameters, we switched the signs of either I x or I y or both for each planet, excluding the innermost planet (for which I x is set to zero and I y is limited to positive values only, without the loss of generality).This resulted in 4 N 1 plconfigurations for each system, where N pl is the number of planets in the system.For each system, we checked the χ 2 value difference with respect to the original best-fitting solution.Only in two systems do the best-fitting solutions stand out among the others at the level of 2σ: the solutions for KOI-834 and KOI-707.This shows that the signs of I x and I y (and with them the sign of b) may be inverted in most cases.The two systems mentioned above merit further study, as our solutions constrain their full 3D orbital geometry.
We note an essential difference in the estimation method between our catalog and the catalog of Shahaf et al. (2021).Shahaf et al. (2021) used individual transit duration measurements and sought significant trends in the duration values for specific KOIs.In other words, they searched for TDVs that could be observed within the time span of the Kepler mission.In this work, we calculate the linear trend in the impact parameter per KOI using the orbital elements' evolution in time, elements that affect all transit variations combined.For many KOIs, we predict a finite TbV trend, although this trend may not be directly seen in the impact parameter measurements.Therefore, our TbV catalog is not a direct empirical observation of TbVs, but a projection of the orbital dynamics on the TbVs.This approach enables predicting the future dynamics of the systems in general: not only discovering faint, slowly varying, and periodic TbVs in existing data, but also predicting the detectability of such TbVs not detectable in the time span of the Kepler data.
Having said all that, the results here represent the largest compilation of planets with such variations to date.It is expected that future population analyses, similar to Millholland et al. (2021), would benefit from this expanded catalog of TbVs.
Incorporating forces out of the plane in the dynamical model has an additional contribution to the masses' determination.The traditional method of using the TTV only can yield mass estimates, but these are sometimes degenerate with eccentricity (Lithwick et al. 2012).Breaking this degeneracy requires either significant eccentricities that would generate a phase shift between the planets' TTVs or additional effects, such as the synodic chopping effect (Deck & Agol 2015) or higher-order TTVs (Hadden & Lithwick 2016).In Figure 6, we show an example of a nonlinear TbV pattern in addition to the linear secularly driven TbV drift.We define b resid to be the residual impact parameter variations after the removal of the best-fitting linear trend.The residual TbVs are approximately periodic, with a period that is the TTV-derived super-period, indicating that these are TbV manifestations of the same planet-planet interaction seen in the TTVs.Since AnalyticLC takes into account the various transit variations simultaneously, no additional mechanism is needed to capture this variability.The TbV information adds to the constraining power of the global light-curve approach.In panel (c) of the same figure, we show the distribution of b resid s , which is defined as the standard deviation of b resid over time of the best-fitting parameters for the 241 planets included in our adopted dynamical solutions in both this work and Paper II.It shows that the typical value of b resid (which is related to near-MMR interactions) is 10 −4 to 10 −3 , roughly an order of magnitude smaller than the typical linear change in b (which is related to secular interactions) accumulated in a year.This shows that the light curve contains information content regarding both secular and near-resonant interactions, where the former is more prone to detection for the typical observation time span of a few years.Important for MMR-originating transit variations other than TTV are the TDVs used to extract the properties of the nontransiting planet in the KOI-142 system (Nesvorný et al. 2013).
We now turn to discussing long-term TbV effects.In Figure 7, we show an overview of the 130 values of planets displaying long-term linear TbV with a significance of better than 2σ.In the left panel, we show a scatter of these planets' TbV rates against P min , the orbital period of the innermost planet in each system, which is a probe of the amount of inward migration the system has experienced.We see that rapid TbV rates occur in systems where the innermost planet's orbital period is less than 20 days.Specifically, there are eight planets displaying TbV rates faster than 0.05 yr −1 , seven of which reside in systems that contain planetary companions with less than a 5-day period (the exception is KOI-988.02).
We note that Figure 7 shows the dependence of b  on the minimal observed orbital period, not on the period of the planet undergoing the TbV.Some positive correlation between the duration variation rate T  and the orbital period was proposed by Shahaf et al. (2021, their Figure 5), at least for the significantly constrained T  values.These are related to our measurements of b  but are not directly comparable because the transit duration depends on additional factors such as orbital velocity, and the impact parameter depends on geometry alone.Note that b  and T  arise from different dynamical phenomena: b  arises mostly as a result of out-of-plane secular interactions, while T  can also arise as a result of in-plane apsidal motion that changes the skyprojected velocity of the planet.
Are systems with short-period planets expected to display stronger TbVs than other systems?If TbVs are probes of mutual inclination, then there is some supportive evidence for that.Millholland & Spalding (2020) have shown that the formation of ultra-short-period (USP) planets, which is a result of inward migration, would require some initial mutual inclination in order to enable the conservation of angular momentum along the migration process.Though we do not focus on USP planets here, the same argument implies that systems with short-period planets, which likely migrated to this configuration, would require some initial mutual inclination within the system.A more quantitative, dynamical analysis of the relation between b  and P min is left for future work.
To further investigate the population of planets displaying TbV, we define The dimensionless quantity n b is the characteristic number of orbits required for the planet to change its impact parameter by unity (given a constant change rate), i.e., a currently transiting planet will be torqued out of transit after typically n b transits.In Figure 8 we plot this number against the orbital period of the TbV-displaying planets.We find that n log b varies roughly inversely with P, that is, that the measured b  values are only weakly correlated with the orbital period (panel (b)).If this finding holds also for larger periods (which is not known), then planets of orbital periods of decades would typically be torqued out of transit after just a few orbital periods.This would imply not only that long-period planets have a lower geometric transit probability (Ragozzine & Holman 2010) but also that they tend to move out of transit in a small number of orbits, limiting the utility of the transit method.Further study is required in order to shed light on this point.We resist fitting the n b (P) trend because we find the fit parameters sensitive to the removal of a small number of points from the sample.
In addition to the TbV analysis, we also examine the direct output of the dynamical fitting: the mutual inclinations among the planets.In Figure 9 we show the empirical distribution of mutual inclinations obtained from all our runs that passed all our validity tests, including a summation of all configurations obtained after changing the signs of I x and I y in systems whose fits permitted this inversion (as explained above).This figure shows that the mutual inclination among the transiting planets, i m , distributes in a lognormal shape, indicating that systems with more transiting planets tend to have somewhat smaller mutual inclinations.Both the shape of the distributions and the  6), who proposed that the statistics of planetary system properties can be explained by an AMD-based model (Laskar & Petit 2017).Their proposed model was later supported by Millholland et al. (2021), based on the statistics of transit duration (Shahaf et al. 2021).We note that in many systems, both in this work and in Paper II, we propose the existence of a nontransiting external perturber.This means that the inclination values obtained in those systems may be suspect.Nonetheless, the qualitative similarity to the results presented by He et al. (2020, their Figure 6) is compelling because their method relied on detection statistics of the Kepler population, while ours relies on dynamical modeling in individual systems.The fact that both methods yield similar findings supports the proposed AMD-based model.

Description of Individual Systems
In this section, we provide details for each of the 23 systems with a dynamical solution.We review previous literature on planetary masses and highlight the dynamical features of the system.The systems are ordered by their KOI index.The comparison to former literature masses is given in Figure 10 for all the KOIs.We include both past estimates of the true planetary mass and estimates of planetary nominal mass (Hadden & Lithwick 2014;Xie 2014, hereafter HL14 and X14, respectively), which is, in many cases, an upper  limit on the true mass.For the values of Hadden & Lithwick (2017, hereafter HL17), we use their high-mass prior (which is appropriate for comparing with our results, which are obtained from using the same prior as theirs).
From our sample of 101 planets, for 46 we did not find any previous literature mass value.For 41 out of those, our mass detection is significant to more than 3σ.For 10 planets, we obtained a median mass estimate lower than 2 m ⊕ , of which 6 are lower than 1 m ⊕ : KOI numbers 248.04, 520.04, 571.02, 841.04, 841.05, and2038.03. For three of those (248.04, 520.04, 2038.03)we found the results to be consistent with the upper limits given by HL17 and JH21, and for the other three (841.04,841.05, 571.02) we did not find any previous mass constraints.Out of 101 planets in systems with a valid solution, for 95 we provide mass constraints better than 3σ.
Kepler-79 (KOI-152).-Allour runs converged to solutions of similar masses and to nearly circular orbits.One of the runs yielded a solution significantly better than all others regarding fit quality; it indicates small mutual inclinations of a few  2014), within ∼1.5σ of the results of HL17 in their high-mass prior case, and consistent with the results of JH21.The solution suggests a clear distinction between the densities of the two inner planets (∼1-1.5 g cm −3 ) and the two outer planets (∼0.15, 0.25 g cm −3 ).In each subpanel, the planets in the system are sorted (left to right) by their orbital period.The numbers above are the suffix of the planets' KOI designation; for instance, the first object in the plot is .Dotted vertical lines separate different planets within a system.Bold underlined KOI numbers are used for systems with at least one planet without a previously reported mass constraint, as in Figure 3.The bold underlined KOI suffix represents a specific planet without a previously reported mass.The legend shows the colors and markers representing different literature sources.Arrows pointing down represent upper limits on mass.).-This system, which contains six transiting planets, has been examined in many studies.Lissauer et al. (2011) analyzed the TTVs known at that time to estimate the planetary masses and studied the system's dynamical stability.These masses were then estimated again (Lissauer et al. 2013) using more (14 quarters) of the Kepler data, and later again by Borsato et al. (2014).HL14 estimated the nominal masses and, later (HL17), the true masses of the planets in this system.Our results agree well with these studies, apart from the fourth known planet from the star (Kepler-11 e, KOI-157.03),where our adopted solution suggests a mass of m 12.5 0.9 0.8 -+ Å while Lissauer et al. (2013) suggest m 8 2.1 and HL17 suggest m 7.2 1.0 1.1 -+ Å .The period ratios in this system suggest that each planet TTV is affected by more than one specific frequency; this, along with the relatively massive planets and the wealth of information encoded in six transit signals, makes the solution distinct with multiple studies in agreement.
Kepler-487 (KOI-191).-Wefound no literature mass for the four KOIs in this system.Until recently, the inner two were designated as candidates and the outer two as confirmed.Currently,  is designated as confirmed in NExSci following the analysis of Valizadegan et al. (2022).
This system is of dynamical interest for several reasons.First, of the four transiting planets it contains, the third one (planet b, orbital period of 15.35 days) displays a transit depth that suggests a Jupiter-size object, not a frequent scenario in compact Kepler systems.Second, this system harbors a USP planet (Winn et al. 2018).This type of planet has probably gone through a runaway migration process that requires strong dissipation mechanisms (e.g., Millholland & Spalding 2020, and references therein).Third, unlike many Kepler compact systems, the planets here are not near any resonance.
We provide one solution, with masses of  .We note that this solution involves a large eccentricity of the innermost 0.7-day-orbit planet.Such a situation would be surprising, as short-orbit planets should have their orbits circularized on short timescales.
A possible explanation we suggest is that, although our pipeline converged to a four-planet model that fits the data, the system contains more than four planets.A nontransiting fifth planet could possibly explain the high-amplitude TTV (∼13.4 minutes) of the outermost transiting planet (planet c, KOI-191.04),without the need for high-eccentricity orbits.
We therefore provide the four-planet solution of the system but note that it may require further investigation in order to understand its dynamical structure.For KOI-232.04 (Kepler-122 e ), HL17 give an upper mass limit of 0.6 m ⊕ based on their high-mass prior.Given the radius of order 2.5 R ⊕ , our solution is more physically plausible.For KOI-232.05 (Kepler-122 f) HL17 give m 5.3 2.4 3.6 -+ Å , a mass significantly lower than ours.With a radius of ∼2.3 R ⊕ , our solution implies an abnormally high density.An alternative explanation for the possible overestimate of Kepler-122 f density is that this system harbors another external nontransiting companion.Such a scenario could explain the TTVs of Kepler-122 e without the need for a high mass for Kepler-122 f and at the same time explain the significant TbVs displayed by the planets in this system (see Table 4) without the need for the few-degree mutual inclinations suggested by the dynamical solution (see Table 1).
Kepler-49 (KOI-248).-HL17analyzed the interaction between KOI-248.01 and KOI-248.02(Kepler-49 b and Kepler-49 c, respectively).Taking the joint errors, our solution agrees with HL17 values and differs from the values provided by Jontof-Hutter et al. (2016) by ∼2σ.A surprising result of our best-fitting solution is the implied mutual inclination among the planets.The solution suggests a mutual inclination of the innermost planet from the other planets of ∼25°.This tilt is unexpectedly high but may have a real dynamical origin.This innermost planet is at a period ratio of 2.8 with the second planet, not close to any first-order mean motion resonance (MMR), while the second and third planets are close to the 3:2 MMR.The short orbital period of the innermost planet (2.58 days) approaches the 1-day boundary of the orbits of USP planets (Winn et al. 2018).Such planets are thought to have undergone rapid inward migration, during which angular momentum conservation requires that they possess initial orbital inclination (Millholland & Spalding 2020).We highlight this system, as it might include a planet in the process of becoming a USP planet.Future observations of the orbital period of the innermost planet could reveal whether this is the case.
Kepler-26 (KOI-250).-Theplanets in this system were confirmed based on Steffen et al. (2012), who analyzed their TTVs and have shown that the TTVs are anticorrelated and hence the planets are dynamically interacting.Our solution suggests masses of 3. for KOI-250.02(Kepler-26 c), and 5 0.9 0.8 -+ for KOI-250.04 (Kepler-26 e).The masses of Kepler-26 b and Kepler-26 c agree with those given by Hadden & Lithwick (2016), Jontof-Hutter et al. (2016), and HL17.The masses of Kepler-26 d and Kepler-26 e agree with the upper limits given by JH21.Our solution involves eccentricity of order 0.1 for the innermost 3.54-day planet and inclinations of about 10°; these are relatively high values, and we suggest that further observations would be needed in order to confirm them.
Kepler-169 (KOI-505).-Wefound no literature values for the masses in this five-transiting-planet system.We provide two solutions; although one of them is of better fit quality by ∼2σ, we select the other one as the adopted solution owing to the more plausible density it suggests for the innermost planet.Both solutions suggest small eccentricities of a few percent and inclinations of a few degrees, which are also expressed as 2σ TbV signals of the four innermost planets.The adopted solution suggests a monotonic mass order, from ∼1 m ⊕ of the innermost planet to ∼7 m ⊕ for the outermost planet.
Kepler-176 (KOI-520).-Alladjacent pairs in this fourtransiting-planet system are near the 2:1 MMR.HL17 estimated the masses of KOI-520.01 (Kepler-176 c) and KOI-520.03(Kepler-176 d): m m 3.8 , 5.9 Å for KOI-520.04 (Kepler-176 e), consistent with the results of HL17.The outer planet's mass and density are exceptional; its derived mass is one of the lowest known of any exoplanet, and with a radius of about 1.36 R ⊕ , it possesses a density of 0.44 g cm 0.11 0.11 3 -+ -, a surprisingly low value for a small planet below the radius gap.The eccentricities are a few percent.The innermost and outermost planets display significant TbVs, explained in this solution by mutual inclinations of 8°-9°.Based on these unusual characteristics, we highlight this system for further observations, and we suspect that there is an additional planet(s) in the system that affects the observed dynamics.Such a planet might explain the low mass and eccentric orbit of the outermost planet, KOI-520.04,suggested by our obtained solution.
Kepler-186 (KOI-571).-Wefound no literature values for the masses in this five-transiting-planet system.We provide one solution with significant mass constraints of >3σ for four out of the five planets, with masses of m 3.5 0.   is surprisingly low for a 1.52 R ⊕ planet, as it is difficult to retain a substantial atmosphere with such a low mass.In addition, the solution includes a high eccentricity value of 0.2 for the innermost, 3.88-day-orbit planet, .We hence consider this solution as questionable and conclude that further observations would be required to constrain the parameters of this system.
Kepler-616 (KOI-593).-Wefound no literature mass for this four-planet system.We provide two solutions, which are consistent in most parameters and differ only in some of the eccentricity components.Our adopted solution suggests masses of m 6. Å for KOI-707.02(Kepler-33 f); those masses differ from the results of HL17 by 1.5σ-2σ.This difference might partially arise from the TbVs of the outermost planet, which is not taken into account in a TTVonly analysis.Our solution suggests small eccentricities of a few percent but significant inclinations of up to 10°; an alternative explanation for the observed TbVs in this system could be an external, nontransiting companion.This system contains a fingerprint of the SMMR effect we mentioned in Paper I: planets Kepler-33 d and Kepler-33 e are slightly inside the 3:2 MMR, with a super-period of roughly 393.76 days, and planets Kepler-33 e and Kepler-33 f are slightly inside the 4:3 MMR, with a super-period of approximately 321.7 days.The relative proximity of these super-periods yields a TTV pattern arising from the SMMR effect, for which we estimate an amplitude of a few minutes at a timescale of roughly 4.81 yr.In addition, this system is one of two (along with KOI-834) in which inverting the impact parameter signs of the best-fitting solutions yields a fit quality worse by ∼2σ (see Section 3.3).This gives us confidence that the impact parameter signs are meaningful, and therefore these two systems deserve further investigation regarding their 3D structure.
Kepler-221 (KOI-720).-Forthis system of four transiting planets, we found in the literature a nominal mass for planet c (HL14).Å for KOI-812.03(Kepler-235 e).The solution suggests eccentricities of up to 0.1 and inclinations as high as 10°(though with large error bars).In addition, the solution suggests that the innermost transiting planet KOI-812.01, which is above the radius gap, has a similar mass to the second planet (KOI-812.04);it is a surprising result because with similar masses it would be more likely for KOI-812.01 to lose its atmosphere than for KOI-812.04.Given all of the characteristics of the solutions above, we conclude that more data would be required to constrain the planetary masses in this system.It's entirely possible that there is a nontransiting companion in this system, and by not taking it into account, we overestimated the mass of .-HL17analyzed the interaction among the four outer planets out of the five transiting planets.These four planets (c, d, e, f) form a near-resonant 2:1 chain.We provide different solutions, suggesting masses of 5-10 m ⊕ for the three innermost planets, a giant-planet mass of order 70 m ⊕ for the fourth planet, and a mass of 10-15 m ⊕ for the outermost planet.All solutions suggest small eccentricities of a few percent.Most of them, including the adopted one, indicate inclinations of a few degrees.The high-prior mass values of HL17 for planet c are insignificant; our results for planets d and f are consistent with HL17 up to the uncertainty.Our solution for planet e suggests a higher mass than in the solution of HL17; we find our solution physically plausible, as the planetary radius is about 9 R ⊕ .This system is one of two (along with KOI-707) in which inverting the impact parameter signs of the best-fitting solutions yields a fit quality worse by ∼2σ (see Section 3.3).This gives us confidence that the impact parameter signs are meaningful, and therefore these two systems deserve further investigation regarding their 3D structure.
Kepler-27 (KOI-841).-Twoplanets are confirmed in this five-transiting-planet system (KOI-841.01is Kepler-27 b, KOI-841.02 is Kepler-27 c), while the others are considered candidates.HL17 estimated the masses of planets b and c and the mass of the innermost KOI-841.03.We provide a solution that suggests masses of m 2. The two outermost planets, which are of long periods, are not significantly constrained.The suggested masses of planets b and c are slightly higher than the values given by HL17.The indicated eccentricity values of the innermost planet, of order 0.1, are high for a 6.54-day-orbit planet; hence, this dynamical solution is suspicious.An alternative explanation could be the existence of nontransiting companions in the system, which could also explain the TbV signals of KOI-841.01 (Kepler-27 b), KOI-841.02(Kepler-27 c), and KOI-841.05 without resorting to the significant inclination of KOI-841.05 suggested by the solution.
Kepler-245 (KOI-869).-Wefound in the literature nominal masses for the two outermost planets given by HL14.We provide three dynamical solutions with similar fit quality.The adopted one is selected based on its coplanar structure, likely in a four-transiting-planet system (Ragozzine & Holman 2010).Our adopted solution suggests masses of m 3. Å for KOI-935.03.All of these masses are within 1.5σ of the HL17 masses.All planets display significant impact parameter variations; this is probably the source of the relatively high mutual inclinations suggested by the solutions.A possible alternative explanation for the impact parameter variations would be an external nontransiting companion that applies torques on the transiting planets.
Kepler-58 (KOI-1336).-Steffen et al. (2013) have shown that the TTVs of KOI-1336.01(Kepler-58 b) and KOI-1336.02(Kepler-58 c) are anticorrelated and hence reside in the same system, based on the data available at the time.Later, HL14 estimated the nominal masses of KOI-1336.01(Kepler-58 b) based on the TTV amplitude of KOI-1336.02(Kepler-58 c).HL17 provided masses for those planets with large error bars.We provide three solutions, all with similar fit quality.The adopted solution is selected based on best matching with N-body integration.Our adopted solution suggests masses of m 6. Kepler-84 (KOI-1589).-Thisfive-planet system has been analyzed by HL14 for nominal masses, and later by HL17.We provide two solutions; the adopted one is slightly better regarding fit quality.The solutions agree on the planetary masses, apart from slightly different values for KOI-1589.03(Kepler-84 e).The adopted solution suggests masses of m 5.8 1.7 Å for KOI-1589.05(Kepler-84 f).Our solution agrees with the upper limits given by HL17 on KOI-1589.04and KOI-1589.05and with the values given for planet c.Our values are smaller than those given by HL17 for KOI-1589.01and KOI-1589.03.KOI-1589.02displays a significant TbV signal (∼4σ).This system displays a small-amplitude TTV pattern owing to the SMMR effect we mentioned in Paper I: Kepler-84 b and Kepler-84 c are slightly inside the 3:2 MMR, with a superperiod of roughly 273.12 days, and Kepler-84 c and Kepler-84 e are outside the 2:1 MMR, with a super-period of approximately 211.84 days.The relative proximity between these super-periods yields a TTV pattern at a timescale of roughly 2.587 yr, for which we estimate the amplitude to be of order 1 minute.
Kepler-416 (KOI-1860).-In the literature, we found only a nominal mass estimate of KOI-1860.01(Kepler-416 b), provided by HL14.The four transiting planets in this system form an 8:4:2:1 near-MMR chain.We provide two solutions of similar fit quality, which are consistent in most parameters and differ mainly in the mass of the third planet and its orbital inclination with respect to the other planets.We select the adopted solution as the one that suggests a more coplanar structure and has a slightly better fit quality (by an amount equivalent to ∼0.5σ).Kepler-338 (KOI-1930).-For this four-planet system, we found only nominal masses for planets KOI-1930.KOI-1930.03,where the density of the innermost planet is significantly larger than the densities of the other planets in the system.Kepler-341 (KOI-1952).-Wefound no literature mass estimates for this four-planet system.We provide two solutions that agree on most of the fitted parameters, where the adopted one (chosen by its much better agreement with an N-body integration) suggests masses of m 3. Å for KOI-2038.04(Kepler-85 e), consistent with the masses given by HL17.Our solution suggests a roughly coplanar, circular structure with no significant observed TbVs.In this case, due to the large number of different solutions, our confidence in the designation of a specific solution as adopted is not large; in order to constrain the system parameters, additional data would be required.The planets in this system form a chain of three 3:2 MMRs.Kepler-85 c and Kepler-85 d have a super-period of roughly 130.8 days, and Kepler-85 d and Kepler-85 e have a superperiod of approximately 136.4 days.This proximity of the super-periods gives rise to a nonpairwise TTV pattern, which we called SMMR in Paper I. Such an effect also occurs as a result of the proximity of other super-periods in this system; for example, planets b and d (innermost and third from inside) are affected by their proximity to the 2:1 MMR, with a superperiod of roughly 114.2 days.This super-period is close to the 136.4-day super-period of planets d and e.This proximity of super-periods generates an SMMR with a timescale of approximately 700 days that affects the b-d-e triplet.Our best model estimates a total TTV with amplitude 7-8 minutes owing to these SMMR effects.

Summary and Future Prospects
This work completes our previous one (Paper II), in which we analyzed light curves of two-and three-transiting-planet Kepler systems by fitting a model based on the analytic approach AnalyticLC, described in Paper I. In the current work, we fit a model to light curves of Kepler four-, fiveand six-transiting-planet systems.The model takes into account the 3D interactions and therefore is sensitive to different types of transit variations in addition to typically modeled TTVs (Hadden & Lithwick 2014, 2017;Jontof-Hutter et al. 2021, and others).The usage of the entire light curve simultaneously yields, in many cases, tighter constraints on planetary mass, as shown in Figure 5, which allowed a very high yield of determined masses: of the 101 planets in systems with good solutions, significant mass was determined for 95 of them.This figure also shows that the contribution to new planetary masses in our two works together is most significant for small TTV amplitude planets.This is probably due to the difficulty of obtaining a good fit for systems with weak dynamical interactions that yield small TTVs.The global fit approach, which integrates the dynamical information from the entire light curve (TTVs, MMR-related TbVs, secular effects, SMMR, etc.), succeeds in extracting additional information from the data.
In addition to planetary masses, these works investigated impact parameter variations.Because long-term TbV results are attributed to orbital nodal precession, these signals are probes for mutual inclination within the system.We note here that in addition to these long-term variations, other TbV patterns arise from near-resonant terms and are automatically taken into account in our model; however, in the catalog we compiled here, we specifically address linear-trend TbVs.In Paper II and this work, we compiled a list of TbV-displaying KOIs.We show the general spread of the TbV rates in Figure 7 and find that the typical rate (or, at least, the typical rate we can detect with a significance better than 2σ) is of order 0.01 yr −1 .This limited sample of KOIs shows that rapid TbVs appear only in systems that host a planet orbiting the star with an orbital period of less than ∼20 days; this may be explained by the exchange of angular momentum in the process of planetary migration.This preliminary conclusion requires a larger sample to gain significance and a more detailed quantitative dynamical explanation.
We compute the number of orbits for a planet to cease transiting and find that this quantity typically decreases for increasing orbital periods; the trend is shown in Figure 8.This dimensionless timescale is essential in planning long-term observation campaigns to search for Earth-like planets.
We foresee two types of future continuations of this work.One is extending our knowledge of the characteristics of individual planetary systems by using an analytic model based on AnalyticLC.This work and Paper II concentrated on Kepler photometric data; however, large data sets of long-timespan RV observations are natural candidates for such an analytic approach, which is particularly efficient for long time spans.Another related direction would be to combine Kepler data with those from TESS (Ricker et al. 2010) and the upcoming PLATO observations (Rauer et al. 2014), or RV data and astrometric data from Gaia (Gaia Collaboration et al. 2018).A second research avenue is a statistical analysis of the number of detected TbVs and their magnitudes to better understand the nature of planetary systems as a population.Previously, studies showed that an AMD-based model is more consistent with the number of detected TDVs than a model based on a bimodal distribution of planetary systems (Millholland et al. 2021).Using our significantly updated catalog of TbVs may allow further progress.
In Paper I, we published the method and code implementation AnalyticLC, used in this study as the fitting model.We report here a few improvements to the code, mainly that it now enables fitting a joint model of photometry and RV more robustly than before.We also added the ability to incorporate the stellar quadruple moment in the dynamical calculations.The code is available at https://github.com/yair111/AnalyticLC. Note.The numbers are printed in a short version to allow the table to fit the page; the precise values are given in machine-readable format.All parameters are given in MKS, apart from the epoch for which they are specified, which is given in days, with respect to Barycentric Kepler Julian Date (BKJD).The value of the gravitational constant used in our integrator is G = 6.674279896547770 × 10 −11 .The axis system is such that y-z forms the plane of the sky and x points toward the observer (transit occurs when x > 0).
(This table is available in its entirety in machine-readable form.)

Figure 1 .
Figure 1.An example for the N-body matching test performed on the best-fitting solution of the planetary system of Kepler-79 (KOI-152), demonstrating the ability of AnalyticLC to generate a model consistent with full N-body integration of a four-planet system.The mismatch is quantified as 10 Nbody 2 c ~, (σ Nbody = 0.0321)-a systematic error much smaller than the statistical error arising from the data uncertainty.(a) TTV pattern of the N-body model (o) and the AnalyticLC model (x) for the four planets in this system (b in blue, c in red, d in yellow, e in purple).The symbols are on top of each other at this scale.(b) "Residuals": times of midtransit mismatch between the N-body-generated model and AnalyticLC, which is of order a few seconds for the innermost planet and of order a minute for the three outer planets.(c) Manifestation of the mismatch to terms of relative flux.The mismatch standard deviation within the in-transit points is 0.7, 10, 25.5, and 9 ppm for the four planets correspondingly, innermost to outermost.The typical Kepler short-cadence data (which is almost all of the data) uncertainty for this star (∼1000-1200 ppm).

Figure 2 .
Figure2.An overview of the N-body match with AnalyticLC for all our runs, including both this work and the former work of the two-and threetransiting-planet systems, including runs that both did and did not pass our acceptance criteria.This map demonstrates the accuracy limits of Analy-ticLC.The horizontal axis is the largest magnitude of all free eccentricities and inclinations in the system at the best-fitting solution derived by AnalyticLC.The vertical axis is the minimal closest approach in the system in units of mutual Hill radii.The color of the points represents the value of σ Nbody , with points at which σ Nbody > 1.5 (which do not pass our acceptance criterion) also marked with black edges.

Figure 3 .
Figure3.Orbital period, radii, and densities of the planets with mass estimates from this work.KOI numbers for which at least one planet does not have a previously reported mass value or upper limit are shown in underlined bold text.The size of the circles represents the absolute planetary size, and the color scale indicates our adopted density.Only planets with densities estimated to the significance of more than 4σ are color filled; others are open.For some of the planets (e.g., in the systems KOI-593 and KOI-1589) the relative masses are well constrained, but the absolute mass (and hence density) is not well constrained owing to a high uncertainty on the stellar mass and radius, and hence these planets are shown here open.For reference, the legend shows the size and density of Earth, Jupiter, and Neptune.The short vertical black lines indicate the locations of the closest first-order MMRs to the observed period ratio of adjacent planets.The black plus signs similarly indicate the locations of second-order MMRs-note that many planets are close to these MMRs.Both firstand second-order MMRs are indicated only if they are close to the observed period ratios.We do not show second-order MMRs that are a multiplication of a first-order MMR (e.g., 4:2 and 2:1).

Figure 4 .
Figure 4. Planetary masses and radii obtained from this work.This figure shows the overall spread of masses and radii and the good agreement of planet masses with former literature values.(a) Mass-radius diagram.Each blue error bar is related to a single planet.The red lines are normalized histograms of the masses and radii obtained by summing up the pdf's of all points.As the sample is small and radii are well determined, the radius-weighted histogram appeared more as a collection of discrete values, so it was smoothed using a Gaussian kernel of width 0.3 R ⊕ to better show population-wide trends.The gray contours are constant-density curves.(b) Comparison of masses obtained in this work with literature values: JH21 values (blue) and upper limits (red); HL17 default mass prior, three of which show upper limits only (yellow); and high-mass prior (purple).We omit planets for which literature values are within less than 2σ of zero.For six objects labeled with their KOI numbers, the results obtained in this work disagree with the results of HL17 by more than 2σ; these are discussed in more detail in the text.The black line shows the identity function.

Figure 5 .
Figure 5. Overall statistics of mass estimates in terms of fractional mass error and TTV magnitude for the planets in this study and our former study (Paper II).(a) Distribution of relative error in planetary mass for the combined results of this study and our former study (Paper II; blue) and past literature (orange).The literature sources taken into account here are detailed in Figure 10 and in Judkovsky et al. (2022b, Figure 10).(b) Distribution of σ TTV , the standard deviation of the TTV signal, in the population of planets for which we obtained a dynamical solution and provided masses in this work and in Paper II.

Figure 6 .
Figure 6.An example for a nonlinear TbV pattern, which is superimposed on a secularly driven TbV linear pattern.(a) The evolution of the impact parameters of the planets in the system KOI-935 (Kepler-31) for the sample median parameter values, shifted by the initial b value, b initial , to make the trends visible at the same scale.(b) b resid , the residuals of this evolution after the removal of the linear trend (which is attributed to secular interactions).The thick lines show the residuals for the sample median parameter values, and the narrow lines represent the dynamics for 10 randomly selected parameter sets from our sample to highlight the nonlinear TbV pattern in this system seen in all solutions.The outermost pair, with orbital periods of roughly 42.63 and 87.65 days, reside near but just outside the 2:1 MMR, causing a super-period of roughly 1570 days, which is clearly seen in the nonlinear TbV signal.The largest nonlinear part of the TbV has a semiamplitude of 4 • 10 −3 in the best-fitting solution.Such a small signal by itself might not have been detected by direct methods that aim to directly measure the impact parameter for each individual event; it was found here by integrating the information of all transit variations together (timing, duration, depth) in a single dynamical model.(c) The distribution of bresid s in the best-fitting parameters in the adopted dynamical solutions for 241 planets included in this work and in Paper II, and the linearly accumulated TbV over 1 yr for the same sample (using the median value of b  for each planet).The dashed lines indicate the median values of the two distributions; these differ by a factor of 10.

Figure 7 .
Figure 7. Distribution of TbV magnitude in our 130-planet sample for which b  is predicted with a significance level larger than 2σ, listed in the Appendix and in Paper II.(a) The absolute value of TbV rate against the shortest-period companion for the 130 planets in our sample that display TbV to better than 2σ.(b) The empirical cumulative distribution function of the TbV magnitude in this sample, showing that the typical value is of order 10 −2 yr −1 , with red circles indicating the median and the percentiles equivalent to 1σ and 2σ in a normal distribution.

Figure 8 .
Figure 8. Measured impact parameter variation rates.(a) Timescale for transits to disappear (namely, number of orbital periods for b to change by 1) against an orbital period of 130 TbV-displaying planets.The radii of the circles represent the absolute planetary radii.(b) Absolute values of b  vs. orbital period.The values of b  are spread around ≈10 −2 across the range of observed orbital periods.The correlation of b  | | with P is weak.

Figure 9 .
Figure9.Distribution of mutual inclinations among all pairs in all our runs that passed the dynamical validity criterion σ Nbody < 1.5, and after shuffling over all configurations of the signs of I x and I y .Each solid curve is obtained by averaging all empirical pdf's from all the runs passing the criterion above calculated over uniform bins in i m and shown in a logarithmic scale.The dashed lines indicate the best-fitting Gaussian for each curve: blue for systems with two or three planets, and red for systems with more than three planets.The vertical lines indicate the 99% confidence intervals of the expectation values of the fits (also listed in the legend), which are 2°.49 ± 0°.02 for systems with two or three planets and 2°± 0°.04 for systems with more than three planets.This gives a qualitatively similar picture to that ofHe et al. (2020) in both the distribution shape and its dependence on multiplicity.
degrees.The outermost transiting planet (KOI-152.04,Kepler-79 e) is almost grazing, with an (absolute value of the) impact parameter of roughly 0.96 and a planet-to-star radius ratio of about 0.026.The two intermediate planets (KOI-152.02and KOI-152.01,namely Kepler-79 c and Kepler-79 d) display significant TbVs of order 0.01 yr −1 , and the innermost planet (b) displays significant, somewhat slower TbVs of order 0.003 yr −1 .The best model TTV nicely matches the TTV data of Holczer et al. (2016).Our adopted solution suggests masses of m for KOI-152.04 (Kepler-79 e).These are within ∼1.5σ of the results of Jontof-Hutter et al. (

Figure 10 .
Figure10.Comparison of the masses obtained from this work with past literature by KOI.In each subpanel, the planets in the system are sorted (left to right) by their orbital period.The numbers above are the suffix of the planets' KOI designation; for instance, the first object in the plot is.Dotted vertical lines separate different planets within a system.Bold underlined KOI numbers are used for systems with at least one planet without a previously reported mass constraint, as in Figure3.The bold underlined KOI suffix represents a specific planet without a previously reported mass.The legend shows the colors and markers representing different literature sources.Arrows pointing down represent upper limits on mass.
for KOI-593.04; the three innermost transiting planets have a significant mass detection of more than 3σ.The solution suggests eccentricities of a few percent and inclinations of a few degrees.Kepler-33 (KOI-707).-Thisfive-planet system has been analyzed by HL14,Hadden & Lithwick (2016), and HL17.Our solution suggests masses of m 4 .03, where the innermost planet has a density slightly higher than Earth and the other planets have densities lower than Earth.The two outermost planets display significant TbVs; these may relate to the 5°-10°of mutual inclinations suggested by the solution; an alternative explanation could be an external nontransiting companion.Kepler-235 (KOI-812).-Wefound no literature mass for this four-planet system.Our solution suggests masses of m 5 . No significant TbVs are observed.Kepler-31 (KOI-935).-Thisfour-transiting-planet system has been analyzed byFabrycky et al. (2012), who used TTVs from the first eight Kepler quarters to conclude that the transit signals arising from KOI-935.01 (Kepler-31 b) and KOI-935.02(Kepler-31 c) arise from the same host star.HL17 estimated the masses of KOI-935.01 (Kepler-31 b), KOI-935.01 (Kepler-31 c), and KOI-935.01 (Kepler-31 d); KOI-935.04 is currently designated as a candidate.We provide three solutions and select the adopted one owing to fit quality; it is better than the others by Δχ 2 ∼ 230, equivalent to ∼7σ.Our adopted solutions suggest masses of m 01 (Kepler-338 b) and KOI-1930.04(Kepler-338 e) in HL14.Hence, our masses for KOI-1930.02(Kepler-338 c) and KOI-1930.03(Kepler-338 d) are the first estimates we are aware of.Our solution suggests masses of m 5 .04 (Kepler-341 e), as well as a coplanar structure with eccentricities of order a few percent.Kepler-85 (KOI-2038).-In the literature, we found a nominal mass estimate for planet Kepler-85 c by HL14 and mass estimates for all planets by HL17.We provide four different solutions, with the adopted one suggesting masses of m 4

Table 6
Instantaneous Coordinates, Velocities, and Orbital Elements of the Best-fitting Parameters for the Adopted Solutions in This Work and in Paper II