Edge-of-the-Multis: Evidence for a Transition in the Outer Architectures of Compact Multi-Planet Systems

Although the architectures of compact multiple-planet systems are well-characterized, there has been little examination of their"outer edges", or the locations of their outermost planets. Here we present evidence that the observed high-multiplicity Kepler systems truncate at smaller orbital periods than can be explained by geometric and detection biases alone. To show this, we considered the existence of hypothetical planets orbiting beyond the observed transiting planets with properties dictated by the"peas-in-a-pod"patterns of intra-system radius and period ratio uniformity. We evaluated the detectability of these hypothetical planets using (1) a novel approach for estimating the mutual inclination dispersion of multi-transiting systems based on transit chord length ratios and (2) a model of transit probability and detection efficiency that accounts for the impacts of planet multiplicity on completeness. Under the assumption that the"peas-in-a-pod"patterns continue to larger orbital separations than observed, we find that $\gtrsim35\%$ of Kepler compact multis should possess additional detected planets beyond the known planets, constituting a $\sim7\sigma$ discrepancy with the lack of such detections. These results indicate that the outer ($\sim100-300$ days) regions of compact multis experience a truncation (i.e. an"edge-of-the-multis") or a significant breakdown of the"peas-in-a-pod"patterns, in the form of systematically smaller radii or larger period ratios. We outline future observations that can distinguish these possibilities, and we discuss implications for planet formation theories.


INTRODUCTION
Short-period, sub-Neptune-sized planets are one of the most prevalent types of exoplanets in the local Galaxy. Discovered in abundance by NASA's Kepler mission , these planets are frequently found in multiple-planet systems with orbital periods ranging from days to months (Lissauer et al. 2011Rowe et al. 2014;Fabrycky et al. 2014). Over the last decade, the architectures of these compact-multiple planet systems ("compact multis") have been characterized in detail. The observed planets have low eccentricities (e.g. Van Eylen & Albrecht 2015;Xie et al. 2016), low inclinations (e.g. Fang & Margot 2012;Fabrycky et al. 2014), and tight orbital spacings that are generally not near sarah.millholland@mit.edu * NASA Sagan Fellow resonances (Lissauer et al. 2011;Fabrycky et al. 2014;Winn & Fabrycky 2015). Within a given system, the period ratios, planet radii, and planet masses tend to be significantly more uniform than would be expected by random chance -a set of patterns known collectively as the "peas-in-a-pod patterns" or "intra-system uniformity" (Weiss et al. 2018, Millholland et al. 2017; for a review, see Weiss et al. 2022).
The statistical properties of the observed compact multis have in turn revealed corresponding properties of the underlying distribution of planetary systems (i.e. including undetected planets) through the help of forward models. For instance, the fraction of stars with multiple planets interior to ∼ 1 AU is estimated at ∼ 30% − 70% Mulders et al. 2018;Zink et al. 2019;He et al. 2019), with the occurrence increasing for cooler stars (Howard et al. 2012;Mulders et al. 2015;Yang et al. 2020;He et al. 2021). There is evidence that the eccentricities and inclinations correlate with the intrinsic planet multiplicity He et al. 2020;Zhu & Dong 2021) and that the inclinations are both relatively small overall and not bifurcated into dynamically cool and dynamically hot sub-populations . Finally, forward models have shown that the patterns of intra-system uniformity in planetary sizes and orbital spacings must be present to a strong degree in the underlying population (He et al. 2019(He et al. , 2020Gilbert & Fabrycky 2020;Mishra et al. 2021), indicating that these patterns cannot be explained by detection biases and require astrophysical origins.
There is an additional feature of compact multis that, relative to the properties mentioned above, has received comparatively little attention. That is their inner and outer edges, which can be summarized statistically as the distributions of orbital periods of the innermost and outermost planets in the systems. The edges of compact multis are not well understood from either an observational or theoretical perspective. However, these features (specifically those of the underlying population of systems) are crucial for a complete characterization of short-period tightly-packed systems.
The inner edges are the simpler of the two from an observational perspective, since the geometric and detection biases are comparatively minimal. A simple calculation based on a sample of Kepler systems with four or more observed planets (to be defined later in this paper) shows that the distribution of observed innermost planet periods has a median of 3.9 days with the 16th-84th percentile interval equal to [2.2, 7.0] days. After accounting for biases, Mulders et al. (2018) showed that the underlying distribution of inner edges likely peaks at slightly longer orbital periods of about ∼ 10 days. This is consistent with being a signature from the protoplanetary disk inner edges (e.g. Millan-Gabet et al. 2007), which may limit where planets are able to form in situ (e.g. Lee & Chiang 2017) or may act as a planet trap that halts inward migration (e.g Terquem & Papaloizou 2007;Izidoro et al. 2017).
The outer edges are considerably more difficult to understand because they are heavily impacted by geometric and detection biases. Specifically, the transit probability and transit signal-to-noise ratio decrease with increasing orbital period as p trans ∼ P −2/3 and SNR ∼ P −1/3 , respectively. Moreover, transit surveys require a minimum number of transits for detection (three for the Kepler pipeline), which, together with the observation baseline and duty cycle, limits the largest detectable orbital period. This upper limit is P 500 days for the Kepler prime mission. For a sample of Kepler systems with four or more observed planets (the same sample as discussed in the previous paragraph), the distribution of outermost planet periods has a median of 40.6 days with the 16th-84th percentile interval equal to [16.3, 76.9] days. Naively speaking, this distribution already seems inconsistent with expectations, since almost all of the outermost periods are more than an order of magnitude smaller than the upper limit of ∼500 days. However, such comparisons are not very meaningful until the relevant geometric and detection biases have been thoroughly accounted for.
In this work, we examine the outer edges of compact multis discovered by the Kepler prime mission. 1 Our main goal is to understand whether the outer edges of the observed systems are consistent with sculpting purely from geometric and detection biases acting upon an underlying population of systems that extend out to (and potentially beyond) the largest periods probed by the Kepler photometry. Alternatively, if we find that the observed systems truncate at smaller orbital periods than required by their detection, then this could indicate either that (1) the underlying planets in compact multis are present only within some restricted range of orbital separations and/or (2) there are other significant changes in the planet properties (e.g. periods, radii) at larger orbital separations. To explore this, we will consider the existence of hypothetical planets orbiting beyond the outermost planets in Kepler multis, and we will estimate the number of these hypothetical planets that we would expect to be detectable.
It is important to note that this investigation of the outer edges of compact multis has consequences beyond demographics; specifically, the outer edges are signatures of planet formation and dynamical evolution. Just as the inner edges of compact multis are likely relics from the disk inner edges, a possible truncation in the outer systems could be the result of disk migration traps (e.g. Zawadzki et al. 2022), planet formation in pebble rings (e.g Chatterjee & Tan 2014), or dynamical perturbations from exterior giant planets (e.g. Pu & Lai 2018), among other possibilities. We will review such theoretical considerations in Section 7. However, it is helpful to keep this motivation in mind from the outset. This paper is organized as follows. We begin with a heuristic demonstration of the detectability of potential transiting outer planets, as a motivation for the more detailed calculations in subsequent sections (Section 2). We then investigate the geometric aspects of the prob-lem and develop a new method for estimating a system's mean inclination and mutual inclination dispersion (Section 3). Based on these estimates, we define a model for the transit probability, as well as a detection efficiency model (Section 4). We then validate these approaches on populations of simulated planets (Section 5). Moving onto observed systems, we estimate the number of hypothetical additional outer planets that we would expect to be transiting and detectable in Kepler multiple-planet systems (Section 6). Finally, we discuss implications of our results for planet formation theories (Section 7).

HEURISTIC CALCULATION
We begin with a heuristic exploration of the detectability of potential transiting outer planets in Kepler multi-planet systems. This offers a simplified preview of later detailed calculations. We consider a sample of 64 Kepler systems with four or more transiting planets. (The details of the sample selection will be provided in Section 6.1.) We imagine that each system has an additional transiting outer planet with properties consistent with the "peas-in-a-pod" patterns. Thus, we set the radii and periods of these "hypothetical planets" to R p,hypo = R p,out and P hypo /P out = P out /P 2nd−out , where P out and R p,out are the period and radius of the outermost observed planet, and P 2nd−out is the period of the second outermost observed planet.
Next, we compare R p,hypo to the minimum planet radius, R p,min , that the hypothetical planets would require in order to be detected with signal-to-noise ratio, SNR > SNR min = 7, which is approximately equal to the detection threshold of the Kepler pipeline. We define the SNR as (Christiansen et al. 2012) (1) Here, t obs is the time interval over which data were collected (∼4 years in most cases), f 0 is the duty cycle, and CDPP eff is the effective combined differential photometric precision, where T ∝ P −1/3 is the transit duration (see equation 4 for the full expression). In Section 4.2, we will use a more advanced definition of the SNR, but this definition suffices for the purposes of this illustration. Given that SNR ∝ P −1/3 R p 2 , we can calculate the minimum planet radius, R p,min , for which SNR > SNR min = 7 at P = P hypo using a scaling to SNR out , the SNR of the outermost observed planet in the system,  Figure 1 shows the architectures of Kepler systems with four or more planets, along with dots indicating the radii R p,hypo and R p,min at period P hypo . For the purposes of the visualization, we include only systems for which all planets in the system have R p < 8 R ⊕ , leaving out 5 systems. Remarkably, in all 59 systems, R p,hypo > R p,min . In other words, the hypothetical planet's nominal radius is larger than the minimum radius required for SNR > 7 in every single case. This result can be seen by the purple outlines around the yellow dots in Figure 1 and the R p,hypo /R p,min ratios on the right-hand-side. If we increase P hypo , then R p,hypo > R p,min remains to be generally true but not in every case. For instance, if we set P hypo = 300 days in all systems, then R p,hypo > R p,min in 45 out of 59 systems.
To summarize, this brief experiment has demonstrated the following: If additional transiting planets exist beyond the known transiting planets, and if they have properties in line with expectations from the "peas-in-apod" patterns, then we would generally expect them to be detectable. In the remainder of the paper, we explore this simple premise in greater detail. We consider similar experiments involving hypothetical outer planets, but we consider both transit probabilities and detection probabilities and formulate detailed expectations of the frequency of such planets we would expect to be transiting and detectable. In the next section, we begin with the key geometric aspect of the problem: mutual orbital inclinations. One of our subsequent objectives is to quantify the likelihood that planets orbiting beyond the known planets in a multi-transiting system are also transiting. In order to do this, we must first obtain estimates of the mean inclination and inclination dispersion. Here we develop a new procedure called the "transit chord ratio method", which uses the ratios of the transit chord lengths of pairs of planets in the same system (Steffen et al. 2010) as a means of constraining the planets' inclinations. This has been done at the population level (e.g. Fang & Margot 2012;Fabrycky et al. 2014) but not the level of individual systems. One benefit of working with ratios of transit chord lengths (rather than each planet's individual transit duration) is that it significantly reduces the impact of the degeneracy between the stellar density and impact parameter, which are difficult to accurately measure with Kepler long cadence photometry (Petigura 2020).
It is important to note that transit chord ratios are sensitive specifically to sky-plane inclinations, the inclinations between the planetary orbit planes and the sky plane. Accordingly, they cannot directly constrain the dispersion of the inclinations between the planetary orbits and the invariable plane (the plane perpendicular to the system's total angular momentum vector). Moreover, since the transit chord lengths are symmetric on the two stellar hemispheres bisected by the b = 0 chord, the method cannot determine whether two planets are transiting on the same hemisphere or opposite ones. Despite these limitations, the method is still capable of providing useful estimates of the mean inclination and inclination dispersion of multi-transiting systems, as we will soon demonstrate.
In the limit R p R a, a planet's transit duration is given by where R is the stellar radius, b is the dimensionless impact parameter, and v mid is the sky-projected orbital velocity at mid-transit. Here we will make the simplifying assumption of circular orbits, 2 such that v mid = na, where n = 2π/P is the mean-motion and a is the semimajor axis. The ratio of transit chord lengths (T v mid ) of a pair of planets (j and k with j < k) in the same system is where on the right-hand side of the first line, we have made the substitution bR = a cos i, with i the inclination between the planet's orbit plane and the sky plane.
In the final line, we expressed a/R in terms of the orbital period and stellar density, ρ . Since T and P are well-measured from transit data, the observed ratios (ξ j,k ) obs = (T k /T j )(P k /P j ) −1/3 can provide constraints on the inclinations (Fabrycky et al. 2014). The uncertainties on (ξ j,k ) obs can be calculated through standard error propagation using the reported uncertainties on T and P . Consider the ratio of transit chord lengths of each planet in a multi-planet system with respect to that of the innermost planet, ξ 1,k , where k = 1, 2, ..., N . For a perfectly coplanar system with constant sky-plane inclination i, data points of (P , ξ 1,k ) would trace out a smooth and monotonic curve. Examples of these curves are depicted in Figure 2 and are calculated by setting the inclinations to i in equation 5. For a set of observed (ξ 1,k ) obs values in a system, we can identify a model curve, (ξ 1,k ) mod , that best fits the observations. The spread of the (ξ 1,k ) obs measurements around the best-fit curve indicates the degree of sky-plane inclination dispersion about the mean sky-plane inclination.
When fitting the model curve, we allow (ξ 1,k ) mod to be modified by a vertical scaling factor, γ, which accounts for the fact that i 1 may not be near i. In other words, (ξ 1,k ) mod need not pass through unity at P = P 1 , since this can otherwise bias the fit. We use non-linear least squares to find the best-fit values of i and γ that minimize the sum of the squared residuals of (ξ 1,k ) mod − (ξ 1,k ) obs . We require i ≤ 90 • in the fit, which is not a loss of generality due to the symmetry of the transit chord lengths on the two stellar hemispheres bisected by the b = 0 chord. The best-fit values of i and γ uniquely yield i 1 , which we then use to solve for the other planets' inclinations from equation 5. Finally, we measure the sky-plane inclination dispersion as the root mean square of ∆i k = i k − i, We illustrate the transit chord ratio method using the TRAPPIST-1 system (Gillon et al. 2017) as an example. With well-determined parameters for all seven planets (Agol et al. 2021), the TRAPPIST-1 system allows us to compare the constraints on the mean inclination i and dispersion σ i from the transit chord ratio method to those obtained from independent techniques. We adopt the values of ρ and each planet's T and P from Agol et al. (2021). Figure 2 shows the observed ratios, (ξ 1,k ) obs = (T k /T 1 )(P k /P 1 ) −1/3 , along with the model curve with best-fit parameters, i = 89.781 • and γ = 0.993. The observed values are tightly distributed around the best-fit curve, yielding a dispersion equal to σ i = 0.105 • . The calculated estimates of i and σ i agree well with those derived from the photodynamical model in Agol et al. (2021), where the mean and dispersion of the i values are 89.783 • ± 0.032 • and 0.053 • ± 0.034 • .
The TRAPPIST-1 system is exceptionally close to coplanar. By experimenting with synthetic systems with larger inclination dispersion (such as some of the systems that will be described in the next section), we identified that when σ i 0.4 • , the fit sometimes forces i to Transit chord ratios for planets in the TRAPPIST-1 system. The black dots with error bars indicate the observed ratios, (ξ 1,k ) obs = (T k /T1)(P k /P1) −1/3 . The solid curves represent coplanar models, (ξ 1,k ) mod , with different inclinations, ranging from i = 89.5 • to i = 90 • in steps of 0.0625 • . The black dashed line indicates the model curve with best-fit parameters i = 89.781 • and γ = 0.993. nearly its maximum possible value, 90 • . This occurs when several planets have |b k | < |b 1 |, which leads to the corresponding ratios (ξ 1,k ) obs > 1. The model attempts to capture these points in the fit (along with those with (ξ 1,k ) obs < 1) by forcing i to 90 • , which corresponds to a horizontal line of (ξ 1,k ) mod vs. P . To remedy these cases, we utilize a second estimate of i that is derived by first estimating the inclination for each planet directly from its measured transit duration (equation 4) and then taking the mean over all planets, i circ . For cases where the transit chord ratio method yields σ i > 0.4 • and i > 89.99 • , we replace the original value of i with (i + i circ )/2. Whenever we replace i, we also recalculate σ i using equation 6.

Validation using simulated planetary systems
We can examine the accuracy of the transit chord ratio method using simulated planetary systems, for which we can compare the derived mean inclination and inclination dispersion to the system's true values. We consider simulated planetary systems from the SysSim (short for "Planetary Systems Simulator") forward modeling framework 3 (Hsu et al. 2018(Hsu et al. , 2019He et al. 2019He et al. , 2020He et al. , 2021 models, the parameters of which are determined using a calibration to summary statistics of the observed Kepler planet population (including the observed distributions of multiplicities, orbital periods and period ratios, transit depths and depth ratios, etc.).
It is important to distinguish between two different planetary populations that are generated within a Sys-Sim forward model. First, there is the "physical catalog", the underlying population of planetary systems, which is drawn directly from the statistical model. Second, there is the "observed catalog", which is obtained by passing the physical catalog through a simulated Kepler detection pipeline and determining which planets would be detected and labeled as planet candidates during the automated vetting process. Once the best-fit parameters of the forward model have been identified (see e.g. Hsu et al. 2019;He et al. 2019 for details), the observed catalog closely resembles the population of observed Kepler systems.
In this study, we will consider simulated planetary systems from the latest edition of SysSim: the "maximum AMD model" (He et al. 2020). This model is based on the argument that a system's long-term orbital stability is closely related to its angular momentum deficit (AMD), the difference between the total orbital angular momentum of the system and what it would be if all orbits were circular and coplanar (e.g. Laskar & Petit 2017). The key assumption of the model is that all systems have the critical AMD for stability. Although we do not expect this to be true for all systems in reality, the model offers an intuitive framework for assigning system orbital properties and has been shown to reproduce many aspects of the Kepler data (He et al. 2020), including Kepler planet Transit Duration Variations (TDVs), which probe mutual orbital inclinations . In this work, we will use pairs of physical and observed catalogs generated from the maximum AMD model. Each catalog pair corresponds to a different set of model parameters sampled from the posterior distributions derived by He et al. (2020).
We utilize ten physical/observed catalog pairs and consider only the systems with four or more detected planets ("4+ systems"), equaling 374 systems total. The observed catalogs contain estimates of T and P of the simulated observed planets with realistic measurement precision. They also contain the "true" values of the stellar density (i.e., assuming perfect measurement precision). In order to simulate realistic uncertainties of the stellar density measurements, we utilize the fractional uncertainties, f ρ = σ ρ /ρ , of Kepler host stars from the Gaia-Kepler Stellar Properties Catalog (Berger et al. 2020a), where we consider only the subset of host stars with four or more observed planets. The median f ρ of this sub-sample is 9%. For each ρ ,true value in the SysSim catalog, we draw a random fractional uncertainty, f ρ ,rand . We then draw a new value of the density as ρ ∼ N (ρ ,true , ρ ,true × f ρ ,rand ). It is important to note that the transit chord ratio method is not very sensitive to ρ because the process of using ratios of transit chord lengths allows this dependency to partially cancel out.
We use our method to estimate i and σ i in each Sys-Sim 4+ system using T , P , and ρ values as inputs. Figure 3 shows the comparison between the derived values and true values of i and σ i . There is good agreement overall, with relatively tight clustering of the calculated vs. true values around the one-to-one line. The median and 16th and 84th percentiles of the distribution of dif- In addition, the offset of the calculated σ i values from the true values is weakly correlated with the true i. This occurs because systems with i far from 90 • are more likely to have all planets transiting on one hemisphere of the star, which yields σ i values that are more accurate but biased high due to the effects of the unmodeled eccentricities. On the contrary, systems with i close to 90 • are more likely to have some planets transiting on separate hemispheres, which is not captured within the transit chord ratios and leads to underestimated σ i values. As for the i estimates, the median and 16th and 84th percentiles of the distribution of differences in i are 0.03 •+0.29 • −0.28 • . We also compare (but do not plot) the true values of σ i with σ i,inv , the dispersion of the inclinations with respect to the invariable plane. We find that median(σ i,inv /σ i ) ≈ 1.6.
Before proceeding to utilize i and σ i to calculate transit probabilities, we introduce an additional scaling factor to σ i . Recall that σ i is the sky-plane inclination dispersion of the observed transiting planets. By definition, σ i does not contain information about the skyplane inclinations of the non-transiting planets. In order to approximately rectify this, we consider each SysSim system with four or more detected planets, and we calculate the sky-plane inclination dispersion, σ i,all , of all planets in the system (i.e. including undetected planets) using the true values of the inclinations. We then compute the ratio of σ i,all and the calculated σ i of the detected planets (as defined in equation 6). We find that median(σ i,all /σ i ) = 1.2. We thus define a scaled sky-plane inclination dispersion, σ i = 1.2σ i , that now approximately reflects the spread of inclinations of the full system. We will use σ i whenever we're referring to the dispersion measured from the detected planets and σ i whenever we're working in the context of transit probabilities, as in the next section. Once the mean inclination, i, and scaled inclination dispersion, σ i , have been estimated for a given system, it is possible to calculate the probability that a hypothetical planet with an arbitrary period transits the star. This requires an assumption that the inclinations of individual planets in the system are well-described by some distribution with mean and standard deviation equal to i and σ i . While there are multiple distributions one could consider, we find that a normal distribution, The calculation of the transit probability first requires specifying the bounded range of inclinations for which a planet is transiting. Again assuming circular orbits, a planet transits with |b| < b max if its inclination falls between i min < i < π − i min , where The probability, p trans , that a planet with period P transits with |b| < b max is thus We illustrate the transit probability calculations using two examples of SysSim systems with four detected planets. Figure 4 shows the inclination range within which a planet would transit with |b| < b max = 1, as well as the transit probability assuming that i ∼ N i, σ i . The transit probability is unity out to a period large enough such that the transiting region crosses over the region of higher probability according to i ∼ N (i, σ i ). There is close agreement between the two transit probability curves associated with the calculated (solid purple) and true (dashed purple) values of i and σ i . As in these example SysSim systems, each system has its own unique curve of transit probability versus P given the set of parameters ρ , i, and σ i . One can thus evaluate the transit probability of observed or hypothetical planets at arbitrary periods.

Detection probabilities
Even if a planet is transiting, it is not guaranteed to be detected. Here we outline the calculation of the detection probability, which depends on the statistical significance of the transit signal. The Kepler transiting planet search (TPS) pipeline required a planet to be transiting at least three times and to have a statistical significance exceeding a threshold of 7.1σ in order to  be detected (Christiansen et al. 2020). The statistical significance of the transit signal is computed with the multiple event statistic (MES), which cannot be specified exactly without applying the TPS pipeline to a given lightcurve but can be approximately computed using the one-sigma depth function (OSDF). The OSDF is a data product provided by Kepler DR25 (Thompson et al. 2018) that quantifies the transit signal that would be expected to result in an MES equal to unity for a given target star, orbital period, and transit duration, after averaging over the epoch of transit (Burke & Catanzarite 2017). Thus, the expected MES is simply where δ is the transit depth. Kepler DR25 provided large tables of OSDF values for each target star as a function of 14 transit durations and ∼ 10 4 orbital periods; we use the downsampled versions of these tables provided by Hsu et al. (2019). Given a target star and a planet's P and T , we use bilinear interpolation to the star's OSDF table and extrapolate whenever P or T are beyond the range.
We can now define a detection efficiency model. We build upon the model from Hsu et al. (2019), which combines the probability of a transiting planet being detected and passing vetting, such that it is labeled a planet candidate. The model is calibrated to the Kepler DR25 pixel-level transit injection tests (Christiansen 2017) and corresponding robovetter results (Coughlin 2017). The probability of the planet passing through both detection and vetting is (Hsu et al. 2019) where N tr is the number of valid transits observed by Kepler, and α Ntr , β Ntr , and c Ntr are N tr -dependent parameters provided in Table 1. The number of transits can be approximated as N tr = floor(t obs f 0 /P ), where t obs and f 0 are the target-specific data span and duty cycle. We apply additional modifications to our model to account for biases that can occur specifically in multiplanet systems, as shown by Zink et al. (2019). For example, the Kepler TPS pipeline masks each existing transit detection before searching for the next highest SNR candidate, reducing the available photometry for lower SNR candidate searches. These completeness issues are relevant since the theoretical planets considered in this study have the longest periods in their systems, with the fewest number of available transits. To ensure our estimates are conservative, we consider the case where these additional biases impact all the outermost planets.   Table 1. We acknowledge that a more accurate measure of multiplicity completeness would involve an injection/recovery test of multi-planet systems through the entirety of Kepler TPS pipeline, a procedure which is beyond the scope of this paper. Our estimates provide a first-order approximation, and the results in the following sections exceed any expected changes from a more refined completeness assessment. Before returning the final probability of detecting a planet, the probability in equation 10 is multiplied by an additional factor called the window function (WF), or the fraction of unique transit ephemeris epochs that permit three or more transits to be observed as a function of orbital period (Burke & Catanzarite 2017). This factor is necessary and complex because it provides the probability that a planet has N tr ≥ 3 given existing data gaps within the limited span of Kepler photometry, whereas equation 10 assumes N tr ≥ 3 to be the case. Similar to the OSDFs, Kepler DR25 provided tabulations of WFs for each star as a function of T and P . We use the downsampled versions provided by Hsu et al. (2019). We use bilinear interpolation to calculate the target star's WF as a function of P and T , extrapolating whenever these values are beyond the range.

VALIDATION TESTS WITH SIMULATED PLANETARY SYSTEMS
The previous section outlined the ingredients necessary to calculate an expected distribution of detections of outer planets orbiting beyond observed planets in multi-planet systems. Before applying this to the Kepler compact multis (Section 6), we will first test our methods on sub-populations of simulated planets from the SysSim maximum AMD model (He et al. 2020). We will essentially "re-simulate" the transit and detection observations using the methods developed in Section 4 as a means of validating that our approach results in accurate total numbers of transiting and detected planets. We will specifically examine outer planets in SysSim systems with four or more detected planets ("4+ systems") and verify that our calculations are in line with the true properties of the synthetic systems.
We begin with the same ten pairs of SysSim physical catalogs and observed catalogs as used in Section 3.2. Recall that a single physical/observed catalog pair represents an entire Kepler population of planetary systems orbiting FGK dwarfs. We consider the 4+ systems only.
In each system, we identify the innermost undetected planet beyond the detected planets, which we will call the "next outer undetected planet". For instance, in each example system in Figure 4, the fifth planet from left is the next outer undetected planet. Not all systems have such a planet; in ∼ 25% of systems, the outermost planet is a detected planet, so the next outer undetected planet is undefined. We focus on the remaining ∼ 75% of systems for which this planet is defined. The true number of detections in this sub-sample is zero by definition, so when we calculate the expected number of detections, we should get a result that is consistent with zero.
For each physical/observed catalog pair, we calculate the distribution of the number of next outer undetected planets that we would expect to be transiting and detectable. We use the following steps: 1. For each next outer undetected planet, we estimate the transit probability, p trans , according to equation 8, where i and σ i are either the true values or the calculated values from the transit chord ratio method.
2. We approximate each planet's transit duration and transit depth (if it were to transit) using scalings to the duration and depth of the outermost observed planet, T next = T out (P next /P out ) 1/3 δ next = δ out (R p,next /R p,out ) 2 .
3. Using the period, P next , transit duration, T next , and transit depth, δ next , we estimate the OSDF, WF, and MES for each planet using the targetspecific data tables and the procedure described in Section 4.2.
4. We use the MES and WF to estimate the detection probability, p det&vet , using equation 10.
5. Using the calculated probabilities, we draw two Bernoulli random variables, with the X trans indicating whether the planet is transiting and X detect indicating whether the planet is both transiting and detected.
6. The sum of Bernoulli random variables across all the systems in a given catalog pair gives a number of "successes" (transits or detections), 7. Finally, by repeating steps five and six 1000 times, we create distributions of N trans and N detect for each catalog pair (with 1000 values in each distribution). Figure 5 shows the distributions of N trans and N detect of the next outer undetected planets across ten SysSim catalog pairs. The first thing to note is that the distributions of N trans and N detect , which use calculated values of i and σ i for the p trans calculation (blue histograms), agree well with the corresponding distributions that use the true values of i and σ i (yellow histograms). Thus, we can trust that using calculated values of these key parameters will provide accurate distributions when we proceed to working with the Kepler multis.
A second observation from Figure 5 is that the distributions of N trans and N detect are in reasonable agreement with the true values of N trans and N detect . (Recall that the true value of N detect is always zero because these experiments are conditioned upon predictions for the first outer undetected planet in each system.) The distributions are consistent with the true values to within ∼ 1 − 3σ, with the average being ∼ 2.2σ for N trans and ∼ 2.4σ for N detect . However, the distributions are always biased high. This is by construction of our exploration. Because the next outer undetected planets were missed within the actual SysSim simulations, they are more likely to be non-transiting than other planets in the same system. Thus, our assumption of i ∼ N (i, σ i ) in p trans (equation 8) partially breaks down and overestimates the transit probability. Regardless, the fact that the distributions agree with the true N trans and N detect to within 3σ ensures that we can reliably use these methods on the observed Kepler multis, as long as any results we find are significant by much more than 3σ.

DETECTION EXPECTATIONS FOR HYPOTHETICAL OUTER PLANETS IN KEPLER MULTI-PLANET SYSTEMS
In the previous section, we used populations of simulated planetary systems to validate our calculated expectations of the number of planets that would be transiting and detected. We now extend these methods to Kepler high-multiplicity systems. We aim to answer the questions: If the Kepler multis hosted additional outer planets with similar properties as the known planets, how many of them would we expect to be transiting and detectable? Moreover, are these expectations reconcilable with the lack of detections beyond the outermost known transiting planets?

Sample of Kepler multi-planet systems
We begin by defining our sample of Kepler multiplanet systems. We use the Kepler DR25 KOI catalog (Thompson et al. 2018; NASA Exoplanet Archive 2022) as our starting point, using all planets with "confirmed" and "candidate" dispositions. Where possible, we replace the stellar parameters and planet radii in the DR25 catalog with parameters from the Gaia-Kepler Stellar Properties Catalog (Berger et al. 2020a,b). In addition, we apply a small set of quality cuts. We consider only planets smaller than 16 R ⊕ with fractional radius uncertainties less than 100%. To avoid stars with large systematic radius errors, we discard targets for which Furlan et al. (2017) found a companion star that contributed more than 5% of the light in the photometric aperture. After these cuts, we are left with 64 Kepler systems with four or more observed transiting planets ("4+ systems").
We now apply the transit chord ratio method (Section 3) to the transiting planets in the Kepler 4+ systems to estimate the mean sky-plane inclination, i, and the skyplane inclination dispersion, σ i (from which the scaled dispersion σ i = 1.2σ i is subsequently derived). Figure  6 shows the distributions of i and σ i for the Kepler 4+ systems and the SysSim 4+ systems examined in Sec- SysSim systems. This offers another validation of the i and σ i estimates, which will be used in transit probability calculations in the next section.

Set-up of hypothetical planet experiments
We now construct several experiments consisting of "hypothetical planets" orbiting beyond the observed planets in Kepler 4+ systems. In general, the hypothetical planets are designed to emulate the continuation of the "peas-in-a-pod" architectures (Weiss et al. 2018;Millholland et al. 2017) by positing the existence of one additional outer planet in each system. We first assign the periods, P hypo , and radii, R p,hypo , of the hypothetical planets. We consider three separate approaches, outlined below and denoted "fixed", "random", and "targeted" sampling.
Fixed sampling: The most straightforward approach in assigning P hypo and R p,hypo is to consider that they perfectly follow the "peas-in-a-pod" patterns (Weiss et al. 2018;Millholland et al. 2017), or the statistical tendency for systems to exhibit intra-system uniformity in planetary radii and period ratios. Accordingly, we Figure 7. Distributions of the number of "hypothetical planets" that we would expect to be transiting, Ntrans, and detected, N detect , if they existed at the outer edges of the 64 Kepler systems with four or more known transiting planets. The left column corresponds to the fixed sampling approach with R p,hypo = Rp,out and P hypo /Pout = Pout/P 2nd−out . The right column corresponds to the random sampling approach with R p,hypo /Rp,out and (P hypo /Pout)/(Pout/P 2nd−out ) sampled from empirical distributions. The dashed histograms in the top panels represent Ntrans, while the solid histograms in the bottom panels represent N detect . The different colors represent ten separate iterations and have no significance otherwise.
can assign R p,hypo = R p,out P hypo /P out = P out /P 2nd−out , where P out and R p,out are the period and radius of the outermost known planet, and P 2nd−out is the period of the second outermost known planet. (This is identical to the initial experiments shown in Section 2.) Random sampling: The patterns of intra-system uniformity of radii and period ratios have an inherent degree of scatter (e.g. Millholland & Winn 2021). Thus, a second and more physical approach to assigning P hypo and R p,hypo is to utilize empirical distributions of radius ratios and ratios of period ratios from our sample. We assign where R p,i+1 /R p,i is randomly drawn from the observed distribution of radius ratios of adjacent planets, and (P i+1 /P i )/(P i /P i−1 ) is drawn from the observed distribution of ratios of period ratios of adjacent planet pairs. However, in an effort to preserve covariances of the dis-tributions with respect to other variables (e.g. correlations between R p,i+1 /R p,i and P i+1 /P i ), our random draws are not taken from the distributions as a whole but rather from pre-defined subsets. We bin the distribution of R p,i+1 /R p,i in five bins of P i and four bins of P i+1 /P i with an average of about 40 points in each bin. The bins are chosen so as to capture the large-scale correlations within the distribution, and the overall results are not sensitive to them. When sampling from the distribution of R p,i+1 /R p,i , we identify the appropriate bin and only sample from the corresponding subset of the distribution. Similarly, we bin the distribution of (P i+1 /P i )/(P i /P i−1 ) in six bins of the inner period ratio, P i /P i−1 , and we sample from the distribution by first identifying the appropriate bin and only drawing from the corresponding subset of the distribution. Targeted sampling: Finally, a third approach of assigning P hypo and R p,hypo is designed to allow us to systematically vary these parameters and observe the resulting changes to the distributions of N trans and N detect . Specifically, we sample R p,hypo /R p,out ∈ [0.2, 1] and (P hypo /P out )/(P out /P 2nd−out ) ∈ [1, 4], thus probing configurations in which the hypothetical planet is either smaller or further away than the "peas-in-a-pod" expectation. We target P hypo and R p,hypo one at a time using this sampling method, with the other parameter assigned using the "fixed sampling" approach.
In addition to the targeted sampling of P hypo and R p,hypo , we also extend our targeted sampling experiments to two other parameters. First, we consider assigning hypothetical planets to only a fraction of systems, rather than all systems. We sample this fraction in the range f hypo ∈ [0.1, 1]. This allows us to investigate how many systems can have additional outer planets and still be consistent with the observations. Lastly, in order to explore the effects of potentially biased estimations of the sky-plane inclination dispersion, we consider a constant scaling factor, γ σ i , such that σ i,new = γ σ i σ i . We sample the scaling factor in the range γ σ i ∈ [1, 10]. We reiterate that we use targeted sampling on only one of the four parameters (P hypo , R p,hypo , f hypo , and γ σ i ) at a time, with all others assigned using the "fixed sampling" approach (which, in the case of the latter two parameters, means f hypo = 1 and γ σ i = 1).
With the sampling approach defined and P hypo and R p,hypo assigned for each hypothetical planet, we proceed to calculate the distribution of the number of these planets that we would expect to be transiting and detectable. We use the same seven steps as outlined in Section 5, where the transit probabilities and detection probabilities are calculated for each hypothetical planet, and repeated trials of summed Bernoulli random variables yield distributions of N trans and N detect . 6.3. Results from fixed sampling and random sampling Figure 7 shows the results of the fixed sampling and random sampling approaches. Ten iterations of the distributions are shown, indicating the degree of variability resulting from the sampling of P hypo and R p,hypo combined with the seven steps outlined in Section 5. Immediately, we can see that the distributions of N detect are inconsistent with zero detections in both the fixed sampling and random sampling approaches. With the fixed sampling approach, the mean and standard deviation of the composite distribution of N detect (consisting of the ten separate distributions pooled together) is 27.55 ± 3.39 detected planets, or ∼ 8.1σ greater than zero detections. This indicates that ∼ 44% of the 64 Kepler 4+ systems would be expected to show detections of the hypothetical outer planets, if they existed.
As for the distributions resulting from the random sampling approach, there is more variability between different iterations due to each iteration having a distinct set of P hypo and R p,hypo values. The median and standard deviation of the composite distribution of N detect is 22.11 ± 3.25 detected planets, or ∼ 6.8σ greater than zero detections. This indicates that the hypothetical planets would be detectable in ∼ 35% of the 64 Kepler 4+ systems.
Both the fixed sampling and random sampling approaches indicate the same general result: If we posit the existence of additional planets orbiting beyond the known transiting planets with properties dictated by the expected "peas-in-a-pod" architectures, then these planets would be detectable in roughly ∼ 35% − 45% of systems, yielding a ∼ 7 − 8σ discrepancy with the lack of detections of such planets in the observed systems. However, it is important to note that the distributions of N detect are always positive by definition, so the more meaningful comparison is that with the corresponding results for the SysSim simulated systems (Section 5; Figure 5), where we found that the N detect distributions were biased high but always consistent with zero to within 3σ.
A final point is that this result is a lower limit on the expected number of detections since we only assigned one hypothetical planet per system. If we considered more than one outer planet per system, this would potentially generate an even larger number of expected detections. For instance, one possible experiment would be to add multiple sequential hypothetical planets to each system, assigned as P hypo,i = P out (P out /P 2nd−out ) i and restricted to P hypo,i < 500 days. We find that this would involve a total of 327 hypothetical planets, compared to the 64 in this study.

Results from targeted sampling
The discrepancy between the expected number of detections and the lack of such detections in the real systems indicates that one or more of our assumptions are breaking down. These assumptions and their related implications are enumerated below: 1. We assumed that each system has a single additional planet orbiting beyond the known transiting planets. However, perhaps not all systems have a such a planet. This would imply that the architectures of compact multis are truncated at a detectable orbital period in at least some fraction of systems.
2. We assumed that P hypo /P out ≈ P out /P 2nd−out . However, perhaps the period ratios are on average larger at sufficiently large orbital period. 5 Figure 8. Distributions of the number of "hypothetical planets" that we would expect to be transiting, Ntrans, and detected, N detect , if they existed at the outer edges of the 64 Kepler systems with four or more known transiting planets. All columns indicate the results of the targeted sampling approach, with the columns corresponding to (from left to right) sampling of f hypo (the fraction of systems with hypothetical planets), (P hypo /Pout)/(Pout/P 2nd−out ), R p,hypo /Rp,out, and γ σ i (the scaling factor of the inclination dispersion). These parameters are sampled from 10 evenly-spaced values within the ranges indicated at the top of the columns. The dashed histograms in the top panels represent Ntrans, and the solid histograms in the middle panels represent N detect . The colors indicate the value of the targeted parameter, as shown in the bottom panels, where the 1σ and 3σ intervals around the medians of the distributions are indicated with the darker and lighter errorbars. The two curves in the bottom panels are artificially offset vertically such that the errorbars can be seen clearly.
3. We also assumed that R p,hypo ≈ R p,out . However, perhaps the radii of planets are on average smaller at sufficiently large orbital period.
4. Finally, we also assumed a mutual inclination dispersion for each system based on our estimation from the transit chord ratio method. Although we validated the method in Section 3, it is possible that the σ i estimates are systematically low or that the true dispersion increases at large periods.
The targeted sampling approach described earlier in this section allows us to investigate these four possibilities by systematically varying f hypo , P hypo , R p,hypo , and γ σ i , which directly map to the four possibilities outlined above. All other parameters are held fixed while the targeted parameters are varied. Figure 8 shows the results of the targeted sampling. The bottom panels summarize the distributions by showing the medians and 1σ and 3σ intervals as a function of the targeted parame-ter. This allows us to visualize how much the targeted parameter must be varied in order to have the N detect distribution be approximately consistent with zero detections to within 3σ.
As shown in the figure, the N detect distribution is consistent (or, in some cases, nearly consistent) with zero to within 3σ when one or more of the following conditions are met: f hypo 0.2, (P hypo /P out )/(P out /P 2nd−out ) 3, R p,hypo /R p,out 0.5, or γ σ i 6. These thresholds are both approximate and conservative. The conditions on P hypo and R p,hypo indicate that the systems would require a significant deviation from the typical "peas-ina-pod" architectures. Moreover, the condition on f hypo indicates that only ∼ 20% of systems can host an additional outer planet with the expected architecture if the population is to be consistent with the observations. However, we note that the four conditions mentioned above could be relaxed if multiple parameter variations were considered simultaneously (e.g. both larger period ratios and smaller radii).
The condition on the σ i scaling factor, γ σ i , indicates that our inclination dispersion estimates would have to be very off (by more than a factor of ∼ 6) before the distribution of N detect is consistent with zero detections at 3σ. However, this raises questions about how large σ i can actually be before it becomes too unlikely that the observed planets are still co-transiting. In other words, perhaps we can rule out σ i values that are as large as those required to significantly shift the N detect distribution. We investigate this using synthetic trials in which we systematically increase γ σ i and determine how frequently the number of transiting planets equals the observed number. For 15 equally-spaced values of γ σ i , we consider each Kepler 4+ system, and we reassign each observed planet in the system a random skyplane inclination according i ∼ N (i, γ σ i σ i ). We then determine how many planets are transiting, N sim , and compare to the observed number of transiting planets, N obs . We calculate the fraction of all 64 systems for which N sim = N obs . Finally, we repeat these random trials 100 times. Figure 9 shows the results of these simulations. The fraction of Kepler 4+ systems for which N sim = N obs rapidly decreases as a function of γ σ i . Note that the fraction is not unity for γ σ i = 1 because the process of resampling the inclinations makes it unlikely to recover all of the exact same planets in transit, given that the observed sample is conditioned upon the planets transiting. With γ σ i as large as ∼ 6 (the value required earlier), less than 15% of systems would still have as many transiting planets as the observed systems. We can thus safely rule out the possibility that we have underestimated σ i by a factor as large as required to make the N detect distribution consistent with zero. However, this still leaves open the possibility that the edge-of-themultis discrepancy is explained by a substantial increase in the mutual inclinations at long periods (e.g. Shahaf et al. 2021).

Constraints on the truncation location
The hypothetical planet experiments from the previous section indicate that the outer edges of the observed Kepler 4+ systems likely cannot be sculpted by geometric and detection biases alone. Rather, these systems provide evidence either for an average truncation (i.e. occurrence rate drop-off) in the underlying architectures and/or some breakdown of the "peas-in-a-pod" patterns at larger orbital separations. We have not yet discussed where this transition occurs, since our hypo- i , i n v = 1 i , i n v = 2 Figure 9. Results of simulations that examine the impact of varying γ σ i , the scaling factor of σ i , on the transit multiplicity of the observed planets in the Kepler 4+ systems. The y-axis shows the fraction of the 64 Kepler 4+ systems for which the simulated number of transiting planets, Nsim, is equal to the observed number of planets in the system, N obs . The gray lines show the individual results from 100 simulations for each value of γ σ i . The yellow (purple) lines show the median and the 1σ (3σ) intervals. Benchmark values of the dispersion of the inclinations with respect to the invariable plane, σi,inv ≈ 1.6σi (Section 3.2), are also labeled.
thetical planet experiments were primarily designed to assess whether it exists rather than where it is. Even so, these experiments can still provide us with some insight. Figure 10 shows the distributions of P out , P hypo assuming the default "peas-in-a-pod" expectation (P hypo /P out = P out /P 2nd−out ), and P hypo assuming larger period ratios (P hypo /P out = 3P out /P 2nd−out ), which corresponds to the approximate condition for which the N detect distribution would be consistent with observations according to our targeted sampling experiments (Section 6.4). The medians of the three distributions are 40.6 days, 78.9 days, and 236.6 days, respectively. The distribution of P out shows a steep drop-off at around ∼ 100 days. Meanwhile, the third distribution rises beyond ∼ 100 days and peaks around ∼ 300 days. Since the third distribution corresponds to a set of hypothetical planets that would be consistent with zero detections to within 3σ, this indicates that the average edge must occur somewhere within the range P edge ≈ 100 − 300 days, or a edge ≈ 0.5 − 1 AU.

Theoretical interpretations
Whether our results are consistent with an average truncation in the underlying architectures or some breakdown of the "peas-in-a-pod" patterns, these findings have important implications for theories of the formation and dynamical evolution of compact multis. The theoretical interpretations will differ depending on which effect is dominant, but the same physical mechanism could, in principle, lead to both an occurrence rate decrease and a change in the architectures at larger separations. Nevertheless, we discuss outer truncation and a "peas-in-a-pod" breakdown separately below.

Outer truncation
With regard to outer truncation, one possibility is that, if compact multis form via orbital migration, they may often experience migration traps that prevent inward migration of some planets and cause large gaps in the systems at ∼ 100 − 300 days (Zawadzki et al. 2022; see also e.g. Paardekooper et al. 2010, Coleman & Nelson 2014, Bitsch et al. 2015, Izidoro et al. 2017, Carrera et al. 2019. Recently, Zawadzki et al. (2022) outlined that migration traps should occur in a timeevolving range of planet mass and semi-major axis (e.g. M p ∼ 0.5 − 2 M ⊕ , a 2 AU at ∼ 5 Myr) where the corotation torque dominates over the Lindblad torque and planets thus migrate outwards instead of inwards. Planets in this mass and semi-major axis range are prevented from reaching the inner regions of the disk, whereas more massive planets avoid the traps and end up on close-in and dynamically-cold orbits. Although Zawadzki et al. (2022) studied this theoretical framework in the context of possible explanations for the "Kepler dichotomy" (e.g. Lissauer et al. 2011;Johansen et al. 2012), it is relevant that their simulations naturally yield planetary systems that are bifurcated into two clusters with a large gap at ∼ 100 − 300 days, since this is precisely consistent with our observational findings.
Another physical mechanism that may cause outer architecture truncation is a proposed channel of compact multi formation called "inside-out planet formation" (Chatterjee & Tan 2014. In this framework, planets are formed sequentially from successive gravitationally unstable rings of ∼cm-m sized pebbles that drift inwards via gas drag and build up at the pressure maximum associated with the dead-zone boundary, which separates the inner and outer regions of the disk where the magnetorotational instability is and is not active. 6 Once a planet forms out of the pebble ring, it may migrate or isolate itself from the accretion flow, after which the dead-zone boundary retreats and the process repeats. The process is limited by the finite extent of the dead-zone boundary's retreat and by the mass reservoir that remains available for planet formation after the initial onset (Hu et al. 2018). Accordingly, this model predicts that planet formation will halt at some orbital separation, creating a break in the outer regions of compact multis.
In contrast to disk effects at the formation epoch, the outer edges of compact multis may also be limited by exterior perturbing planets. Recent results indicate that distant giant planets (with a 1 AU, M p 0.5 M Jup ) are common in systems with inner super-Earths Bryan et al. 2019). These distant giant planets are often eccentric and may be mutually inclined with respect to the inner planets (Dawson & Murray-Clay 2013;Masuda et al. 2020). Several authors have shown that distant giant planets can dynamically excite an inner system of super-Earths, provided that the distant perturber exerts a stronger gravitational influence on the inner planets than their mutual gravitational coupling (e.g. Lai & Pu 2017;Pu & Lai 2018;Denham et al. 2019;Spalding & Millholland 2020;Tamayo et al. 2021). The strength of the distant giant planet's perturbation relative to the inner planets' coupling depends on the separation between the outer perturber and the inner system (among other parameters). Accordingly, planets in compact multis that extend out to sufficiently large orbital separations may be dynamically excited or even destabilized by a distant giant. However, since distant giants are present in only a fraction of compact multi systems, their presence likely cannot be the sole explanation of our findings.

Breakdown of the "peas-in-a-pod" patterns
A different set of theories may be relevant in the case where our results are explained by smaller and/or more widely spaced planets at larger orbital separations. One possibility is that the dominant channel of compact multi formation -whether it is in situ accretion (e.g. Hansen & Murray 2012Chiang & Laughlin 2013;Dawson et al. 2015), distant assembly followed by inward migration (e.g. Terquem & Papaloizou 2007;Cossou et al. 2014;Coleman & Nelson 2014;Izidoro et al. 2017;Carrera et al. 2019), or some variation thereof (e.g. Chatterjee & Tan 2014) -might naturally produce smaller planet masses at larger orbital separations ( 0.5 − 1 AU) due to limitations of solid material or growth timescales.
For instance, in the planet formation simulations by Izidoro et al. (2017), in which sub-Neptune systems are produced through inward migration and disruption of resonant chains, the resulting systems appear to have a gap at around ∼ 0.5 − 1 AU, with the planets beyond the gap being less massive ( 1 M ⊕ ). This is also seen in other works such as Zawadzki et al. (2022) and is related to the slower migration timescales for the smaller mass planets. These small planets can also encounter migration traps where the co-rotation torque is dominant (as discussed in Section 7.2.1), whereas the larger mass planets avoid the traps and migrate inwards. This process naturally leads to a bifurcation of the system, with larger sub-Neptunes inside of ∼ 0.5 − 1 AU and smaller planets beyond that.
Stepping back from the specific details of planet formation physics, one can consider the general process of energy optimization at the assembly epoch, which also leads to the prediction that the "peas-in-a-pod" patterns should break down at larger orbital separations. Adams et al. (2020) (see also Adams 2019) calculated the lowest energy states available to forming planetary pairs subject to conservation of angular momentum, constant total mass, and fixed orbital spacing. They found that a configuration with approximately equal masses (i.e. "peas-in-a-pod") is the most energetically favorable when the total mass in the planets is less than some critical threshold, m c ∼ 40 M ⊕ . Above this threshold, the optimum state is one in which most of the mass is in one planet. The critical threshold m c decreases with increasing a, which indicates that mass uniformity is no longer energetically favorable at sufficiently large orbital distances ( 0.5 − 1 AU). Thus, the breakdown of the "peas-in-pod" patterns might simply be a consequence of planets optimizing the available energy.
To summarize, a variety of physical processes may be responsible for an outer architecture truncation or transition. It is possible that the same physical mechanism (e.g. migration traps) can cause multiple effects (e.g. occurrence rate decrease and smaller planets beyond ∼ 0.5 − 1 AU). We note that the list of theories above is by no means exhaustive, and it is beyond the scope of this paper to investigate all possibilities in detail. This topic is ripe for future theoretical investigations.

CONCLUSION
The inner and outer edges of compact multi-planet systems are fundamental signatures of their formation and evolution. The distribution of the innermost planet orbital periods of the underlying distribution of compact multis (i.e. after accounting for observational biases) peaks around ∼ 10 days (Mulders et al. 2018), consistent with being a relic from the protoplanetary disk inner edges (e.g. Terquem & Papaloizou 2007;Lee & Chiang 2017). The outer edges are significantly more affected by observational biases and thus have received no prior in-depth investigations, as far as we're aware. However, these biases are well-understood, meaning that a robust characterization is possible.
In this paper, we presented evidence that the outermost planets in the observed Kepler high-multiplicity systems truncate at smaller orbital periods than expected from geometric and detection biases alone. We showed this using experiments of the detectability of "hypothetical planets" orbiting beyond the outermost observed planets in Kepler high-multiplicity systems. These experiments were first demonstrated heuristically in Section 2 and then developed robustly in the remainder of the paper.
We considered one hypothetical outer planet per system for 64 Kepler systems with four or more observed planets, with each hypothetical planet's properties dictated by expectations from the "peas-in-a-pod" patterns (Weiss et al. 2018;Millholland et al. 2017). Using models of the transit and detection probabilities, we estimated that 22.11 ± 3.25 (approximately ∼ 35%) of these 64 planets would be transiting and detectable, constituting a ∼ 7σ difference with zero detections. This is significantly different than analogous results using simulated planetary systems, which were consistent with zero to within ∼ 2 − 3σ. We thus identified a strong discrepancy between expectations based on hypothetical planets and the lack of additional outer planet detections in the observed systems. Crucially, these results are a lower limit on the true extent of the discrepancy, since