An upper limit for the growth of inner planets?

The exotic range of known planetary systems has provoked an equally exotic range of physical explanations for their diverse architectures. However, constraining formation processes requires mapping the observed exoplanet population to that which initially formed in the protoplanetary disc. Numerous results suggest that (internal or external) dynamical perturbation alters the architectures of some exoplanetary systems. Isolating planets that have evolved without any perturbation can help constrain formation processes. We consider the Kepler multiples, which have low mutual inclinations and are unlikely to have been dynamically perturbed. We apply a modelling approach similar to that of Mulders et al. (2018), additionally accounting for the two-dimensionality of the radius ($R =0.3-20\,R_\oplus$) and period ($P= 0.5-730$ days) distribution. We find that an upper limit in planet mass of the form $M_{\rm{lim}} \propto a^\beta \exp(-a_{\rm{in}}/a)$, for semi-major axis $a$ and a broad range of $a_{\rm{in}}$ and $\beta$, can reproduce a distribution of $P$, $R$ that is indistinguishable from the observed distribution by our comparison metric. The index is consistent with $\beta= 1.5$, expected if growth is limited by accretion within the Hill radius. This model is favoured over models assuming a separable PDF in $P$, $R$. The limit, extrapolated to longer periods, is coincident with the orbits of RV-discovered planets ($a>0.2$ au, $M>1\,M_{\rm{J}}$) around recently identified low density host stars, hinting at isolation mass limited growth. We discuss the necessary circumstances for a coincidental age-related bias as the origin of this result, concluding that such a bias is possible but unlikely. We conclude that, in light of the evidence that some planetary systems have been dynamically perturbed, simple models for planet growth during the formation stage are worth revisiting.


INTRODUCTION
Exoplanetary systems are diverse, with the constituent planets exhibiting a wide range of masses, semi-major axes and eccentricities (see e.g. Winn & Fabrycky 2015, for a review). Since the discovery of the first 'hot Jupiter' orbiting 51 Pegasi at a separation of ∼ 0.05 au (Mayor & Queloz 1995), it has been clear that planet formation and evolution processes do not exclusively, or even typically, result in Solar system-like architectures. A wide variety of processes to explain this diversity have been suggested, pertaining both to the primordial protoplanetary disc of dust and gas from which the planets form (Dullemond et al. 2007;Armitage 2011) and the subsequent evolution of the planetary system. These mechanisms are often strongly dependent on parameters that are poorly constrained empirically, such as disc viscosity (e.g. Bitsch et al. 2013a) or the initial conditions for exoplanets undergoing dynamical evolution (e.g. Barnes & Greenberg 2006;Carrera et al. 2019). It is therefore ★ andrew.winter@uni-heidelberg.de challenging to demonstrate the relative importance of particular mechanisms in reproducing the observed exoplanet distribution.
A key step towards understanding the origin of the observed exoplanets and system architectures is separating out processes that govern isolated planet-disc evolution from those that may perturb the planetary system during or after formation. An illuminating example is the dichotomy between the single Kepler planets and those found to be in multiple systems (Lissauer et al. 2011). The single systems have higher eccentricities than those in multiple systems (Xie et al. 2016), which may be the result of dynamical scattering and ejection or high mutual inclination of companions (Zhu et al. 2018). Similarly, whether dynamical perturbation drives the inward migration of hot Jupiters remains a topic of debate (Dawson & Johnson 2018). Their orbits may be the result of planet-planet (e.g. Weidenschilling & Marzari 1996;Johansen et al. 2012) or binary (e.g. Nagasawa et al. 2008;Belokurov et al. 2020) dynamical perturbation, or disc induced migration (e.g. Lin et al. 1996;Heller 2019). Complicating things further, stars form in clusters or associations (Lada & Lada 2003), and planet formation is subject to numerous influences from its formation environment (Parker 2020) that may leave an imprint on the planet population (Winter et al. 2020b;Dai et al. 2021;Adibekyan et al. 2021;Longmore et al. 2021). Given these considerations, efforts to perform population synthesis for the observed population may be confounded by the fact that we are not observing a continuous population. This work is motivated by the growing evidence that a subset of planetary systems may have been perturbed from their 'natural' formation configuration, and may therefore be confusing the interpretation of the exoplanet population as a whole.
In the remainder of this work, we explore the possibility that a simple model for the growth of planets may explain a subset of the observed exoplanet population that has not been perturbed. In Section 2 we outline the simplest theoretical expectation for a planet growing in a low viscosity disc. We then consider this limit in the context of two samples of planets that are the least likely to have undergone dynamical perturbation: Kepler multiples (Section 3) and the recently identified RV planets hosted by stars with low positionvelocity density (Section 4). In this way, we demonstrate that an upper limit in planet growth may explain the properties of unperturbed systems. We discuss whether such a limit is consistent with other observational constraints in Section 5, and draw conclusions in Section 6.

THEORETICAL MOTIVATION
A wealth of theoretical works have investigated the role of planetdisc interaction in inducing the migration of planets (Kley & Nelson 2012). The main empirical motivation for such theory is the existence of massive planets at short orbital periods (e.g. Mayor & Queloz 1995;Dawson & Johnson 2018). However, the rate and stopping mechanisms (e.g. Rice et al. 2008;Ida & Lin 2008;Matsumura et al. 2009), or even the direction (e.g. Michael et al. 2011), of migration are all sensitive to underlying assumptions. The outcome of migration is thus dependent on many parameters that have relatively poor empirical constraints, and affect both the initial conditions and disc physics. Mounting evidence suggests that some combination of environment (Brucalassi et al. 2016;Shara et al. 2016;Winter et al. 2020b;Dai et al. 2021), multiplicity (Belokurov et al. 2020;Hirsch et al. 2021), or orbital instability (Weidenschilling & Marzari 1996;Johansen et al. 2012), leaves an imprint on planetary systems. It is therefore possible that the torques exerted on a planet during formation do not lead to efficient migration (i.e. from several tens to ∼ 0.1 au). In this case, complex formation models applied to explain the population of observed exoplanets may not be necessary. To explore this possibility, it is instructive to consider the limit of no migration; a simple model in which the planet simply grows by accreting in its vicinity.
The 'isolation mass' is the mass to which a planet can grow by clearing out an annulus of material within the protoplanetary disc. This mass is determined by the local surface density Σ pl at the location of the planet (or planetesimal). While the isolation mass is usually considered in the context of initial accretion of solid material (dust), the same principles can apply to gas accretion if angular momentum transport in the disc is inefficient. This is the simplest picture of planet growth, where gas transport across the initial gap opened up by the nascent planet is negligible (i.e. the case of low viscosity). If the growth of a planet is not limited by the opacity of accreting material, then runaway growth is possible (e.g. Pollack et al. 1996). In this case, the isolation mass can be estimated by assuming that planetesimals initially follow circular orbits, and accrete all material approximately within their Hill radius H . The resultant isolation mass estimate is: where pl is the the semi-major axis of the forming planet. The annulus radius is approximated Δ ≈ H , where the constant ≈ 2 √ 3 (Lissauer 1993;Armitage 2007a). One then needs to make assumptions about the functional form of the surface density, which we will assume has the form (see discussion in Section 5.1.1): where subscript '0' values refer to normalisation constants in units of the associated variables. In this case the isolation mass is determined by the equation: with no analytic solution for general . For a viscous disc, the radial dependence of the disc surface density is dependent on the temperature profile. Since we are assuming little migration, we expect to be in the regime of low viscosity, where heating is radiative rather than viscous. For a constant viscous -disc that is not viscously heated, the mid-plane temperature is ∝ −1/2 (Kenyon & Hartmann 1987), for radius within the disc. In equilibrium, the corresponding surface density in the intermediate region between the inner and outer edges of the disc is Σ ∝ −1 ( = 1 - Hartmann et al. 1998 where we have chosen * ,0 = 1 , 0 = 1 au. However, equation 3 cannot apply throughout the disc, since as the semi-major axis pl becomes large, the fractional feeding zone radius Δ / pl also increases. When Δ → pl , the approximation must break down, and the feeding zone radius is no longer proportional to the Hill radius. We show Δ / pl as a function of pl in Figure 1, illustrating that we expect the small feeding zone approximation to break down in the range pl ∼ 1−10 au. Because we have assumed that surface density scales linearly with stellar host mass * , Δ / pl is independent of * . If planet growth is limited by the isolation mass we would expect an upper limit following equation 3 for planets with semi-major axes pl 1 au, transitioning to a shallower relationship at some larger pl .

UPPER MASS LIMIT FOR KEPLER MULTIPLES
Due to its well-characterised detection efficiency, data from the Kepler mission has provided invaluable constraints on the occurrence rates of exoplanets within ∼ 1 au of their host star (e.g. Borucki 2016). However, inferring occurrence rates from the full sample suffers from a potential limitation similar to that of population synthesis models; the implicit assumption that a single continuous sample of exoplanets are observed. In particular, the observed dichotomy of single and multiple systems in the Kepler sample (Lissauer et al. 2011) is strong evidence that the observed exoplanets form through multiple distinct channels. The higher eccentricities of single planets (Xie et al. 2016;Mills et al. 2019;Bach-Møller & Jørgensen 2021) suggests that an internal (Carrera et al. 2019) or external (Dai et al. 2021;Rodet et al. 2021;Longmore et al. 2021) perturbation has driven scattering and high eccentricities, either ejecting planets or yielding high mutual inclinations (Zhu et al. 2018). By contrast, multiple systems discovered by transit must have similar inclinations and are therefore most likely to have formed without a dynamical perturbation. To quantify the relative occurrence rates of planetary systems forming in isolation, it is therefore instructive to consider only the Kepler multiple systems.
Here we consider the probability density function (PDF) of exoplanet properties as a function of mass and semi-major axis. Our approach is inspired by that of Mulders et al. (2018, hereafter M18), with a number of important differences. First, we are not interested in absolute occurrence rates, only the relative distribution of detected exoplanets. We therefore need only consider confirmed exoplanet detections, as opposed to all Kepler targets. Related to this point, we restrict ourselves to Kepler multiple systems for the reasons outlined above. In addition, we do not assume that the two dimensional PDF is a separable function of period and radius. Specifically, we will ask whether an upper limit in exoplanet masses is statistically favoured over the separable broken power-law model adopted by M18. Following from this, the metric we adopt for comparison between observations and models accounts for the covariance of the period-radius distribution, rather than assuming that period and radius distributions can be compared independently. Finally, our PDF is defined in terms of mass and semi-major axis rather than radius and orbital period. This is significant for the functional form of our PDF, but we convert to period and radius using empirical constraints for comparison with the observed distribution (Section 3.3.3).
In the remainder of this section we describe the sample and cuts made in Section 3.1, then the model for the PDF in Section 3.3. The (relative) completeness of the Kepler sample is described in Section 3.2, and the models we evaluate are outlined in Section 3.3. The parameter space exploration is described in Section 3.4. Our results are summarised in Section 3.5.

Sample
Throughout this section we use the exoplanet and stellar host properties from the Kepler Data Release 25 (DR25), which is a uniform processing of the full four year dataset (Q1-17 -e.g. Twicken et al. 2016;Berger et al. 2018). These data were downloaded from the NASA Exoplanet Archive (NEA - Akeson et al. 2013) on the 29 th November 2020. As discussed above, we are interested in recovering the distribution of exoplanets in systems that are the least likely to have undergone a perturbation subsequent to formation. We therefore consider only systems with two or more confirmed transiting exoplanet detections, implying low mutual inclinations. This choice has the added benefit that false positives are unlikely. Since we expect an upper mass limit to be dependent on stellar mass, we wish to consider a narrow range of estimated host star masses 0.5−1.5 . We further exclude evolved stars with surface gravity log < 3.5 (where is in cgs units). The properties of the final host star sample are shown in Figure 2. It comprises 432 systems, of which 285 have two detected planets and 92 have three. The total number of planets in these systems is 1085.

Survey completeness
In order to draw comparisons between relative exoplanet frequencies in our simple model and those of the Kepler planets, we need an estimate of the relative fraction of planets at a given mass and semimajor axis that would be detected by the recovery pipeline -i.e. the relative false negative frequency. For a fuller discussion of the completeness considerations, see M18; we briefly summarise as follows.
The most obvious bias in the transit detection efficiency is the geometric probability that a planet will transit. For a circular orbit, this probability is: corresponding to reduced detection efficiencies at long orbital periods. At short orbital periods 3 days further incompleteness factor harm is due to the behaviour of the harmonic fitter, applied to remove periodic stellar activity in the Transiting Planet Search (TPS -Jenkins et al. 2010;Tenenbaum et al. 2012;Christiansen et al. 2013). We interpolate the false negative rates between 0.5−3 days estimated by Christiansen et al. (2015), adopting values of harm : 0.4 at 0.5 days, 0.6 at 0.8 days, 0.9 at 1.5 days, and 0.99 at 2.0 days, linearly interpolating between these values in this range. The completeness for a given data release from the Kepler satellite is limited at long periods by the requirement of three transits across the baselines -i.e. a maximum orbital period of two years for the entire four year dataset. In addition, for small planet radii pl , the chance of detection is limited by the signal-to-noise ratio. Factoring in these variable detection efficiencies, the completeness of the Kepler transiting sample can be quantified as a function of the Multiple Event Statistic (MES) used in the TPS algorithm. We    . Total detection efficiency for the host star Kepler-121. The colour bar shows the fraction of planets at a given orbital period (0.5−730 days) and planet radius (0.3−20 ⊕ ) and random inclination that would be considered a detection as a consequence of geometric effects and processing by the pipeline summarised in Section 3.2. Black points mark the measured periods and radii of Kepler-121 b and c. calculate the signal-to-noise detection efficiency ( S/N ) for each star in our sample using K PORT 1 (Burke & Catanzarite 2017). Finally, even detectable exoplanets are not necessarily vetted as planet candidates by the Kepler Robovetter (Coughlin et al. 2016). We use the parametric approximation for the vetting efficiency by M18, based on the simulated data products 2 and following Thompson et al. (2018): with = 0.63, = 0.19, break = 53 days, = −0.07, = −0.39.
The overall detection efficiency for a given host star is then given by: We show an example of the total detection efficiency for the host star Kepler-121 in Figure 3.

Broken power-law model
M18 made the simplifying assumption that the logarithmic PDF that describes the relative number pl of the Kepler planets at a given period and radius is a separable function of period and radius: This PDF is written as the product of two broken power-laws: for which the constants , , break , , and break are free parameters. M18 also estimated the absolute occurrence rate (normalisation), or number of planets per star. However, we are here only interested in the relative planet numbers, since the motivation for this work is that some subset of planets are a distinct population in the first instance. This leaves six parameters for which to fit.

Mass limited model
We consider an alternative to the above model: In this case, the PDF of planet masses follows a similar broken power-law as before: which is empirically motivated (Howard et al. 2012;Petigura et al. 2013), and possibly reflects two physical regimes in the growth of planets. However, the overall two-dimensional PDF is no longer a separable function of semi-major axis (period) and mass (radius). The PDF lim ( pl ; pl ) instead represents a 'cusp', or upper limit in planet mass pl that is dependent on their semi-major axis pl . Unlike the (double) broken power-law model, the model presented here does not explicitly treat the observed break in planet frequency at orb ∼ 10 days (Howard et al. 2012). However, the PDF is not separable into components that are functions of only pl (or pl ) and pl (or orb ). Thus a break in the PDF of pl also produces a break in the PDF of pl .
If disc viscosity is low and type II planet migration is inefficient, then in situ growth after gap opening eventually leads to the exhaustion of material that the planet(esimal) can easily accrete. As discussed in Section 2, such a mass can be approximated by estimating the quantity of disc material that orbits within the planet's Hill radius. Planet mass should therefore be limited by the isolation mass: In deriving the isolation mass in Section 2, we neglected the deviations in the disc surface density at the inner edge. Given that the Kepler planets are close-in, we might expect the surface density to steepen in this regime. We therefore write the mass limit: where in ∝ in , the inner edge of the disc. Strictly, is not a free parameter since = 1.5 is a consequence of the assumptions discussed in Section 2. We therefore initially consider only the normalisation constant 0 and in as free parameters, although we subsequently also fit for .
Defining the mass limit as in equation 13, we can impose a smooth decline in the probability of a planet having mass pl above this limit: where ln ( ; is the log-normal PDF in log-space. The dispersion also implicitly accounts for any dispersion between the mass and radius conversion. We have therefore introduced three additional free-parameters 0 , in and that describe the hypothesised upper limit in planet masses, for a total of six free parameters with , and break .

Conversion to observable quantities
The radius and period of a given planet are the quantities that can be inferred directly from transit measurements, while our model is a prediction for the mass and semi-major axis distributions. We therefore need to convert the physical quantities predicted by the model to the observable ones. To convert semi-major axes to orbital period, we simply adopt a stellar mass * = 1 . Converting mass to radius, however, requires appealing to empirical constraints.
The threshold of the high mass regime is chosen such that the relation inferred for the intermediate mass regime, intersects the shallow high mass relationship at pl = 12.1 ⊕ (Bashi et al. 2017). The PDF in log pl , log orb space is then: normalised over the domain of interest. Strictly speaking, the physical relationship is not a well-defined mapping from mass to radius. This is because the rocky (small radius) and the gaseous (larger radius) populations overlap in mass in the range pl ∼ 3−30 ⊕ (Otegi et al. 2020). At planet masses in this range, a planet may retain a (tenuous) gaseous envelope or this envelope may be photoevaporated (Owen & Jackson 2012;Owen & Wu 2013). The transition is responsible for the 'Fulton gap', an absence of planets with radii pl ∼ 2 ⊕ (Fulton et al. 2017). In principle, one could account for the overlap by estimating the ratio of rocky to gaseous planets. However, this would complicate the model, especially in the absence of any strong constraints for the ratio. We therefore adopt equation 16 in converting masses and radii in our fiducial model.
We investigate the influence of our choice of mass-radius relation by considering an alternative relation in which planets are predominantly rocky ( pl ∝ R 1 ) up to a greater mass than assumed in equation 16 ( pl > 3.11 ⊕ ). However, such an exercise requires revisiting the functional form of the mass-radius relation. While strictly equation 16 is not differentiable at the hard transitions between mass-radius relations, the coordinate transform in equation 18 is numerically calculable over our grid because the transitions occur at infinitesimal ranges in planet radius space. However, here we cannot allow a hard transition ( pl ∝ R 1 → pl ∝ R 2 ) because there exists a (non-negligible) discontinuity in the relation as a function of radius. This would result in an undefined PDF over the gap when applying the co-ordinate transformation (equation 18).
To account for the above numerical concerns, and in order to explore the influence of an alternative mass-radius relation, we adopt: where such that rg pl ( pl ) is continuous and differentiable over the rockygaseous transition. Here rg and Δ rg are parameters that set the location and width of the transition from rocky to gaseous massradius relations respectively. In principle, we could allow rg and Δ rg to vary as free parameters. However, we will find that the best Fiducial mass-radius relation Alternative relation Figure 4. The factor applied to the two-dimensional PDF in period-radius space (equation 18) using the two different planet mass-radius relations we adopt, as a function of planet radius and in arbitrary units. For the fiducial relationship (equation 16, solid line), this factor is proportional to the inverse of the power-law index for the local relation R 1,2,3 ( pl ) in radius space. By contrast, the alternative relation (equation 19, dashed line) incorporates a smooth transition between rocky and gaseous relations at higher planet mass (radius). fitting fiducial model with six free parameters is indistinguishable from the observed data by our comparison metric. It is therefore not appropriate to attempt fitting two additional parameters. Instead we fix rg = 10 ⊕ and Δ rg = 0.1, resulting in a transition between rocky and gaseous planets at approximately the mid-way mass in the overlap in the R 1 and R 2 regimes. While this choice is semiarbitrary, it serves the purpose of investigating the influence of our choice of mass-radius relation on our fitting results. In Figure 4 we show the factors applied to the two-dimensional PDF due to the adopted mass-radius relations.

Sampling
To assess the likelihood of a given model described by parameters Θ, we need to compare to the observed sample . To achieve this, we first define the two dimensional model PDF across a grid in orb − pl space. We choose a grid log-uniformly spaced between 0.5−730 days and 0.3−20 ⊕ at resolution 200 in period and 150 in radius. We draw a sample of planets˜a ll from this grid. To ensure synthetic quantities are not fixed to a grid cell value, we offset each period and radius by an amount drawn from a uniform distribution between plus/minus half a grid cell. We thus obtain a period and radius for each synthetic planet. To directly compare the model results with the observations, we must now account for the Kepler selection biases. We select a subset of these planets by Monte Carlo sampling from˜a ll with probabilities proportional to detection efficiencies at , . The specific detection efficiencies D for each synthetic planet are drawn at random from the pre-computed grid of the observed sample of host stars. This grid has the same resolution as the model PDF. We then draw a random number ∼ U (0, 1) for each synthetic planet, and consider it to be 'detected' for < D /maxD i . This yields a new sample˜d et . Unlike M18, we are aiming only to quantify relative occurrence rates, and therefore do not need to reproduce the absolute number of 'detected' planets. We are therefore free to normalise the D by the maximum detection efficiency maxD i , in order to maximise the size of the detected sample |˜d et |.
In order to assign a likelihood (Section 3.4.2 below), we require total size of the comparison sample to be sufficient such that no significant stochastic differences between distributions occur between realisations for fixed Θ. We aim for a sample size of 3000, and therefore randomly draw this many from˜d et to give our comparison sample˜. We must therefore choose the size of our initial sample |˜a ll | such that |˜d et | > 3000 for the majority of Θ. We find that |˜a ll | = 1.2 × 10 5 is nearly always sufficient. In rare cases where it is not, we draw˜with the same size by allowing multiple drawings of the same synthetic planet . This should not significantly change our test statistic except in cases where the fraction of detections are so low that any such Θ would in any case be highly disfavoured.

Estimating likelihood
We require a statistical test to measure the agreement of the synthetic sample of periods and radii˜obtained from our model with parameters Θ with the observed data . M18 estimate the likelihood using Fisher's method, by which one can combine likelihood estimates on independent comparison quantities to estimate an overall likelihood (Fisher 1925). However, performing comparison tests independently for radius and period neglects any possible covariance between the two. We address this problem by performing a single statistical test quantifying the (dis-)similarity of˜and in orb − pl space.
One can demonstrate differences between two two-dimensional datasets by generalising the Kolmogorov-Smirnov (KS) test. The approach for this was proposed by Peacock (1983), the philosophy of which was to replace the maximal difference between cumulative distributions with the analogous difference across all quadrants in a dataset defined at each data point. While such a metric is not distribution-free, unlike the one-dimensional case, it was shown by Fasano & Franceschini (1987) that it is nearly distribution-free if the two distributions have a similar (Pearson's) correlation coefficient . Press et al. (2007, see Section 14.8) used the results of Fasano & Franceschini (1987) to infer an approximate formula for the test statistic in two-dimensions. The corresponding probability of obtaining two samples from the same distrubution is: is the effective sample size: and is the complement of the cumulative distribution function of the  Figure 5. Example of the two-dimensional PDF in orbital period and radius resulting from the mass limited model. Red points are the observed periods and radii of the Kepler multiples in our sample, faint yellow points are all planets drawn from the model PDF (the value of which is indicated by the colour bar), and blue points are those that are 'detected' from this synthetic sample. The parameters are the median values from the posterior distributions obtained in the F -PL model (in Table 1), with fixed power-law index = 1.5. The solid green line is the upper limit on the planet masses (defined by mass normalisation 0 and inner semi-major axis in ), with one deviation ( ) above this line indicated by the dashed green line.
KS distribution (P KS ). Since the value of can be defined at vertices drawn from either˜(˜) or ( ), in practice the mean = (˜+ )/2 is used. Equation 21 is a good approximation for 20 and KS,2D < 0.20. The former is no problem in this context. The latter simply indicates that two samples are not significantly different, and in our case that there are not sufficient constraints from the data to choose between model parameters. In fact, in Section 3.5.1 we do find that the upper limit model, for a finite range of Θ, is not significantly different from the observed distribution. We discuss the implications in that section.
We apply the Markov chain Monte Carlo (MCMC) implementation E (Foreman-Mackey et al. 2013) to explore the multidimensional parameter space. For observed data , Bayes' theorem states that the likelihood of a model with parameters Θ is: the probability of obtaining the data from the model. In this case, the -value from a two-dimensional KS test is the estimated probability of obtaining the observed data given that it comes from the same parent distribution as the synthetic data (i.e. the model). Hence ( |Θ) ≈ KS,2D , as defined in equation 21. Adopting uniform priors within a given range, we have: Equation 26 defines the likelihood we adopt for the MCMC parameter space exploration in Θ.

Best-fitting parameters
The posterior distributions obtained from the MCMC exploration of our mass limited model are presented in Appendix A, and the results are summarised in Table 1. Initially we run a model for which we fix the power-law index = 1.5, which would be expected from the predictions of the isolation mass limited growth (Section 2). This model (F -PL) with parameters as in Table 1 is illustrated in Figure 5. Good agreement is evident across the entire parameter space, with the possible exception of at large radii pl ∼ 12 ⊕ . This is where the PDF value is strongly sensitive to the assumed massradius relation at high masses. For convenience, we have chosen a hard transition at pl = 12.1 ⊕ to a nearly flat relationship (equation 16). Smoothing this transition would yield a less abrupt increase in the PDF just above the transition radius, and thus reduced radii for a subset of the synthetic detections at the largest periods and radii. However we do not attempt to address this failing in the model because such planets account for a small fraction of the total. Indeed, by the metric of the two-dimensional KS test the observed period-radius distribution is statistically indistinguishable from that of the model, with KS,2D ≈ 0.178.
This KS test -value means that we are data-limited in determining shape of the posterior distribution rather than modellimited. This means that more restrictive constraints on the model parameters requires either alternative empirical constraints on the physical parameters (priors) or a larger sample size. It also means that adding parameters cannot improve the model by any information criterion that introduces a penalty for degrees of freedom (e.g. Bayesian/Akaike information criterion -BIC/AIC).
The analogous posterior results for the alternative mass-radius relationship we discuss in Section 3.3.3 (equation 19) are labelled F -PL-in Table 1. The KS test result indicates a worse fit for the F -PL-model with respect to the F -PL model. This may be because our choice of mass-radius relationship in the former model is not an empirical fit. However, we are mainly interested in whether an alternative prescription for the transition between the rocky and gaseous planets alters our results. In this regard, we find marginally different best-fitting parameters describing the planet mass-function. However, the main parameters of interest are those that describe the hypothesised upper mass limit. These are consistent with the values obtained in our fiducial (F -PL) model. We conclude that the form of the inferred upper limit is not strongly dependent on the choice of mass-radius relation.
While a fixed = 1.5 is consistent with the observational con-    Figure 7 but for the BPL-PR model that is a separable PDF in period and radius. Such a model is similar to that of M18, with best-fit parameters given in Table 2. straints, we can ask the degree to which this index is statistically favoured. We may therefore attempt a further parameter exploration allowing to vary. However, given the model outcome is indistinguishable from the data for fixed = 1.5, allowing an extra parameter will always yield over-fitting and poor constraints on parameters. In particular, we find that allowing the inner edge in to vary (particularly in 0.1 au) results in strong degeneracy between the parameters describing lim , and poor constraints as a result. We therefore fix log in /1 au = −1.4 before repeating the MCMC exploration with ≡ arctan as a free parameter. This choice is semi-arbitrary: physically we expect an inner edge at a significantly greater extent than the stellar radius, and the posterior distributions are not strongly dependent on the adopted value of in for in 0.1 au.
The results of fitting for are summarised in Table 2, with the posterior distribution shown in Appendix A. We first find that the parameters describing , the broken power-law describing the planet mass PDF, are not significantly affected by allowing a variable . We also find that the power-law and normalisation are consistent with those obtained by fixing = 1.5. The relatively poor constraints on these values are due to the strong covariance between the parameters describing the upper limit (see Figure A3). In particular, large is associated with large and 0 . We are limited in that Kepler is only sensitive to planets with short orbital periods, such that the range in period (semi-major axis) space is prohibitive. Nonetheless, the constraints in the full parameter volume space are such that further constraints on one of these parameters would allow much better constraints on the others.

Comparison to broken power-law model
Fitting the double broken power-law model yields best-fitting parameters shown in Table 2, with similar results as those found by M18. An MC realisation for the best fitting parameters is shown in Figure 6. The results of the two-dimensional KS test result for the maximum likelihood parameters yield ML KS,2D ≈ 4 × 10 −4 , indicating a worse fit than the upper limit model. The upper limit model is therefore the preferred description of the two-dimensional distribution in orb − pl space.
In order to draw more qualitative comparisons between the two classes of model we compare the two-dimesional PDF that would be inferred if data was drawn from each model to that inferred directly from the real data. We first draw a single Monte Carlo realisation mod of each model. We then estimate the two-dimensional PDFs from the datapoints for both model and observations ∈ , where = 1, . . . , | | and = mod or obs , the observed data. We achieve this for each dataset by using a Gaussian kernel density estimate (KDE): where is the covariance matrix and G is a Gaussian kernel with unity variance and zero mean. The shape of this PDF is dependent on the bandwidth or smoothing length scale, ℎ . Here we follow Scott (1992) in adopting: where = 2 is the number of dimensions (period and radius). We apply the S (Virtanen et al. 2020) implementation gaus-sian_kde to compute the KDE here and throughout this work.
The Gaussian KDEs for the Kepler multiples and the fixed power-law upper limit (F -PL) model are shown in Figure 7. Certain regions of orb − pl space are poorly sampled, such that comparisons to the model are of less significance. As a guide, we indicate in the bottom panel the contour along which the geometric mean of the observed and model KDEs drops to √︃ˆm odˆobs = 0.05. Comparisons between the model and observed distributions should be made inside the region bounded by this contour. In general we see good agreement between the two models, including a dearth of planets at large radii and short orbital periods, and a feature at the transition between gas-rich and terrestrial planets. As already suggested, the model slightly overestimates the radii of the most massive planets, which is a consequence of the adopted functional form of the mass-radius relation (equation 16).
The F -PL model can be compared to the results of the double broken power-law model (BPL-PR) in Figure 8. In the latter case, the model significantly overestimates the number of planets with large radii ( pl 4 ⊕ ) and short orbital periods ( orb 10 days) compared with the observed distribution. This can be seen immediately from the KDE, but is made clearer by comparing the bottom panels of Figures 7 and 8. In addition, the feature at the transition from terrestrial to gaseous planets is seen as a blue blob in Figure 8, indicating an underestimate of the observed PDF by the model. These mismatches indicate why the broken power-law model is disfavoured by the data with respect to the upper limit model.

Summary
We have demonstrated that an upper limit in planet masses can explain the distribution of the orbital periods and radii of the Kepler multiple systems. By the metric of a two-dimensional KS test, we show that such an upper limit can yield a sample of planets that are statistically indistinguishable from the observed population. This model represents a significant improvement over a double broken power-law model describing both the period and radius distributions independently. Since the model is not distinguishable from the data, no model that invokes a greater number of free parameters can be favoured by an information criterion that introduces a penalty for degrees of freedom. Hence, this model represents the best possible description of the data in the absence of one or more of the following: (i) A greater sample size (in particular, with a greater range in orbital period).
(ii) Stronger constraints on physical parameters. (iii) An alternative comparison metric.
(iv) A simpler model.
It is unlikely that (i) with such well characterised detection efficiencies as Kepler will become available in the near future. More imminently, (ii) may be possible by appealing to protoplanetary disc physics and demographics, but we do not attempt to do so here. An alternative comparison metric -(iii) -may also be possible, although any such metric should account for the two-dimensionality of the data as we have done here. Finally, while (iv) is possible, our model is already simplified by assuming that difficult-to-treat mechanisms (e.g. migration) do not influence the planet population.

LONGER PERIOD PLANETS
When considering whether an upper limit may also describe a subset of systems at higher mass and longer period, we are not able apply a similar approach to Section 3. This is both because of the massradius degeneracy for massive planets and because Kepler planets are only confirmed if they transit three times during the four year baseline, such that a hard upper limit of two years exists in orbital period. Although samples of ∼ 100s of radial velocity (RV) detected planets with homogeneously analysed stellar properties exist (e.g. Fischer & Valenti 2005), no sample with a single discovery instrument currently compares to Kepler in both the number of planets and the characterisation of the detection efficiency. Constraining a two-dimensional model PDF in the same way is therefore not feasible for the longer period planets. We therefore take an alternative approach.

Motivation and sample
For inner planets, scattering by gravitational interaction with an exterior body most frequently results in angular momentum loss, and therefore the decrease of the semi-major axis (e.g. Heggie . Such inwards migration may result from system instability and planet-planet scattering, which may be the origin of hot Jupiters (e.g. Carrera et al. 2019). These dynamical instabilities may develop in an isolated system, but scattering can also be induced as a result of an external perturbation during the passage of a neighbouring star in an open cluster (e.g. Shara et al. 2016;Cai et al. 2017, and discussion in Section 5.3.1). If such perturbations produce a significant fraction of detected short period massive planets, and there exists a physical upper limit for planet growth as we have inferred in Section 3, we would then expect planets that form in high density environments (HDEs) to preferentially be found above it with respect to those born in low density environments (LDEs). Kinematic signatures of a star's formation environment may persist for ∼Gyr time-scales. This possibility is supported by the chemical-kinematic correlations between neighbouring stars and simulations tracing co-moving pairs of stars (Kamdar et al. 2019a,b;Nelson et al. 2021), as well as recent results indicating that dynamical heating is inefficient for open clusters (Tarricq et al. 2021). It may therefore be possible to identify exoplanet hosts that formed and evolved in LDEs by their kinematics. Winter et al. (2020b) recently delineated exoplanet hosts by position-velocity density and found significant differences between architectures of systems around hosts in the age range 1 − 4.5 Gyr. Systems around lower density hosts may be associated with a higher multiplicity rate among Kepler planets (Longmore et al. 2021), hinting that these planets are less likely to have experienced dynamical perturbation. Here we compare the HDE and LDE population of planets around ∼ 1 mass stars discovered by RV surveys. We restrict the sample to those with masses 0.7−1.5 (as quoted in the NEA) such that they are close to 1 (within ∼ 0.15 dex). While the samples are heterogeneous such that occurrence rates cannot be directly calculated (although see also Dai et al. 2021), the approach for delineation does not yield systematic differences in the observational biases between each sample (i.e. no bias of biases). Therefore, in principle differences between HDE and LDE systems can be studied (although see discussion regarding stellar ages in Section 4.3). To maximise detectability, we consider planets with masses pl sin > 300 ⊕ leaving a sample of 68 HDE and 22 LDE planets. To rule out the possibility of tidal inspiral affecting our sample (see Section 4.3) we additionally exclude all the (massive) planets with semi-major axes pl > 0.2 au, yielding 61 and 20 planets respectively.

High density versus low density environment comparison
We wish to compare both HDE and LDE planet populations to the hypothesised limit inferred from Kepler multiples. For any planet, we can define the logarithmic deviation from this limit: ,Kep ≡ log pl lim (29) such that planets below the limit have ,Kep < 0 and planets above have ,Kep > 0. A caveat to this approach is that we take the results for the fixed power-law index ( = 1.5) fit, since otherwise the location of the limit extrapolated to wider separations is subject to large uncertainties.
We show the result of computing the ratio of planet masses to´1 the Kepler multiple limit in Figure 9. The LDE planets preferentially have masses close to (or below) the extrapolated limit from the Kepler multiples with moderate significance KS = 1.2 × 10 −2 , suggesting that this growth limit for unperturbed planets persists out to larger separations. The distribution of masses above the Kepler limit appears as an almost continuous distribution for the HDE planets, while the two outlying hot Jupiters around LDE host stars also have outer companions (HIP 91258 b - Moutou et al. 2014-and HD 68988 b -Wright et al. 2007). These hot Jupiters may therefore also have reached their present orbits due to a dynamical perturbation, as suggested by the enhanced multiplicity fraction of massive hot Jupiters (Belokurov et al. 2020). Indeed, among lower mass hot Jupiters (50 ⊕ pl 1 J ) that are not associated with enhanced binary fractions, the RV HDE sample includes ten, versus none in the LDE sample.
To mitigate the possibility of tidal inspiral influencing our findings (Hamer & Schlaufman 2019), we show the results for the sample with pl > 0.2 au in Figure 10. We find a similar trend, with slightly greater significance KS ≈ 8 × 10 −3 . This suggests that the findings are not driven purely by an increased frequency of hot Jupiters among the dense population.
To illustrate our findings, we again consider the KDEs for our two samples. We define: whereˆl andˆh are the Gaussian KDE in pl − pl sin space for the planets around LDE and HDE stars respectively. The results of this calculation are shown in Figure 11. It is apparent that a higher fraction of the planets around stars in LDEs fall close to or below the upper limit inferred from the Kepler multiples, and extrapolated out to greater separations (black line). We perform a two dimensional KS test as before, finding marginal significance. However, since we are specifically interested in the deviation from the Kepler-inferred limit, a test across the orthogonal axes is not the best demonstration of the differences in this case. We do not find any significant difference in the orbital eccentricities, as defined in the NEA, between the HDE and LDE samples (both with median ≈ 0.2). Heterogeneous determinations may influence this finding by effectively introducing noise to the distributions.

Possible correlations with stellar age
Dynamical heating of stars yields lower phase space densities over time. Hence a sufficiently large random sample of stars with accurate stellar age determinations should exhibit differences in age when divided into high and low position-velocity density. In this case, potential covariance between stellar age and detected planet properties must be considered (see discussion by Mustill et al. 2021). Winter et al. (2020b) used a sample that had nominal ages from the NEA between 1−4.5 Gyr, such that the HDE and LDE samples did not have significantly different age distributions. However, stellar ages, even when homogeneously determined, are often challenging to accurately assess. For example, using a small set of RV host stars in a slightly larger age range (0.5-5 Gyr), Adibekyan et al. (2021) demonstrated LDE planets have longer orbital periods with respect to HDE planets (cf. Figure 11). However, they also found different age distributions when adopting the theoretical isochrones by Bressan et al. (2012). In general, difficulties in accurately determining stellar ages, particularly when using metallicity or stellar kinematics that may exhibit covariance with birth environment, mean that the possibility of systematic differences as a driver of the apparent difference in exoplanet system architectures must also be considered.

Physical effects
Physically, one way in which age may play a role in determining exoplanet properties is due to the tidal inspiral of hot Jupiters. Hamer & Schlaufman (2019, see also Mustill et al. 2021) found hot Jupiters are more common for dynamically colder (high density) host stars, but concluded stellar age, and not environment, was the most likely origin for this difference. However, this finding does not persist for lower mass short period planets (Hamer & Schlaufman 2020). This implies either a specifically tuned and planet mass dependent tidal dissipation efficiency, or that tidal inspiral is not the origin for the stellar kinematic dependence of hot Jupiter incidence. Arguing against tidal inspiral, hot Jupiter rates are enhanced and not depleted in the dense and old (∼ 4.5 Gyr) open cluster M67 (Brucalassi et al. 2016). In addition, while Adibekyan et al. (2021) are able to distinguish somewhat older ages for LDE stars, and fractionally fewer short period planets, they do not find differences in ages between stars hosting long and short period planets. This suggests a stronger covariance of planet architectures with dynamical properties than with (measured) age.
The tidal dissipation time-scale for orbiting planets is a strong function of semi-major axis (Adams & Laughlin 2006), such that cold or warm Jupiters would not be expected to inspiral during their lifetimes. We have excluded hot Jupiters in our comparison, effectively ruling out tidal inspiral as a factor in depleting the short period planets. We recover a similar finding that massive planets on short period orbits are associated with high phase space densities. This suggests that tidal inspiral is not responsible for our findings.
An alternative age-related influences on a planetary system orbiting at several au would be planet-planet scattering or Kozai-Lidov oscillations (e.g. Naoz et al. 2011). If systems are frequently formed in a state of marginal dynamical stability, the probability of scattering may increase with the age of the system. Scattering would be expected to produce planets on short period orbits. Therefore younger (possibly HDE) planets should be found with longer orbital periods with respect to older planets. This is the opposite of our findings, and therefore not a convincing explanation for the differences between the LDE and HDE planets.
We conclude that age is physically unlikely to be the origin of the differences for massive planets that are not hot Jupiters. If our findings do have a physical origin, then we might also infer that hot Jupiters are an extension of the (scattered) short period planet population. In this case, the enhanced hot Jupiter rate around dynamically colder (high density) stars (Hamer & Schlaufman 2019) could be due to slower dynamical heating of open clusters (Tarricq et al. 2021) rather than tidal inspiral. Some combination of both scattering and inspiral remains possible.

Observational effects
It remains possible that our results are affected by a correlation between age and the detectability of planets. Stellar activity can result in a jitter that can make detection and confirmation of true planet signals challenging. Given uncertain ages, this may influence our findings. The jitter amplitude at visible wavelengths for FGK stars of age ∼ 600 Myr is ∼ 13−16 m s −1 (Paulson et al. 2004;Quinn et al. 2012), and an order of magnitude greater than this for ages ∼ 100 Myr (Paulson & Yelda 2006). In the age range 1−4.5 Gyr, decreases comparatively slowly from ∼ 10 m s −1 to ∼ 5 m s −1 (Hillenbrand et al. 2015).
An age-related bias is more challenging to categorically rule out, both because of the aforementioned issues determining stellar ages and because per-star completeness among RV planet hosts have not been quantified in the same way as for the Kepler stars. However, age-related biases are unlikely to be the origin of our results for several reasons: • The semi-amplitude of the RV signal for a planet on a circular  Figure 11. The distribution of planet semi-major axes ( pl > 0.2 au) and masses ( pl sin > 300 ⊕ ), divided into LDE (blue) compared to HDE (red). The colour bar represents the result of subtracting the HDE KDEĥ from the LDE KDEl,ˆ=l −ĥ. The black line is the extrapolated best-fitting upper limit inferred using the F -PL model in Section 3. The cyan line traces the contour of fixed RV variation semi-amplitude Δ at three times the RV jitter due to stellar activity for a 1 Gyr old ( = 10 m s −1 ) FGK star. High density Low density We have excluded planets with masses < 1 J such that the ratio of the semi-amplitude Δ / > 3 nearly everywhere in the parameter space for stars older than 1 Gyr (see cyan line in Figure 11).
• The transition between HDE and LDE regimes in Figure 11 is coincident with the Kepler inferred limit and not coincident with contours of constant Δ (i.e. pl ∝ 1/2 pl ). • The difference between the HDE and LDE samples appears to be an absence of easily detectable massive planets on short periods in the LDE sample, rather than a lack of difficult to detect low mass planets on long period orbits in the HDE sample (i.e. the red region in Figure 11 lacks blue points, but the blue region does not lack red points). In the absence of a comprehensive list of RV target stars in the age range 1−4.5 Gyr, we cannot directly calculate the absolute detection rates in each sample. In principle, our findings could therefore be due to a greater number of RV target stars in HDEs than LDEs. However, for a given total number of 'easily detectable' planets det , we can estimate the probability of lde of these detections being in the LDE sample, given the fraction lde of stars with a density categorisation being low density (equal numbers of HDE and LDE RV target stars corresponds to lde = 0.5). For equal occurrence rates, the probability of obtaining lde LDE planets or fewer in this domain is: In the red region of Figure 11, which we assume is above some agerelated detectability threshold, we have det = 36 and lde = 2. The total fraction for all stars in the solar neighbourhood varies between ld ≈ 0.3−0.6, for which equal occurrence rates can be rejected with null 3 × 10 −4 . Note that this value is smaller than the previously quoted null probabilities because we are here arguing from absolute numbers across samples, rather than relative numbers within each sample. Uncertainty is introduced to this test because it depends on the true lde within the superset of RV target stars. However, RV surveys preferentially select inactive stars, hence any bias should be towards older, possibly lower density target stars (increased lde ). Our results therefore suggest an absence of short period high mass planets around LDE stars, rather than a bias-related reduction in the fraction of detections around HDE stars along the Kepler limit.
• We find no significant correlation across the sample between the NEA age and the semi-major axis (Pearson = 0.19 and = 0.10) or planet mass ( = −0.10, = 0.38). The age estimates are also similar across both samples (Figure 12), with no obvious excess at young ages in the high density population. We emphasise that these ages are heterogeneously determined, and thus may be subject to large uncertainties (see discussion by Adibekyan et al. 2021).
Despite these arguments, we do not rule out an age bias contributing to our results. Unfortunately, the sample required to do so directly is a large catalogue of RV target stars with characterised detection efficiencies and ages in the range 1−4.5 Gyr, along with tens of detected planets out to several au. While this is not realistic, an alternative may be careful modelling of the differential detection efficiency with age at long periods. This is beyond the scope of this work. Nonetheless, the absence of planets around low density stars above the Kepler upper limit is an intriguing coincidence, and the expected result if such a limit were physical.

DISCUSSION
In this work, we show that Kepler multiples are well-described by an upper limit that resembles the expected isolation mass for a low viscosity disc. Planets around stars with low position-velocity space appear to preferentially follow this Kepler inferred limit, with few planets found above the hypothesised limit compared to high density stars. These findings suggest that planet growth in isolation yields isolation mass planets, and high mass planets on short orbital periods are the product of dynamical scattering.
In this section we discuss the observational constraints from the literature that support or contradict the hypothesis that planet mass without any dynamical perturbation is limited by the local surface density during the protoplanetary disc phase. Assuming mass limited growth, we also discuss the processes that may result in the observed exoplanet architectures.

Disc surface density, viscosity and feeding zones
We first consider whether the surface density that would yield an upper limit comparable to that obtained in Section 3 is consistent with observations. For a solar mass star, the surface density we have assumed (see Section 2) has the form: where we take = 1, appropriate for a radiatively heated disc with constant viscous (Shakura & Sunyaev 1973;Hartmann et al. 1998). From the upper limit in planet masses we have inferred and equation 3, we can estimate Σ 0 ≈ 1.8 × 10 3 g cm −2 . The resultant disc mass is then: for outer disc radius out and formation efficiency form ≤ 1. If form ∼ 1, then this is broadly consistent with masses of young protoplanetary discs (e.g. Eisner et al. 2005;Jørgensen et al. 2009;Andrews et al. 2013;Tychoniec et al. 2020). We can further ask if the value we adopt for the gas disc surface density power-law index = 1 is empirically supported. This value is the consequence of a constant viscous , with mid-plane temperature ∝ −1/2 . Broadly, a = 1 power-law out to some characteristic radius is consistent with the surface density profiles inferred from (sub-)mm continuum observations of protoplanetary discs, which trace the dust (Andrews et al. 2012;Andrews 2015). However, such constraints are somewhat complicated by the fact that we are primarily interested in the gas surface density profile, which dominates the mass budget. A variable dust-to-gas ratio, dust migration and grain growth (Birnstiel et al. 2010;Sellek et al. 2020), as well as related opacity effects (e.g. Rosotti et al. 2019;Ribas et al. 2020) complicate the inference of the gas surface density from such measurements. In addition, detected planets are almost exclusively at semi-major axes < 10 au, while most surface density constraints probe material at several tens of au separations from the host star that may follow a different surface density profile. Nonetheless, = 1 is consistent with a number of observed surface densities in the dust (e.g. Tazzari et al. 2016;Fedele et al. 2017) and the gas (e.g. Huang et al. 2016;Cleeves et al. 2016). Recent resolved studies of the gas temperature profiles have also found agreement with the canonical power-law temperature profile (Calahan et al. 2021).
At separations pl 4 au, if planets are allowed to grow to isolation as following the relationship inferred in Section 3, their masses would exceed 4 J . Such planets may form from fragmentation rather than core accretion (for a recent review see Kratter & Lodato 2016), similar to brown dwarfs (Chabrier et al. 2014;Schlaufman 2018). At these distances one would no longer expect planet formation to be limited by the isolation mass, since growth to this limit would mean consuming the majority of the disc mass (see discussion in Section 2).

Disc induced migration and eccentricity
If the growth of exoplanets inside ∼ 1 au is limited by the isolation mass, this suggests that efficient type II migration (once a gap opens in the disc) does not operate in this regime. While migration models may be able to produce observed systems in conjunction with other mechanisms (Armitage et al. 2002;Ida & Lin 2008, these models are all highly dependent on parameter choices and initial conditions (e.g. disc surface density, viscosity, metallicity and the planet formation position and timescale -Armitage 2007b; Alexander & Armitage 2009; Mordasini et al. 2012). In addition, problems persist in explaining both hot and warm Jupiters by disc induced migration, particularly tuning the required braking mechanisms to reproduce both the mass and period distribution of massive planets (Hallatt & Lee 2020).
While our results suggest type II migration is not the origin of high mass, close-in planets, it does not preclude the action of type I migration throughout the disc or type II at larger separations. There may be special locations within a disc that curtail migration such as dead zones (Matsumura et al. 2009) and ice-lines (Ida & Lin 2008). If migration is inefficient within the ice-line, this could be consistent with the apparent accumulation of planets close to this separation (Fernandes et al. 2019). It is possible that planets that form inside or around the ice-line have their growth limited by the isolation mass, while planets that form further out at several au are still capable of significant disc-induced migration.

Hot Jupiters around young stars?
If the growth of inner planets is in situ, and limited by the isolation mass, it might be expected that hot Jupiters around young stars should be rare. It is challenging to (unambiguously) detect planets around the youngest stars in the plane of a disc by the light curve during transit. Detection and characterisation using RVs is similarly challenging since magnetic activity can induce signals that are orders of magnitude greater than line of sight velocity variations. This means that, even once detected, constraints on planetary properties are often weak.
Several studies have claimed the discovery of exoplanets in open clusters and associations (Quinn et al. 2012(Quinn et al. , 2014Malavolta et al. 2016;Mann et al. 2016Mann et al. , 2017Rizzuto et al. 2020). The existence of hot Jupiters in open clusters has been used to suggest that hot Jupiters either form in situ or migrate though planet-disc interaction. However, hot Jupiters found in these environments are older than ∼ 600 Myr such that tidal circularisation is possible. RV surveys of younger open clusters ( 100 Myr) have so far yielded no hot Jupiter discoveries (Paulson & Yelda 2006;Bailey et al. 2018;Takarada et al. 2020).
For individual very young (few Myr old) stars, follow ups on hot Jupiter candidates have cast doubt on a number of detections. One of the greatest challenges to the upper limit hypothesis we have explored in this work would have been the hot Jupiter detected around the 2 Myr CI Tau (Johns-Krull et al. 2016). This detection would be particularly problematic due to the extended protoplanetary disc hosting an ensemble of gas giants out to ∼ 100 au . It would therefore be difficult to explain the hot Jupiter as a product of a dynamical perturbation. However, Donati et al. (2020) presented observations constraining the magnetic activity of CI Tau, showing that the rotation period is equal to that of the hypothesised hot Jupiter (9 d), and that the field is sufficient to produce a magnetospheric gap extending to close to the corotation resonance. This casts doubt on the origin of the previous detection. Similarly, Damasso et al. (2020) were recently unable to detect the previously claimed hot Jupiter around the ∼ 2 Myr V830 Tau (Donati et al. 2016(Donati et al. , 2017 in HARPS-N data, despite performing a wide range of tests. There remain two examples of hot Jupiter detections around weak-lined T-Tauri stars that have not yet been called into question: one around the 17 Myr old TAP 26 (Yu et al. 2017) and one around the 5−10 Myr K2-33 (David et al. 2016). TAP 26b is detected by radial velocity variations, requiring careful modelling to infer that differential stellar rotation cannot be the origin for the signal. The period is ∼ 10 days, but could not be uniquely determined from the data, and the minimum mass is 1.66 ± 0.31 J . K2-33b is inferred by transits of period 5.4 days, and has a radius of 5.8 ± 0.6 ⊕ . Only upper limit constraints exist on the mass this young planet, and the radius is marginally consistent with the upper limit we infer in Section 3. For neither planet does the available data offer (strong) constraints on the orbital eccentricity. If real, TAP 26b in particular is a challenge to any formation model with no efficient type II migration through the disc. However, since no substantial disc remains and there is currently no suggestion of a stable system of exterior planets, in principle migration through the disc initiated by scattering remains possible. If dynamical scattering occurs during the disc phase, the subsequent orbital evolution of a massive planet can be unpredictable, and one must then consider how growth, migration and planet-disc interactions alter the orbit of such a planet (e.g. Papaloizou & Terquem 2001;Bitsch et al. 2013b;Ragusa et al. 2018).

ALMA 'planets' in discs
Studies of dust in protoplanetary discs have demonstrated that structure is common (e.g. ALMA Partnership et al. 2015;Andrews et al. 2018;Clarke et al. 2018). Many of these structures are dips in the dust surface density in rings that may be consistent with the clearing of material by planets (e.g. Zhang et al. 2018). Although other explanations for these rings are possible (Zhang et al. 2015;Flock et al. 2015;Gonzalez et al. 2015;Okuzumi et al. 2016), many discs also exhibit kinematic deviations from Keplerian velocity that are consistent with presence of protoplanets (Casassus & Pérez 2019;Pinte et al. 2020).
The best constraints on planet growth and evolution in a disc come from the two accreting protoplanets inside the transition disc around PDS 70 (Hashimoto et al. 2012;Keppler et al. 2018;Haffert et al. 2019). The resonant configuration and accretion rates of these planets can be explained by planet-disc interaction and migration (Toci et al. 2020). PDS 70 b and c have orbital separations of ∼ 20 au and 35 au and masses of ∼ 4−17 J and ∼ 2−10 J respectively (Müller et al. 2018;Haffert et al. 2019;Christiaens et al. 2019;Bae et al. 2019). Such planets (or sub-stellar objects) may be massive enough to have formed by fragmentation rather than core accretion (Schlaufman 2018), and outside of the range in which planets can possibly grow to the isolation mass without consuming the majority of the disc mass (if the surface density follows equation 33 with Σ 0 ≈ 1800 g cm −2 ). While it may be possible to constrain planet masses from the geometry of the dust gaps (Rosotti et al. 2016), practically all gaps that have been resolved in discs are at 10 au from the host star. It is therefore not presently possible to constrain masses of forming planets in the inner disc. Recent developments in the processing of interferometric observations may in future be applied to resolve structures in discs down to sub-beam scales to test the hypothesis we have presented here (Jennings et al. 2020;Jennings et al. 2021).

Multi-planet systems in mean motion resonances
Mean motion resonances (MMRs) in planetary systems are thought to result from migration within a gas disc (e.g. Cresswell & Nelson 2006). Systems in MMR are unlikely to form as a result of rapid or high eccentricity migration (Quillen 2006;Pan & Schlichting 2017). Multi-planet systems rarely occupy resonant configurations (Fabrycky et al. 2014), which suggests that any disc migration does not ubiquitously result in resonant capture (see also : Hansen & Murray 2012;Petrovich et al. 2013;Chatterjee & Tan 2014). However, there exist a number of examples of planetary systems that are in resonant configurations, such as those around the M Dwarf GJ 876 (Batygin et al. 2015;Nelson et al. 2016;Cimerman et al. 2018), the solar mass Kepler-36 (Carter et al. 2012), and the very low-mass red dwarf TRAPPIST-1 (Gillon et al. 2017;Luger et al. 2017).
The existence of any resonant systems requires a degree of slow, low-eccentricity migration in some systems. In support of such migration being frequent, resonant chains have also been shown to be easily broken during type I migration (Goldreich & Schlichting 2014;Hands & Alexander 2018). Further, Longmore et al. (2021) suggest that resonant systems may be correlated with low density stellar environments, in which case MMRs resulting from isolated planet formation processes in a disc may be more frequent than previously thought. Hence, it is possible that (nearly) all planets experience slow, low eccentricity migration in the disc phase, with a subset escaping MMR configurations or undergoing subsequent dynamical perturbation. By contrast, Morrison et al. (2020) demonstrate planets can form in resonant chains during oligarchic growth in a depleted gas disc, without much migration. Either scenario might give rise to a fraction of systems in MMR, while not the majority. In either case, the upper limit we suggest in this work only requires that rapid type II migration does not result in massive planets on short period orbits, and does not preclude some degree of slow type I migration.

Circumbinary planets
Circumbinary planets around binaries with periods 10 days have been found to be abundant, with incidence rates comparable to planets around single stars (Armstrong et al. 2014, although one must also factor in the high probability of transit - Martin & Triaud 2015). There are now 14 known transiting circumbinary planets (Kostov et al. 2020), and most of them orbit very close to the innermost stable orbit around the binary. In situ formation is very unlikely (e.g., Paardekooper et al. 2012), and high-eccentricity migration is impossible (due to high-eccentricity orbits being unstable). The orbits of these planets are consistent with formation at larger radius followed by disc-driven migration, as predicted prior to their discovery (Pierens & Nelson 2008). However, the torques exerted on these planets mean that migration in such systems does not reflect the scenario in single systems, and would not necessarily contradict the hypothesised mass limit in this work.

Mechanisms operating on perturbed systems
We have illustrated that there may exist an upper limit for the growth of inner planets. This begs the question: what mechanisms could yield deviation from this mass in the observed exoplanet sample?

Binary/multiple interactions and stellar encounters
High eccentricity migration, the result of dynamical perturbation as opposed to disc migration, is a promising mechanism to produce hot Jupiters and other high eccentricity massive planets. Since a high fraction of stars are multiples (Duquennoy & Mayor 1991;Raghavan et al. 2010, although not necessarily most -Lada 2006, one might expect some fraction of systems to be influenced by a companion. In particular, the exchange of angular momentum during Kozai-Lidov (KL) oscillations can produce extreme eccentricities and inclinations (Nagasawa et al. 2008;Naoz et al. 2012;Naoz 2016). This mechanism may produce hot Jupiters when coupled with tidal circularisation (Section 5.3.2). Although the absence of high eccentricity transit detections in Kepler data marginally opposes this scenario (Dawson et al. 2015), this may be explained by non-steady state flux of planets inwards (early migration) or the considerable uncertainties in the circularisation time-scale.
Investigations into the occurrence rates of hot Jupiters around binary components are often sub-divided into short and long period binaries. Ngo et al. (2016) find that hot Jupiters are associated with double the binary fraction in the separation range 50−2000 au than field stars (see also Fontanive et al. 2019). The opposite trend for binaries with separations 1−50 au (suppressed by a factor ∼ 4 for hot Jupiter hosts) suggests a transition from formation and inward scattering to formation suppression. Moe & Kratter (2019) show that, accounting for selection biases, the correlation of hot Jupiter occurrence with wider binaries only apply to the most massive hot Jupiters that did not necessarily form by core accretion. However, this is not directly comparable to the similar result of Belokurov et al. (2020), who find a higher frequency of tight binaries (separations ∼ 0.1−10 au) among hot with respect to cold Jupiters with mass > 1 J (the majority with masses 1−4 J ). Whether the 4 J threshold, motivated by the metallicity occurrence rate correlation (Schlaufman 2018), is directly related to a correlation with binary occurrence remains an open question. It is possible that lower mass hot Jupiters form at smaller separations and require a tighter binary companion to perturb them, which might reconcile the findings of Belokurov et al. (2020) and Ngo et al. (2016). While the nature of these correlations requires further analysis, it is clear that some fraction of massive planets on short period orbits have experienced perturbation from a binary companion.
Apart from long time-scale oscillations, rapid planet-planet scattering has been proposed as another mechanism to produce high eccentricity (massive) planets (e.g. Chatterjee et al. 2008;Carrera et al. 2019). The close passage of two stars can also gravitationally perturb a (proto)planetary system, inducing instabilities and scattering bodies away from their initial orbital configuration (e.g. Malavolta et al. 2016;Pfalzner et al. 2018a;Cai et al. 2017-see Davies et al. 2014, for a review). Although such close interactions during the protoplanetary stage (within ∼ 3 Myr) occur mainly within small-scale bound stellar multiple systems (Bate 2018;Winter et al. 2018b), if the whole cluster remains bound over 100 Myr then planetary scattering by neighbouring stars is possible at stellar densities 10 2−3 stars/pc 3 (e.g. Malmberg et al. 2007Malmberg et al. , 2011Spurzem et al. 2009;Li et al. 2019;Stock et al. 2020;Wang et al. 2020a,b). Implementing N-body calculations of systems with multiple massive planets, Shara et al. (2016) find that hot Jupiter formation by scattering occurs for ∼ 1 percent of systems at stellar densities typical of known open clusters. Scattering may also be significantly enhanced when accounting for the increased crosssections of binaries (Li & Adams 2015;Li et al. 2020), possibly in conjunction with KL oscillations or von Zeipel-Kozai-Lidov if the mass of inner and outer companion are comparable (Naoz et al. 2011;Rodet et al. 2021). Hot Jupiter formation by encounters and subsequent scattering is directly supported by the enhanced occurrence rates in the dense open cluster M67 (Brucalassi et al. 2016) and the younger, less dense Praesepe (Quinn et al. 2012). Indirectly, multiple studies have found that hot Jupiters are more frequent at high position-velocity density (Hamer & Schlaufman 2019;Winter et al. 2020b;Adibekyan et al. 2021). This finding may be related to the slower dynamical heating of open clusters (Tarricq et al. 2021), although if stellar age estimates are sufficiently uncertain the findings could also be attributed to tidal inspiral of older hot Jupiters (Hamer & Schlaufman 2019). However, tidal inspiral would not explain the preference for planets in the LDE sample at or below the Kepler inferred growth limit among the non-hot Jupiter planets, discussed in Section 4 (see Figure 11). Stellar scattering preferentially influences planets at wide separations, as encounter rate (scattering cross-section) scales with pl for gravitationally focused encounters, or 2 pl for hyperbolic ones. Assuming the planet remains bound post-encounter, it most frequently loses angular momentum, resulting in decreased semi-major axis (Hills 1975;Hut 1983;Heggie & Rasio 1996). We therefore expect encounters to produce a population of planets with masses preferentially above any limit for formation in isolation. Our findings in Section 4 are broadly consistent with this expectation.

Tidal circularisation and atmospheric photoevaporation
Once a massive planet has lost angular momentum by a stellar encounter, the induced high eccentricity can result in numerous close passages with the host star. Further angular momentum loss, this time through tidal interactions with the host star, can then lead to orbital circularisation, inclination damping and shortening of the orbital period (Hut 1981;Rasio et al. 1996;Ivanov & Papaloizou 2004Faber et al. 2005;Jackson et al. 2008). Large quantities of energy can be lost during this tidal interaction, making it a promising mechanism for the formation of hot Jupiters. Once on short period orbits, massive planets may lose mass due to irradiation and photoevaporative heating by the central star (Owen & Jackson 2012). This photoevaporation has been suggested as the origin of the observed 'radius valley' (Owen & Wu 2017;Mordasini 2020).

External photoevaporation
Protoplanetary discs can be photoevaporated by massive neighbouring stars in young (few Myr old) star forming environments (Johnstone et al. 1998;Adams et al. 2004;Winter et al. 2018a). This is both predicted theoretically and found empirically to reduce the disc mass (Maucó et al. 2016;Ansdell et al. 2017;Eisner et al. 2018;Facchini et al. 2016;Haworth et al. 2018;Winter et al. 2019b;van Terwisga et al. 2019) and shorten the disc life-times (Stolte et al. 2004(Stolte et al. , 2010(Stolte et al. , 2015Fang et al. 2012;Guarcello et al. 2016;Concha-Ramírez et al. 2019;Winter et al. 2019aWinter et al. , 2020a. In terms of the planet masses, one would therefore expect to reduce masses of planets that form, by some combination of reducing the local disc surface density and limiting the time-scale for formation. While the principles of photoevaporation are well-known, and there is now ample direct evidence to suggest it operates on discs (e.g. O'Dell & Wen 1994;Kim et al. 2016;Haworth et al. 2020), how this may influence exoplanet architectures remains uncertain. A reasonable assumption is that photoevaporative winds should decrease the mass-budget available for planet formation, and do so preferentially at large separations from the host star (Johnstone et al. 1998). This need not necessarily suppress planet formation; photoevaporation of the gas in the disc may also result in enhanced dust-to-gas ratio, which can make the growth of planets by dust instabilities (e.g. streaming) more likely (Throop & Bally 2005). Nonetheless, if discs are externally depleted in dense environments, and the most massive planets must grow from gas in the outer disc, this might explain the dearth of massive planets and brown dwarfs in globular clusters (Gilliland et al. 2000;Bonnell et al. 2003).

Other mechanisms
Numerous other mechanisms may exist that drive deviation from the hypothesised mass limit. The deposition of short lived radionuclides into the discs due to supernovae or the winds of massive stars (e.g. Martin et al. 2009) can lead to the heating and dehydration of planets (Lichtenberg et al. 2019). Ionisation by cosmic rays that may also lead to more compact discs in some star forming regions (Kuffmeier et al. 2020). Enhanced metallicity in some environments may speed up planet formation and slow down disc dispersal (Ercolano & Clarke 2010). This list is unlikely to be exhaustive.

Ramifications for the Solar System
With the exception of Jupiter, the planets hosted by the Sun are much lower mass than the upper limit for which we have presented evidence in this work. This may be the consequence of premature disc dispersal. The idea that the Solar System has been sculpted by its star formation environment is not a new one (see e.g. Adams 2010, for a review). Evidence ranges from the the abundance of daughter isotopes of short lived radionuclides that may have originated in a supernova nearby the young solar nebula (e.g. Meyer & Clayton 2000;Goswami & Vanhala 2000;Adams & Laughlin 2001;Wadhwa et al. 2007;Gounelle et al. 2009;Banerjee et al. 2016;Boss 2017), to the high eccentricities and inclinations of the trans-Neptunian bodies that could originate from a stellar encounter in a moderate density cluster (Becker et al. 2018;Pfalzner et al. 2018b;Moore et al. 2020;Batygin et al. 2020). By the definition of Winter et al. (2020b) the Sun occupies a high phase space density. It is possible that the Sun originated in M67, members of which have similar ages and metallicities (Yadav et al. 2008;Heiter et al. 2014). If such a cluster initially had ∼ 10 4 members (Hurley et al. 2005), this would be consistent with constraints from the Solar System properties (Portegies Zwart 2009;Adams 2010;Pfalzner 2013;Moore et al. 2020). However, this scenario remains a subject of debate (Jørgensen & Church 2020;Webb et al. 2020).
Since the giant planets in the Solar System do not have high eccentricities, they are unlikely to have been significantly dynamically perturbed. However, if the young Solar System was exposed to strong UV fields in its birth environment, this may have reduced the gas reservoir from which planets could accrete before they were able to grow to the isolation mass. In this case, the feeding zones could remain small such that the young planets would not merge, and we may expect much lower planet masses than the limit we infer in this work. This is a speculative scenario, but a feasible way to produce a system of low eccentricity planets well below any hypothesised mass limit in the growth of planets forming in isolation.

CONCLUSIONS
In this work, we investigate the distribution of exoplanet architectures in light of evidence that some subset of planetary systems have been perturbed from their formation configuration. In particular, we consider whether there may exist an upper limit in exoplanet masses that is disguised by the perturbed planets.
For close-in planets ( 1 au), we present an approach for modelling the period-radius distribution of Kepler multiples, accounting for the two-dimensionality of the data. Around solar mass hosts, a model incorporating an upper limit in the planet masses of the form: can reproduce the observed distribution of periods and radii when coupled with an empirical mass-radius relation. The power-law index = 1.5 is consistent with the data, as would be expected for in situ growth limited by the isolation mass. Imposing this upper limit provides a better fit to the data than a separable PDF in period and radius.
For longer period RV detected planets with pl > 0.2 au and pl sin > 300 ⊕ , those hosted by stars at low position-velocity density have masses that are preferentially below the Kepler inferred limit with = 1.5. This finding cannot originate from tidal inspiral of hot Jupiters ( pl < 0.2 au) and is not easily explained by agerelated observational biases. This suggests that planetary systems around high density host stars are more likely to have been scattered inwards, consistent with the suggestions of higher planet multiplicity among host stars in low density environments (Dai et al. 2021;Longmore et al. 2021).
We conclude that, given the mounting evidence that a subset of planetary systems have undergone dynamical perturbation, attempting to identify planets that have not been perturbed can offer insight into the formation process.

SOFTWARE
We thank the contributors for making the following packages, of which we made use, public: P , M (Hunter 2007), N (Harris et al. 2020

A3 Fixed inner-radius (F -_ )
The posterior as a result of fixing the inner-radius is shown in Figure A3, with the same priors as described above (Section A1), plus a uniform prior for between zero and , where ≡ arctan for power-law . In this case large (steep power-law) is associated with large . We also find that the steepest power-laws ( → /2) are effectively forbidden by the upper limit we place on 0 . To yield better constraints on we require additional data or a more restrictive and physically motivated prior.

A4 Double broken power-law model (BPL-PR)
The posterior distribution for the fitting of the separable model in period and radius is shown in Figure A4. The priors are uniform in the range break = 1−100 days, break = 1−10 ⊕ and power-law indices −10 to 10. The resulting distributions are similar to those inferred by M18. This paper has been typeset from a T E X/L A T E X file prepared by the author.  Figure A1 except for a fixed inner edge (log in = −1.4) and variable arctan (F -_ ).