Low strength and a shared size distribution connect the hot and cold Kuiper belt populations

Domain: earth_science
Published: 2025-11-28
Editorial outcome: major_revision

Abstract
Understanding how the observed size distribution of Kuiper Belt objects emerged from the primordial planetesimal population remains a major challenge in planetary science. In particular, the similarity between the size distributions of dynamically hot and cold populations in the 100–400 km range calls for a unified explanation. Motivated by this observational constraint and by theoretical models of pebble accretion, we model the coupled dynamical and collisional evolution of the trans-planetary disc under the early instability scenario. Here we show that a common planetesimal formation mechanism operating up to about 400 km and modulated by an exponential taper shape, combined with a transition to a shallow power law distribution in the inner disc, can reproduce the current size-frequency distributions of cold and hot populations. This outcome requires assuming that planetesimals were weaker to impacts than previously thought. Our results connect formation models with current observations, thereby constraining the early conditions of the Solar System and the properties of planetesimals.

Introduction
The inventory of dynamically hot and cold Kuiper belt objects (KBOs) has been proposed to be nearly complete for objects with absolute magnitude H r < 5.3, and this has been confirmed by dedicated surveys. The achievement of near completeness in that magnitude range opens new possibilities to better constrain the early Solar System after the gas dispersal phase.

The Outer Solar System Origins Survey (OSSOS), along with following extensions referred to as OSSOS++, provides a comprehensive debiased sample of dynamically hot and cold Kuiper belt (KB) populations (hereinafter HKB and CKB, respectively), together with their absolute magnitude H r distribution. Supported by such surveys, the H r distribution of the cold classical population was found to be accurately approximated by an exponentially tapered function in H r ≃ 5 to 12 range. Later findings indicate that the absolute magnitude distributions of dynamically hot and cold populations in the 5.5 < H r < 8.3 range (corresponding roughly to 100–400 km diameter range) are essentially indistinguishable. Such observational evidence provided by OSSOS++ contradicts the previously prevailing view, when it was widely accepted that the cumulative size-frequency distributions (SFD) of the hot and cold components of the KB were well described by broken power laws of the form N(>D) ∝ Dq, where N(>D) is the number of objects with diameter larger than D, and q is the cumulative slope of the distribution. Specifically, at the large-size end (D > 100 km), the cold component was thought to have a steeper slope (q cold = −7.5) than the hot component (q hot = −4.5). The evidence provided by has posed a challenge to the previous notion that the primordial Kuiper belt displayed different SFDs slopes in its inner and outer regions, implying different formation processes.

The general understanding of the outer solar system dynamical evolution is encompassed under the umbrella of the new Nice model. According to this model, the planetesimals surrounding the giant planets, remaining after their formation, were located in a disc extending about 22 to 30 au, with an estimated mass about 30 M⊕. As dynamical instability ensued, the influence of Neptune played a significant role in depleting that disc. Most of the mass was dispersed at very large heliocentric distance, while some bodies were implanted in the current region of Plutinos, and others eventually found their way to the HKB.

There is a consensus on the hypothesis that the cold population, with a current mass of about 10−3 M⊕, was born in situ and remained primordial, based on the large fraction of binaries still in the CKB and their distinctive red colours. This argument was later reinforced by the New Horizons images of Arrokoth, showing scarcity of craters on its surface, suggesting that collisional evolution did not significantly alter the primordial SFD of that population. Based on measurements by OSSOS++, it has been suggested that the SFD of the cold and hot components may have maintained their primordial shape for objects in the size range of about 100–400 km. Accordingly, the same mechanism responsible for the accretion of bodies in their primordial formation region might be common to both populations.

As far as planetesimal growth is concerned, the streaming instability (SI) model gives valuable insights into the initial phase of planetesimal formation by providing an efficient and fast mechanism for growing planetesimals from dust particles or pebble clumps within protoplanetary discs. The SI is an aerodynamic process driven by the velocity difference between solid particles and gas and the reaction of solids onto the gas. This mechanism led to the formation of high-density regions, providing a favourable environment for gravitational collapse and subsequent aggregation into planetesimals.

Lyra et al. found that considering a pebble size distribution in the accretion process (also refereed as polydisperse pebble accretion) following SI was critical in planetesimal formation. According to this study, such a process was responsible for accretion of up to a few hundred km (about 300–400 km) at heliocentric distances as far as the current location of the CKB. This finding is essential in defining a common formation mechanism for the primordial trans–planetary populations (between 22 and 30 au) and the cold population of bodies up to such size, confirming previous conjectures. According to accretion models, further growth strongly depends on mass conditions at different heliocentric distances. This is the reason why accretion to Pluto (and possibly larger) size bodies was efficient in the inner region between 22 and 30 au. Instead, beyond 40 au, planetesimals should have remained as such, as further accretion was prevented by low mass density conditions in the nebula.

The inventory of large bodies in the Kuiper belt seems nowadays relatively complete for H r < 5.3, while the SFD of fainter (smaller) bodies remains largely incomplete. Valuable insights into impactor SFD and collisional processes are achieved by exploring planetary surface impact craters. Thus, the surfaces of Pluto and Charon offer a unique opportunity to examine the crater population in the outer solar system. Robins et al. provide results that reconcile earlier discrepancies in the estimated crater size distribution slopes obtained from the New Horizons fly-by of both worlds. They estimate that the size distribution of small craters has a shallower slope than for large ones, consistent with previously reported results by Singer et al. (hereinafter S19). Nevertheless, according to the authors, the crater size-frequency distribution slopes are uncertain, in particular at small sizes. In fact, slope values for craters below 10–15 km is −0.7 ± 0.2, and −2.8 ± 0.6 for larger craters; implying a deficit of small KBOs below 1–2 km. Such estimated slopes are very sensitive to the choice of parameters and the craters used in the statistical analysis, and they must be taken as minima in the uncertainty determination. On the other hand, Morbidelli et al. (hereinafter M21) established a connection between their re-examined Pluto and Charon cratering record and Arrokoth’s cratering record, including the dust measurement (down to 2 μm approximately) obtained from the New Horizons space mission by NASA. Their findings revealed that the SFD shows a moderate decline below 2 km, which contrasts with the claim made by S19 of an exceedingly shallow slope. In addition, M21 argue that the shallow portion of the SFD extends to around 30 m before transitioning to a steeper slope to match the dust abundance constraint.

Collisional models are commonly used to study small body populations evolution in the solar system. However, since the adoption of the new Nice model framework, only two studies have examined the collisional evolution of the trans-planetary disc in detail. First, Asteroid-LIke Collisional And Dynamical Evolution Package (ALICANDEP-22) was applied to conduct a systematic study exploring different initial distributions for the trans–planetary disc and different scaling laws that govern the impact response. They got a good match with available observational evidence when it was still thought that the CKB and HKB had different SFD slopes at the bright end. Later, the Boulder code was applied to study the collisional evolution of the primordial KB and other related populations as the destabilised population and the Jupiter Trojans. The main constraints used in the latter include the size distributions of the Trojan asteroids and craters on the giant planets satellites.

The presence of Pluto-sized and larger bodies in the hot population implies that the initial SFD of the inner trans-planetary disc must have extended beyond 400 km. Moreover, theoretical studies of planetesimal accretion suggest that, in the inner disc, a power law tail likely developed with a relatively shallow slope, enabling the growth of several thousand dwarf planets. When combined with an exponentially tapered distribution, characteristic of the cold population, this configuration provides a natural explanation for the properties of the hot and cold Kuiper Belt populations within a unified formation framework.

Here we show that the observed SFDs of both dynamically hot and cold Kuiper Belt populations can be reproduced by modelling their collisional and dynamical evolution from a common primordial distribution. We adopt an initial SFD shaped by a shared formation mechanism up to 400 km, transitioning to a shallow power law tail in the inner disc, and evolve it under the early instability scenario. Our simulations indicate that matching current observations requires assuming low material strength for planetesimals.

Results
Outline of the model
We modelled the KB collisional and dynamical evolution using the ALICANDEP-22 package. This model simulates collisional evolution within the framework of the current interpretation of the Nice model for the dynamical evolution of the outer solar system. Dynamical evolution is carried out through statistical elimination of bodies from an inner disc and implantation into the current KB. The modelled scenario begins with the initial mass concentrated in an inner disc (outside Neptune’s orbit), ranging from 22 to 30 au. We also considered a small concentration of mass between 42 and 47 au, standing for the in-situ formation of the cold classical population. The initial population born in the inner disc undergoes three evolution phases, including collisional and dynamical evolution processes such as excitation, migration, and mass depletion. Each phase consists of the following features. Phase 1: We call this the pre-instability phase, lasting from time t 0 to t i. It denotes the onset of Neptune’s outward migration as it enters the inner disc. This phase results in an increase in eccentricity and inclination of planetesimals. Such dynamical excitation increased relative velocities (see “Methods”) and, therefore, the overall collisional evolution. Additionally, some smooth mass-loss during this period—according to dynamical models—is also included. Phase 2: This period is known as the instability phase. In our model, it lasts from t i to t i + Δt. In this phase, Neptune passed through the inner disc, which led to minor collisional evolution and significant dynamical depletion. During this time, most planetesimals were ejected from the solar system; some were removed inwards, while others were implanted into the excited Kuiper belt populations. Phase 3: This is the post-instability phase, lasting from t i + Δt to the end of the simulations (4500 Myr). At the beginning of this phase, the inner disc was already depleted, and the outer regions corresponding to Plutinos and hot classical belt populations were populated. Negligible collisional evolution took place in this phase.

In ALICANDEP-22, the outer region of the solar system is divided into three zones, representing different Kuiper belt populations: Plutinos, cold classicals, and hot classicals. The collisional evolution of each zone is performed taking into account mutual interaction between different zones. The orbital elements, semi-major axis (a), eccentricity (e) and inclination (i) of the inner zone and the three outer zones for each dynamical phase, are shown in Table 1. For the instability onset we tested t i = 10 to 100 Myr, with Δt = 10 and 30 Myr.

Regarding the response of bodies to impacts, numerical studies tend to shift from 200 m to around 20 m, the size at which scaling laws for the threshold specific energy for disruption (Q* D) have a minimum. This was first suggested to suitably fit the cratering record on Pluto, Charon, and Arrokoth, and to explain the survival of the satellite of the Trojan asteroid Eurybates (Queta). Moreover, a study on collisional evolution of the primordial Kuiper belt, its destabilised population, and the Trojans asteroids, also finds that such scaling laws produce better results to match the Trojan asteroids populations and the crater distributions on giant planets satellites. Further details about the implementation of scaling laws are presented in the Methods section.

The initial SFD is an essential boundary condition in the progress of populations of systems under moderate collisional evolution. In such systems, as in the case of the Kuiper belt, the SFD slopes for large bodies are not modified along the evolution; thus, their features must be primordial. Accordingly, we must reproduce the gradual transition towards larger bodies currently observed in the SFD of the implanted hot component. Therefore, we propose that the initial population in the inner trans–planetary region combines a tapered distribution up to D tr about 100–400 km, with a decreasing power law distribution for D > D tr. Such population contained around of 4000 Pluto–size (and larger) bodies. Two values were tested for D tr: 100 and 400 km. In the case of the primordial cold disc, an SFD identical in shape to the inner disc but with much less mass—and no population of bodies larger than 400 km—was chosen as the initial distribution, following the planetesimal formation mechanisms discussed above. At the small-size end, the cumulative H r magnitude distribution is governed by an exponent α s (see “Methods,” subsection Initial size distribution). In the following section, we show results for primordial distributions with α s values of 0.4 and 0.3, which are our steep and shallow distributions (see “Methods,” subsection Model flexibility to match observational uncertainties, for a detailed exploration of the initial SFD parameters). The total mass of the primordial population is approximately 30 M⊕, assuming an average bulk density of 1000 kg/m3.

Observational constraints
The main constraint in this study comes from the debiased magnitude distribution of the hot and cold populations provided by updated surveys. The OSSOS++ team graciously provided us with the magnitude distribution of each component, making this task easier. Due to the fact that the OSSOS++ survey may overestimate the size distribution of bright bodies, we preferred to use the Minor Planet Center (MPC) database for the distribution of bodies with H r ≲ 5.5. Several pencil-beam surveys have been conducted to identify faint Kuiper belt objects. Bernstein et al. used the Hubble Space Telescope (HST) to measure the faintest Kuiper belt objects for which a direct SFD can be obtained. They detected three CKB bodies and none in the HKB population. Kavelaars et al. debiased the three detections in the cold population using the same procedure as for OSSOS++. This estimation is also used in this work as an observational constraint. In the Methods section, we utilize observational uncertainties to investigate the flexibility of the model while ensuring it remains consistent with the range of possible SFD permitted by the OSSOS++ data, including the debiased estimation by Bernstein et al.

We assumed the average albedo of each population (p r = 0.15 and 0.08 for cold and hot KBOs, respectively) for size conversion. Both tested transition sizes D tr = 100 and 400 km provide satisfactory results; however, we centred our analysis on the case corresponding to D tr = 100 km, because it best matched the constraint of a smooth slope transition in the 100–400 km size range to larger bodies. In the subsection Considerations on the Kuiper Belt Populations, two additional constraints are discussed: the change in the SFD slope at 1–2 km and the Jupiter trojans SFD slope for 5–100 km size bodies.

Model outcome
Figure 1 shows the model outcome of the cumulative SFD for the hot and cold populations after 4500 Myr of evolution. For comparison, the SFD corresponding to the debiased absolute magnitude distribution for each dynamical population is also shown along with the corresponding Poisson 95% confidence range estimated by the OSSOS++. Different scaling laws were tested, showing that the final distributions, for large bodies, are not very sensitive to the very values of the scaling law for shattering (Q* S), provided that we choose a weak one and that the minimum value of the distribution is about 20 m (see “Methods,” subsection Impact outcome).

Hot classical population
The match between the modelled hot classical population and the debiased distribution is excellent, particularly in the size range from about 70 to 400 km. For larger bodies, the modelled SFD follows the MPC database, essentially complete at this size range. Nevertheless, for completeness, we also explore the uncertainty range of the OSSOS++ survey (see Methods, subsection Model flexibility to match observational uncertainties).

The hot classical population was born in a massive disc, though with low relative velocities due to small values of eccentricity and inclination. That circumstance only allowed fair collisional grinding down to about 20 km. However, at sizes below 100 km, a moderate wavy behaviour in the SFD can be noticed, in the case where the weakest scaling law is assumed (solid blue line in Fig. 1a). Collisional grinding appears then subtle, and the very characteristics of each fragmentation scaling law leave minor fingerprints in the evolved SFD, even when delaying the onset of instability from 10 to 100 Myr (see Methods, subsection Time dependence). The shift of the minima in the scaling laws down to some 20 m is crucial in the propagation of the wavy pattern, even at large sizes (about 10 km).

Cold classical population
Agreement of our results with the debiased distribution is very accurate for the cold component across the comparable size range (70 < D < 400 km). That is particularly evident when the assumed initial distribution has a steep slope at the small–size end (α s = 0.4), as is observed today. Further simulations constrain the possible values of α s in the range of 0.27 to 0.42 to be consistent with the error bar of the debiased estimation of the faintest detected objects in the CKB (black error bar in Fig. 1b). The uncertainties of the OSSOS++ also allow a modest variation in the slope of the taper function at large sizes. The threshold values defining the initial SFD are explored in the Methods section.

Negligible collisional evolution took place in this region, and the final distribution is a slightly depleted version of the primordial one. ALICANDEP-22 can keep track of the amount of objects that are disrupted, or just shattered and gravitationally re-aggregated. In this way, we find that—in the cold classical population—bodies larger than about 1–10 km are expected to be primordial and to have experienced only cratering events, even in the case the weakest shattering scaling law is considered in our model. Consequently, their dynamically relaxed initial conditions strongly support that this population was preserved in its primordial state due to negligible collisional evolution. This confirms findings indicating that Arrokoth is most likely a primordial body.

Considerations on the Kuiper Belt populations
It has been pointed out in the literature that the KBO SFD exhibits distinctive features. First, a break at 100 km referred to as the “knee” was found, where the cumulative slope transitions from steeper to shallower. Also, a secondary break at about D = 1–2 km, known as the “elbow,” was suggested, where the slope shifts from a value in the range of −2.8 to −3.2 to an even shallower slope with a corresponding q in the range of −0.7 to −1.2, in terms of cumulative distribution (see Introduction). We find that an exponential taper shape SFD of the primordial population, below 100–400 km, easily reproduces the knee at about 100 km, as shown in Fig. 1, for the hot and cold dynamical population, and in Fig. 2 for the total KB population (Plutinos, HKB, and CKB).

The transition region between the knee and the elbow (2 < D < 100 km) must follow a cumulative SFD slope similar to the one observed in the captured population of Jupiter Trojans (q ~ −2), described by the red dotted line in Fig. 2. We faithfully reproduce such behaviour in all the simulations adopting a steep slope (α s = 0.4) for the current H magnitude distribution. In contrast, beginning the evolution with a shallower slope (α s = 0.3), implies that the final SFD will also be shallower than expected in that transition size range. This is due to the fact that collisional evolution is unable to increase the slope to the one expected for Jupiter Trojans. We show this in Fig. 2, highlighting the need for the chosen slope (α s = 0.4) to meet observational constraints.

The behaviour of the SFD in the vicinity of the elbow deserves closer examination. Following previous studies on crater counts on Pluto and Charon, we consider the impact record estimated for Vulcan Planitia on Charon in S19, as it may best represent the primary population of impactors. Therefore, we obtained the nominal impactor SFD and plotted it along with the outcome of our model (Fig. 2) by suitably scaling it to the Kuiper belt distribution.

The model that best agrees with the predicted impactor population of cratering record estimations features a steep initial slope (α s = 0.4) and the weakest scaling law. This is represented by the solid blue line in Fig. 2. The match is remarkable in the 0.7–10 km size range, approximately. The S19 impactor distribution at smaller sizes instead bends over along a shallower slope. There is limited data availability in the sub-kilometer size range, so uncertainties related to crater count are relevant. Figure 2b shows the elbow size range in detail, including uncertainty estimation within 2σ (for small bodies), using as a reference the number of impactors corresponding to about D = 1.5 km, where the distribution gradually changes its slope. Our best matching case reproduces fairly well the gradual slope change (between 1 and 2 km) confirmed by crater count, while for smaller bodies it is close to the upper 2σ slope uncertainty, for the estimated distribution of impactors. For completeness, the cumulative distribution estimated by M21 is also plotted in Fig. 2b. The slope of the SFD corresponding to M21 model is similar to that found in our study in the 0.1–1 km size range, except for their underestimated population (by a factor about 3) below 1 km. Given the uncertainties in slope estimation of the crater size distribution, mainly due to the smooth transition from large to small craters, S19 and M21 estimations are not necessarily at odds with each other, and our model seem to be compatible with both interpretations within uncertainty limits.

In any case, we find that the only way to reproduce the elbow is to assume a weak strength for bodies in the primordial trans–planetary disc. This finding is in agreement with the 2023 study by Bottke et al. of collisional evolution based on current Trojan asteroids SFD and giant planets satellites cratering data, requiring a similar scaling law to be matched. Thus, despite the relatively small influence of collisional evolution, restricted to an early phase and constrained by the onset of the early solar system dynamical depletion, it is still significant in producing a bump at 1–2 km and propagating the wavy shape up to 10 km in the evolved SFD. Observational efforts are imperative for improving our comprehension of sub-kilometre outer solar system bodies and refining our understanding of the distinctive features characterising such population.

While the SFD estimated by M21 broadly corresponds to our overall final SFD in the approximate size range of 2–100 km, differences show up at smaller sizes, particularly around 30 m, where the slope turns steep in M21. In our model, the transition to a steep slope occurs at about 100 m, rather than at 30 m. Therefore, M21 underestimate the total number of bodies in the 10−2 < D < 1 km range, compared to our results. At the smallest size end (D < 10−3 km), both models tend to converge to the steady state slope. An explanation for such discrepancy may be that their estimation of the total SFD is derived from the impact cratering record of bodies associated with different dynamical populations, such as Pluto and Charon (in the plutinos 3:2 mean motion resonance region), and Arrokoth (CKB). The key difference in our model, compared to M21’s, is the ability to follow the collisional evolution of individual populations, including their mutual interactions, resulting in accurate outcome of the end–state overall populations.

Methods
The collisional model used in this work consists of two main components. One takes care of determining the impact outcome of any collision, the other one handles the collisional evolution, using the output of the first component as input to iterate the evolution of a given initial population over 4500 Myr. The entire code package is named ALICANDEP-22 and was fully described and tested in. We provide here some details that may help understanding the model outcomes described in this work.

Impact outcome
To model the impact outcome of any collision, we use the fragmentation and re-accumulation model by Petit & Farinella, with further improvements. This part of the package computes the number of fragments produced in any possible collision between objects of different size (mass) bins, considering a given scaling law for shattering (Q* S). Q* S is the threshold specific energy required to break a given mass body into a spectrum of single fragments, the largest one having exactly half the mass of the original target. Instead, Q* D is the specific energy required to disperse half of the mass of a given target. This is usually estimated from the outcome of smoothed-particle hydrodynamics codes. It is important to note that Q* D and Q* S are equal in the strength regime (small bodies), where no self–compression nor fragment re-accumulation takes place, but they are different in the gravity regime (large bodies). Studies suggest that KBOs and Jupiter Trojans—likely a captured population of the initial inner disc—may be very weak, with a minimum of the Q* D scaling law at 20 m. The fragmentational algorithm is based on the specific shattering energy (Q* S). This allows the estimation of the fraction of re-accumulated mass by self—gravity, a distinctive feature of the model. We reformulated the usual Q* D expression following de Elia and Brunini as:

Q* D = C 1 D^λ 1 [1 + C 2 D^λ 2] (1)

where C 1, C 2, λ 1, λ 2 correspond to the coefficients. Using the values expected for scaling law for icy bodies corresponding to an average impact speed of 3 km/s.

Then, from the parametric expression eq. (1), corresponding expression for Q* S is obtained.

Q* S = C 1 D^λ 1 [1 + C 2 D^λ 2] (2)

In this expression, changing the values of C 2 and λ 2 the corresponding Q* S can be changed. C 2 and λ 2 are parameters shaping the scaling law in the gravity regime. We use C 1 = 1.55, λ 1 = 0.36, C 2 = 8, and λ 2 = 1.8, to get our nominal Q* S scaling law, in agreement with Marschall et al. Thus, Q* S and Q* D overlap in the strength regime, and differ as expected in the gravity regime (see black and blue dot-dashed lines in (Fig. 3a). In our simulations, we also test two alternative scaling laws, obtained by scaling the nominal one up and down by a factor of 10 (Fig. 3a), to check the sensitivity of collisional evolution on the adopted scaling law.

Initial size-frequency distribution
Our study aims to build a common primordial SFD for the cold and inner discs for bodies up to a given size. The following methodology was used to build up such initial population. Assuming that the cold classical population retains its primordial SFD, as supported by the estimation of the magnitude distribution (N(<H r)), by Kavelaars et al., we build up the cold population starting from the following relationship:

N(<H r) = N c1 10^α s (H r − H 0) exp [−10^β(H r − H B)] (3)

where N(<H r) is the population of bodies with absolute magnitude less than H r. Parameter values are: H 0 = −2.6, α s = 0.4, β = 0.25 and H B = 8.1, from Petit et al., and N c1 is a constant value. To estimate the mass, we assumed an average albedo, p r = 0.15, typical of the cold population, and bulk density ρ = 1000 kg m−3, resulting in a total initial mass M 0 cold = 1.3 × 10−3 M⊕. In this way, we also obtain the primordial SFD for that region with one largest body of about 375 km in size.

To build up the distribution of the inner population, we use the same functional form as for the cold SFD, but scaled up to obtain a suitable larger mass by several orders of magnitude. Such function extends the distribution up to about D = 800 km size bodies, due to the first term in Eq. (4), supported by streaming instability simulations. Such distribution need to be combined with a population of even larger bodies so that the SFD must be extended to larger sizes (dwarf planets and beyond) in the inner disc. The estimation of the population of large bodies is the topic of ongoing research in planet formation. In fact, the evolution of the streaming instability power law initial mass function (IMF) is yet under study. However, some research on pebble and planetesimal accretion suggests that the IMF becomes shallower for larger bodies due to accretion. A slope of α l = 0.14 was estimated for such distribution, implying an upper limit of about 1000 bodies with H r < 2 (dwarf planet size bodies). Nevertheless, such estimation turns out to be clearly smaller than the requirement of up to 4000 large bodies necessary for the grainy Neptune migration mechanism. Here, we adopt a slightly steeper slope (α l = 0.18) at large size range so that grainy migration is preserved and dynamical depletion allows for the right number of dwarf planets currently observed in the Kuiper belt. Therefore, the cumulative SFD for the inner disc is obtained as the combination of an exponential taper function and a power law for bodies larger than D tr. The resulting function is as follows:

N(<H r) = N i1 10^α s (H r − H 0) exp [−10^β(H r − H B)] + N i2 10^α l (H r − H 0) (4)

In order to obtain the SFD in the hot population, we assume an albedo of p r = 0.08, representative of such population. Figure 3b shows the initial SFDs tested in this study for the cold population and the inner region for D tr = 100 km. Simply grafting the SFDs for bodies larger and smaller than D tr does not provide a suitable match to the observed populations when evolved. Instead, summing up the taper and the power law (Fig. 3b), leads to a distribution that provides an excellent agreement between the observed and modelled populations. We test D tr values of 100 and 400 km to look for the best fit to observables. We also test a shallower SFD initial slope (α s = 0.3) for small-size bodies, but in this case, collisional evolution is not able to rebuild the population to the currently estimated SFD.

Impact velocity and intrinsic probability
Abedin et al. estimated that impacts between KBOs can currently happen over a wide range of relative velocities, mostly well below 5 km/s, in agreement with previous estimates. Therefore, we assume a Maxwellian relative velocity distribution around the most probable value (v p) for each modelled population. v p is adjusted to the dynamical conditions governed by eccentricity and inclination at each time step for each population. Collisional probabilities are calculated accordingly, taking into account the corresponding intrinsic collisional probability P int. The intrinsic collisional probability of bodies occupying a given volume depends on the average relative velocity and the volume itself. The average relative velocity can be calculated as a function of the average orbital eccentricity and inclination of the system (v rel / sqrt(e^2 + sin^2 i)), while the volume depends on the range of orbital inclinations and the boundary orbital semi-major axes (see Table 1).

Figure 4 shows the change in v p and P int for the inner disc along phases 1 and 2, corresponding to t i = 10, 30, 50, and 100 Myr. In the outer disc, we assume that once bodies are implanted into the Plutinos and HKB regions, they do not change their orbit elements any longer. Thus, the most probable values of the impact velocity and the mean intrinsic probability for that zone are constant. The cold population increases its inclination during phase 1, implying an increase in v p and decrease in P int. The values of v p and P int for each KB population are listed in Table 2.

Time dependence
Dynamical models give a rough estimate of when the dynamical instability started (t i). In previous studies the dependence of the evolution on the onset of dynamical instability was searched for, reaching up to 700 Myr. In this new approach, based on an early onset, we test t i = 10, 20, 30, 50, and 100 Myr. The nominal duration of the instability time interval, Δt, is set at 10 Myr as a nominal value, but we also explore longer instability time spans up to Δt = 30 Myr, without significant changes in the final distribution. For all runs in this section, we consider the weakest scaling law.

Collisional evolution is limited for two main reasons in this system. On the one hand, as described by dynamical models, the mass of the inner belt has to be reduced to 20 M⊕ by the onset of dynamical instability. On the other hand, the characteristic initial SFD concentrates most of the mass in large bodies (D > 100 km), lowering collisional probabilities.

We explore the effect on the SFD of assuming different values for such onset time, and we find that delaying the onset of instability from 10 to 100 Myr has the only effect of slightly increasing collisional evolution, causing a dip around 10 km (from 3 to 30 km, approximately). In this case, the wavy pattern triggered at small sizes has enough time to extend to larger sizes (Fig. 5). However, no reliable observational evidence for the existence of such dip is currently available, but cannot be completely ruled out, as it lies below observational completeness (about 100 km).

Figure 6 shows changes in mass over time, due to collisional and dynamical effects, giving a snapshot of the outcome associated with varying the onset of the dynamic instability. Our results show that the later the dynamical instability sets in, the more significant the mass fraction lost by collisional grinding. For instance, if instability occurs early, at 10 Myr, only 2% of the primordial mass is lost due to collisions, but if it is delayed until 100 Myr, mass loss may be as high as 15%.

Model flexibility to match observational uncertainties
In this section, we explore the threshold values of the parameters that define the initial distribution, producing an evolved population consistent with observational uncertainties. It is worth mentioning that, for completeness, we also explore the uncertainty range in the SFD at large sizes provided by the OSSOS++ project. Equation (3) defines the initial cold population, characterised by the parameters N c1, α s, and β. Following eq. (4), the parameters for the inner population are N i1, N i2, α s, β, and α l. Our methodology for constructing the initial distribution begins with the cold population. The α s parameter was constrained using the error band for the debiased estimation for the three detections by Bernstein et al. in the cold population (approximately at 15 km). The β parameter was constrained using OSSOS++ uncertainties at large sizes, where these are more significant. Eventually, the normalization parameter N c1 needed adjustment. Then, we build the initial distribution for the inner zone population. Here, the premise is to keep the slopes of the exponential taper component (α s, β) constant, established in the previous step, leaving N i1 and the parameters defining the power law component (α l and N i2) as free parameters. N i1 and N i2 are fixed to ensure the primordial distribution starts with about 30 M⊕, keeping the number of dwarf planets consistent with the grainy migration scenario. After systematically exploring each parameter defining the initial SFD of the cold and inner populations, we identify the values listed in Table 3. In order to limit the number of simulations in this section, the scaling law is set as the one corresponding to the most fragile objects, and t i = 30 Myr.

Figure 7a illustrates the initial SFD in the inner belt and the evolved SFD in the hot population. The shallow beta (0.23) requires a steeper slope (α l = 0.45) of the power law contribution for the evolved HKB distribution to match the upper OSSOS++ uncertainty at larger sizes. Contrarily, steeper β requires a shallower power law slope to match the lower OSSOS++ uncertainty. Figure 7b shows the initial and evolved SFD for the cold population. A slope for small bodies as shallow as α s < 0.27 (0.42) can accommodate the uncertainties for the faintest detections in the cold population. OSSOS++ uncertainties at larger sizes require 0.23 < β < 0.3. Thus, a narrow range for each parameter is allowable to commit the observational uncertainties for the CKB. For reference, the best fit is achieved for α s = 0.4 and β = 0.3, closer to the upper threshold values.

For all the runs listed in Table 3, primordial mass in Pluto-class objects (D > 1800 km) is between 18 and 40% of the initial mass contained between about 1900 and 6000 bodies, in agreement with the grainy migration requirement. Furthermore, mass in bodies larger than 100 km is in the range of about 47–88%. According to our best fit, after 4500 Myr of dynamical and collisional evolution, the Poisson probability of having at least one Pluto-class object in the plutinos population is 0.73. The probability of having at least two dwarf planets (D > 1400 km) in the HKB is 0.80.

A deeper statistical analysis, such as a Bayesian MCMC approach, could be useful in the context of a more detailed investigation, which lies beyond the scope of this work. This method would allow for a probabilistic exploration of the parameter space, making it particularly valuable for population matching within specific size ranges, where more precise estimations of the number of bodies could be directly compared with modelling results.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Benavidez et al. concerns the size distribution of planetesimals in the outer Solar System. The authors assume a primordial size distribution based on the streaming instability for the cold classicals but include an additional power law for large planetesimals in the inner regions. The authors then run a dynamical model for the planets and planetesimals of the Solar System and show that they match well the current populations in the hot and cold components of the trans-Neptunian populations.\nI find that the result presented here is potentially important, since it addresses the connection between the primordial planetesimal sizes and the observed sizes in the outer Solar System. But the narrative and the figures are currently presented in a quite confusing way that makes the manuscript very hard to penetrate. I have some requests here for significant rewrites of the paper for more clarity. I encourage the authors to restructure their manuscript substantially in order to elucidate the importance of their result.\n(1) The details of the model here are very important for judging the veracity of the results and thus can not only be hidden away in the methods section. The main manuscript should contain a description of the dynamical model, what physics is included, what is the time-scale of the dynamical instability, what is the duration of the simulations (4.6 Gyr?) and how the results depend on some of these choices.\n(2) The authors quote Morbidelli et al. (2021) as a key comparison paper. However, how their results differ or compare with M21 is in the end stated only in quite terms. I encourage the authors to make their differential analysis with M21 much more clear. Do these results here support M21 or are they more aligned with Singer et al. (2019)? This is an important debate in the community that this paper can contribute to.\n(3) The caption of Figure 1 is way too short. Please expand here on all the simulation parameters indicated in the legend (QS*, SL1, alpha). What is the time when this is plotted? How does the cumulative number evolve with time? What is the role of collisional grinding relative to dynamical depletion? Also, the authors should demonstrate here or in another plot what happens for different instability time-scales. Can the authors exclude very long instability time-scales from this analysis?\n(4) The caption of Figure 2 is also way too short. What do we learn from the discrepancy between the red dotted line and M21? Does this agree more with Singer et al. (2019)?\n(5) The authors should also vary the number of planetesimals in the initial two populations. Do the models provide any significant depletion of the outer (cold) population?", "reviewer_2": "Review of NCOMMS-23-32884\nThe authors present a concise description of some collisional evolution simulations that they have done which simultaneously accounts for the varied dynamical history of the proto Kuiper Belt. The outputs of their simulations produce size distributions which broadly match observations presented by Petit et al.\nThe authors make a bold claim that they simultaneously reproduce the size distributions of both the cold classical and dynamically excited populations. This is a disingenuous claim that oversells what they have actually found. Such a result is exactly what is to be expected given the input dynamics and the shattering model they have chosen, so as to be broadly compatible with other observations such as the pluto cratering record. That is to say, little evolution is expected for the cold population, and by nature of their adopted strength model, similar results are expected for the hot population. They haven't started from a vastly different size distribution that converges on the current. Their result is better stated that the features in the current size distributions of both populations that match expectations from the streaming instability model are preserved given expectations from the dynamics we know to have occurred, and with an appropriately chosen strength model. This is a much more truthful - and much less exciting - result than what they claim. At best, they have put together the pieces, and got out exactly what one would expect.\nThe authors also make the claim that they start with the same size distribution for both populations, and end up with matches to the current observations. This isn't actually true, as they require that the hot distribution start with a power law feature for large bodies that doesn't exist in the cold population. If they had actually created this feature from a common starting point, that would be a very exciting result, but instead they arbitrarily adopt a curve that matches current observations without any physical justification, but is required as their model doesn't naturally reproduce it. Nor should it as this feature looks accretional in nature rather than destructive. They simply do not start from the same place as the authors state.\nThe authors make no effort to admit the uncertainties in the observed size distribution, but rather only present the nominal measurement. This hides a property of their initial size distribution - it is already within the uncertainty curve quoted by Petit et al for the cold population, and is no further away from the observation than the final outputs of their simulation. Their initial curve would match the observations just as well as the final result, but with no evolution what so ever.\nThere is a potentially interesting result that the authors do not discuss, which also forms the basis of their calculations. Namely the strength scaling laws they require that differ significantly from those of a decade ago. This is a critical piece of their calculations, and cannot be presented in an unreadable future publication, but rather absolutely must be presented and justified here.\nFinally, they seem to avoid one of the interesting results they have come up with, where the simulated distribution of the hot population is incompatible with previous size distribution efforts derived from the pluto cratering record. Rather than address this interesting result, they decide to punt to a presumably future effort.", "reviewer_3": "The paper by Benavidez et al. presents an interpretation of the of the observed size-frequency distribution (SFD) of the Kuiper belt -- by means of a collisional model (they call it \"Alicandep-22\").\nThey do not model the accretion phase, only the fragmentation phase.\nIn particular, they focus on the dynamical instability as predicted by the Nice model (e.g., Morbidelli & Nesvorny 2020).\nSimply put, the authors show that if the proposed instability is early (10-30 My), there is not enough time -- the collisional evolution of the Kuiper belt is almost negligible, especially for D > 100 km; the initial and final SFDs are almost the same (see their Fig. 1). Only for D < 1 km, there is some evolution (Morbidelli et al. 2021). The only substantial difference between the two SFDs is a dynamical depletion of the massive inner disk, by two orders of magnitude (Fig. 1).\nConsequently, this is a more-or-less negative result.\nA question is, whether this result is suitable for the Nature Comm. journal? I'd say that even if the result is negative, it should have important implications, it should address some problem/contradiction. However, even authors themselves conclude with statements like \"as described in recent literature\" and \"recalled in this research\".\nConsequently, I also do not think it is a 'break through'.\nAnother question is, whether this result is entirely new? The paper by Benavidez et al. (2022; MNRAS, 514, 4876) has a very similar contents. It is not explained, what is the difference between the two papers by the same author(s). Moreover, it is also conceptually similar to Morbidelli et al. (2021). It is not explained, why a new model is necessary and what are the differences between the two models.\n--\nMajor points can be summarized as follows.\nAbstract:\nUnfortunately, the abstract does not have a suitable structure for a Nature-family journal. In particular, at the beginning, it is not clear, what is the problem to solve. The authors mention:\ni) hot vs. cold KB, ii) pebble accretion + streaming instability, iii) 30 vs. 40 au, iv) dynamical instability; but not the problem itself.\nAbstract should also report the results, but the authors only mention one assumption. May be, the results is \"to reconcile\"?\nIt is simply not clear...\nline no. 52\nI believe that the description of the streaming instability (SI) is incorrect. Pebbles are accumulated by aerodynamic drag into planetesimals.\n125\n\"We prefer this SFD to shallower ones...\"\n<- I strongly believe the small-size end is more related to the Q^*, not initial conditions, if there is some collisional evolution...\n133\n\"However, slight discrepancies appear at sizes larger than ~ 500 km, mainly due to the Poisson random statistics...\"\n<- Consequently, it is not possible to run just one model, but -let's say- 10 (or even better 100) to determine the mean value and its uncertainty.\nCould the authors compute also the statistics, please?\n147-158\nD < 2 km slope change is not explained at all. What is the reason for it in your model, please? One also cannot say that Morbidelli et al. (2021) is better or worse, instead, one should recall what are the differences (in physics, not only outcomes) between the two models.\n197-212\nThe discussion (Sec. 3) is focussed on the streaming instability (SI) or pebble accretion (PA), but neither of these processes is included in the Alicandep, which only models late phases of evolution. Consequently, it seems there are no (direct) implications of this work for SI or PA (cf. the statement \"Our findings perfectly match this theory, ...\").\n--\nMinor points can be summarized as follows.\nline no. 30\n\"the cold component ... a steeper slope than the hot component.\"\n<- Please, quantify the slopes and their uncertainties.\n33\n\"Explaining the distinct initial SFDs...\" <- It seems to be a premature statement. Here, one should still talk about the observed SFDs?\n36\n\"original disc surrounding the Sun\" 22-30 au\n<- it's not the \"original\" disk\n<- it's the planetesimal disk\n<- it's surrounding the giant planets, remaining after their formation\nMention high mass (30 M_E) of the planetesimal disk somewhere here?\n40-42\nPlease, do not use abbreviations so often:\nHKB = hot KB\nCKB = cold KB\n47\n\"within the same size range\"\n<- Which s. r. do you refer to, please?\n49\nplanetary -> planetesimal?\n55\nexplain \"polydisperse\" somewhere?\n76\n\"Here we contribute to constraining such a population by comparing our outcome to estimations extrapolated from observational restrictions.\"\n<- I'm sorry, but this sentence seems to me too complicated...\n80\n\"comprehensively reassesses\" <- better wording?\n85\ntrans-planetary == inner == primordial == original?\n<- Please, use consistent wording throughout the paper!\n94\n\"able to replicate\" <- better wording?\n95\n\"achievement is accomplished\" <- better wording?\n99\n\"Notably, such maximum size is determined by the locally available initial mass.\"\n<- Please, add a ref. (this is not from your model).\n110\n\"focus our comparison on\" <- better wording?\n123\n\"across all comparable size ranges\"\n<- What s. r. do you refer to?\n174\nuncovered -> realized?\n175\nreplication -> match?\n252\n\"Petit & Farinella and later updates.\" <- Which updates, please?\nThe word \"corresponding\" is misleading, because in Fig. 4, it is not corresponding to your D_min = 200 m.\n303\nunder ... to -> upon ... from\nFig. 1\nIt is useful, when captions contain not only a description, but a substantial statement, like \"Why this Fig. is used here?\".\nWhat is demonstrated here -- possibly, a negligible collisional evolution for D > 10 km bodies?\nFig. 2\nditto\nFig. 3\nditto\nWhich Q^* is preferred? Is Marschall et al. used? Etc.\nFig. 4\nThis figure is not very informative. Will be the evolved SFDs different for different initial conditions, or not?\nBtw., is the Yarkovsky effect accounted for in your Kuiper belt model?\n--\nNormally (A&A, MNRAS, Icarus, ...), I'd suggest a substantial revision of the paper. However, in the case of Nature Comm., referees should be even more 'strict', right? Here, I leave the decision to the editor...\nWith kind regards, Miroslav Broz"}, "author_response": "Response to Reviewer of the paper: “Deciphering the Size Distribution of the Primordial\ntrans-Planetary Disc”\nWe sincerely thank the reviewers for their comments and feedback, which helped us\nsubstantially improve the quality of the revised manuscript.\nWe have considered all reviewers' comments, answered their questions, and made the\nrequired changes. Additionally, we conducted several new simulations to meet specific\nrequirements and address queries raised by the reviewers.\nWe provide below our rebuttal to the reviewers' comments and remarks.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe manuscript by Benavidez et al. concerns the size distribution of planetesimals in the outer\nSolar System. The authors assume a primordial size distribution based on the streaming instability\nfor the cold classicals but include an additional power law for large planetesimals in the inner\nregions. The authors then run a dynamical model for the planets and planetesimals of the Solar\nSystem and show that they match well the current populations in the hot and cold components of\nthe trans-Neptunian populations.\nI find that the result presented here is potentially important, since it addresses the connection\nbetween the primordial planetesimal sizes and the observed sizes in the outer Solar System. But\nthe narrative and the figures are currently presented in a quite confusing way that makes the\nmanuscript very hard to penetrate. I have some requests here for significant rewrites of the paper\nfor more clarity. I encourage the authors to restructure their manuscript substantially in order to\nelucidate the importance of their result.\nWe thank the reviewer for their comments and suggestions throughout the manuscript.\nAction: We have fully addressed this requirement. The results and discussion sections have been\ncompletely reorganized with a new structure.\n(1) The details of the model here are very important for judging the veracity of the results and\nthus can not only be hidden away in the methods section. The main manuscript should contain a\ndescription of the dynamical model, what physics is included, what is the time-scale of the\ndynamical instability, what is the duration of the simulations (4.6 Gyr?) and how the results\ndepend on some of these choices.\nWe appreciate these suggestions. We agree with the reviewer that having a general outline of the\nmodel and some of the main quantitative details helps. In the revised version, we tried to do so,\nkeeping a balance between giving enough insight into the model and avoiding repetition of the\nthorough description and dependence on initial conditions given in Benavidez et al. (2022).\nAction: We included the requested description of the model and postponed the analysis of the\ndependence of some model parameters to the methods section for which we ran some extra\nsimulations and made two new plots.\n(2) The authors quote Morbidelli et al. (2021) as a key comparison paper. However, how their\nresults differ or compare with M21 is in the end stated only in quite terms. I encourage the authors\nto make their differential analysis with M21 much more clear. Do these results here support M21\nor are they more aligned with Singer et al. (2019)? This is an important debate in the community\nthat this paper can contribute to.\nWe agree it may be worth comparing with Morbidelli et al. (2021) and Singer et al. (2019).\nAction: We created a new section, “Considerations in the Kuiper Belt Populations”, comparing\nMorbidelli et al. (2021) and Singer et al. (2019). We included, in Fig 2, the (scaled up) impactor\npopulation estimated by Singer et al., using their eq. S2. Note that we created this population\nusing the data for Vulcan Planitia (VP) in their supplementary material. We used that set because\nthey consider the cratering record on VP to represent the primary impactor population best.\nThis section also discusses the challenges of achieving a good match between impactor estimation\nSFD and uncertainties on crater counting.\n(3) The caption of Figure 1 is way too short. Please expand here on all the simulation parameters\nindicated in the legend (QS*, SL1, alpha). What is the time when this is plotted?\nAction: The caption of the figure was expanded to include all requested information.\nHow does the cumulative number evolve with time? What is the role of collisional grinding\nrelative to dynamical depletion? Also, the authors should demonstrate here or in another plot what\nhappens for different instability time-scales. Can the authors exclude very long instability time-\nscales from this analysis?\nAction: We performed additional simulation for this request and the results are presented in a\nnew subsection in Methods, “Time dependence”, where we show the effect on the SFD of\ndelaying the dynamical instability onset to (t_i=10, 30, 50 and 100 Myr) considering Dt=10 Myr,\nand also tested a longer Dt=30 Myr. For all runs in this section, we considered our weaker scaling\nlaw. As requested, we included a new plot (Fig. 7) for the SFD at different key times: t=0, t_i\nand\nt=4500 Myr for all cases. We also complemented our analysis with a plot of mass evolution\ncorresponding to different instability onsets (Fig. 8).\nIn more detail, we also attach an extra plot (not included in the manuscript) showing the SFD\nevolution at different times for the case that best matches the constraint (weaker scaling law,\nalpha=0.4, t_i=30 Myr). In this extra plot (Fig. A), we included the SFD at t=0, t_i=30, t=40 and\nt=4500 Myr. It demonstrates that little collisional grinding occurs in the post-instability phase.\nWe invite the reviewer to read the specific section that addresses their questions.\nFigure A: Evolution of the SFD for the total population at different times.\n(4) The caption of Figure 2 is also way too short. What do we learn from the discrepancy between\nthe red dotted line and M21? Does this agree more with Singer et al. (2019)?\nAction: The comparison with Singer et al. (2019) is now addressed in the section \"Considerations\nfor the Kuiper Belt Population\". Please see the discussion in the new manuscript.\nAlso, the caption of Fig. 2 was extended.\n(5) The authors should also vary the number of planetesimals in the initial two populations. Do\nthe models provide any significant depletion of the outer (cold) population?\nWe would politely (partially) decline this request. Let's explain why. Regarding the cold\npopulation, as this outer zone is not strongly affected by dynamical effects, it only loses mass by\ncollisional evolution. Then, if we increase the initial mass, the simulation will end with more mass\nthan expected in this zone due to the low level of collisional grinding. This is due to the relaxed\ndynamical conditions that imply low relative impact velocities.\nRegarding the hot population, we ran some extra cases to show the results of assuming a greater\ninitial mass (45 M_Earth) in the primordial inner belt (Fig. B). Our results show that the shape of\nthe SFD does not change. The reason is that most of the extra mass is depleted by dynamical\neffects. In fact, dynamical models obtain the best results when assuming five giant planets\ninitially, including one ice giant, with a mass comparable to that of Uranus and Neptune, that was\nejected into interstellar space by Jupiter. Dynamical models obtained optimal results when the\nejected planet was placed into the external 3:2 or 4:3 resonances with Saturn and the mass of the\ndisc was 20 M_Earth (Nesvorny and Morbidelli, 2012). Therefore, we consider that successful\nruns must agree with such dynamical constraints and keep about 20 M_Earth at the onset of\ndynamical instability. Therefore, mass depletion at early phases must accomplish that condition,\nwhich prevents extra mass to influence the shape of the distribution collisionally.\n1015\n1010\n105\n100\n10-3 10-2 10-1 100 101 102 103 104\nD (km)\nFigure B: SFD evolution of the inner disc at the onset of dynamical instability increasing the\nnominal mass from 30 to 45 M_Earth, black and blue lines, respectively.\nrebmun\nevitalumuC\nInitial Inner pop. M =30 M\ni\nSFD at t=30 Myr\ni\nReviewer #2 (Remarks to the Author):\nReview of NCOMMS-23-32884\nThe authors present a concise description of some collisional evolution simulations that they have\ndone which simultaneously accounts for the varied dynamical history of the proto Kuiper Belt.\nThe outputs of their simulations produce size distributions which broadly match observations\npresented by Petit et al.\n(1) The authors make a bold claim that they simultaneously reproduce the size distributions of\nboth the cold classical and dynamically excited populations. This is a disingenuous claim that\noversells what they have actually found. Such a result is exactly what is to be expected given the\ninput dynamics and the shattering model they have chosen, so as to be broadly compatible with\nother observations such as the Pluto cratering record. That is to say, little evolution is expected\nfor the cold population, and by nature of their adopted strength model, similar results are expected\nfor the hot population. They haven't started from a vastly different size distribution that converges\non the current. Their result is better stated that the features in the current size distributions of both\npopulations that match expectations from the streaming instability model are preserved given\nexpectations from the dynamics we know to have occurred, and with an appropriately chosen\nstrength model. This is a much more truthful - and much less exciting - result than what they\nclaim. At best, they have put together the pieces, and got out exactly what one would expect.\nThe authors also make the claim that they start with the same size distribution for both\npopulations, and end up with matches to the current observations. This isn't actually true, as they\nrequire that the hot distribution start with a power law feature for large bodies that doesn't exist\nin the cold population. If they had actually created this feature from a common starting point,\nthat would be a very exciting result, but instead they arbitrarily adopt a curve that matches current\nobservations without any physical justification, but is required as their model doesn't naturally\nreproduce it. Nor should it as this feature looks accretional in nature rather than destructive. They\nsimply do not start from the same place as the authors state.\nWe are sorry that we were not able to introduce our main finding suitably. The current model is\nwidely based on Benavidez et al. (2022, B22), where we fully describe our model and explore a\nwide range of initial conditions. In that paper, as in precursor ones (Campo Bagatin and\nBenavidez, 2012, CB12), we show how only given combinations of size distributions and initial\nmass can reproduce observables. We obtained such conclusions by starting from many different\ninitial populations (please, see both papers for details), but only a limited set of parameters can\nmatch current populations. What we do now is to explore the remaining parameter space, starting\nfrom the limited set of initial conditions that allowed CB12, and especially B22, to reproduce\navailable observables. In this case, we look for initial conditions compatible with the Lyra et al.\n(2023) results on pebble accretion at different heliocentric distances. We believe that it is a worthy\noutcome that the same accretion mechanisms -therefore the same shape for the initial size\ndistribution lead -through moderate collisional evolution and significant dynamical depletion- to\ncurrent observables. This match also needs fine-tuning of collisional response. As predicted by\nLyra et al. (2023), accretion of large bodies (>350 km) is possible in the inner region by ice\naccretion but not in the cold region due to low mass density. Therefore, a suitable size distribution\nhas to be assumed there. As little collisional grinding occurred at that end, the current SFD must\nbe a depleted projection of the initial one.\nWe explain how we constructed the initial SFD in the main text and the Methods section. In\nparticular, we try to explain that the size distribution is the same in both the inner and cold\npopulations up to the planetesimals size expected to form at the corresponding heliocentric\ndistance, according to what stated by Lyra et al. (2023). And this is, in fact, one of the main\nmessages of the paper itself. The same formation process can explain the size distribution of\ndifferent populations beyond Neptune, and we show that limited collisional evolution and a\nsuitable start of the dynamical instability led those populations to the current situation. In fact,\nour SFDs are the same for the inner and cold populations up to that size (~100-400 km,\nrespectively). Larger bodies can only accrete in the region where density is high enough.\nAction: We made improvements to the abstract to convey our results accurately and we made\nseveral changes in the structure and in the main text of the manuscript to better explain the\nimplications of our results.\n(2)The authors make no effort to admit the uncertainties in the observed size distribution,\nbut rather only present the nominal measurement. This hides a property of their initial size\ndistribution - it is already within the uncertainty curve quoted by Petit et al for the cold\npopulation, and is no further away from the observation than the final outputs of their\nsimulation. Their initial curve would match the observations just as well as the final result,\nbut with no evolution what so ever.\nWe must acknowledge the reviewer's valid point. Regardless, we can only work with the available\ndata we have access to. However, an extra constraint came from the expected elbow at 1-2 km\nand the SFD slope between 5<D<100 km, which should agree with the Jupiter trojans slope as\nwe discussed in section: “Considerations for the Kuiper Belt Population ''.\n(3)There is a potentially interesting result that the authors do not discuss, which also forms the\nbasis of their calculations. Namely the strength scaling laws they require that differ significantly\nfrom those of a decade ago. This is a critical piece of their calculations, and cannot be presented\nin an unreadable future publication, but rather absolutely must be presented and justified here.\nWe thank the reviewer for this suggestion. This \"unconventional\" scaling law was first pointed\nout by Morbidelli et al. (2021) and then applied to model the Jupiter Trojans by Marschall et al.\n(2022). Therefore, we started experimenting with this scaling law to reproduce Marschall's results\nwith a simplified version of our code adapted to one single population. In this process, we decided\nto adapt our fragmentation code by having Q*_D as a reference to constrain our Q_S in order to\nreplicate Marschall et al. (2022) results. In this way, we suitably adjusted parameters C'_2 and \\\nlambda'_2 in eq. 2 (Methods section). Thus, we achieve results similar to the Boulder code, even\nif we start with a somewhat different approach. Then, we scaled up and down this scaling law\n(dot-dashed line in Figure 3) to perform the simulations in current work. Such previous\n\"calibration\" of the scaling law is part of a Master thesis dedicated to studying the collisional\nevolution of Jupiter Trojans. We evaluated the possibility of including that result in the present\npaper, but we concluded it would be better to publish it in a separate article. The main reason for\nthat decision is that we wanted to focus on the main conclusions of this first paper. We did much\nwork comparing different collisional code outcomes and carried out other additional applications\nof the model applied to Jupiter trojans population, what is a matter of a subsequent paper (in\npreparation). We consider that the relevant point regarding to the scaling law is it weakness and\nthe location of the minimum and this are now fully justified and explained in the new manuscript.\nIn addition, a scaling law similar to this unconventional one was also found in a recent paper on\ncollisional evolution by Bottke et al. (2023). Even if this paper does not use as a constraint the\nKuiper Belt population and rather focuses on reproducing Jupiter Trojans SFD and the cratering\nrecord on giant planet satellites, they tested a large number of scaling laws for Q*_D, finding that\nthey obtained the best fitting results with a similar scaling law to that used in Marschall et al. that\nis the one that we used to calibrate our Q*_S in our whole work. This is extra support for letting\ngo the traditional scaling laws for trans-Neptunian bodies, with a minimum at D~200 m and\nmoving to a minimum at D~20 m.\nAction: We have strengthened the argument about using such unconventional scaling law to\naccomplish this request. This matter is also justified in the new manuscript, including a recent\nreference to Bottke et al. (2023) and provides more details of the used scaling law (Q*_S) to make\nit easier to replicate our results.\nFinally, they seem to avoid one of the interesting results they have come up with, where the\nsimulated distribution of the hot population is incompatible with previous size distribution efforts\nderived from the pluto cratering record. Rather than address this interesting result, they decide to\npunt to a presumably future effort.\nThe referee is likely making reference to Singer et al. (2019), as no specific paper is mentioned.\nThis comment triggered a new analysis of both M21 and S19 SFDs in order to compare them with\nour results.\nAction: This is now addressed in section \"Considerations on the Kuiper Belt Population\".\nReviewer #3 (Remarks to the Author):\nThe paper by Benavidez et al. presents an interpretation of the of the observed size-frequency\ndistribution (SFD) of the Kuiper belt by means of a collisional model (they call it \"Alicandep-\n22\").\nThey do not model the accretion phase, only the fragmentation phase.\nIn particular, they focus on the dynamical instability as predicted by the Nice model (e.g.,\nMorbidelli & Nesvorny 2020).\nSimply put, the authors show that if the proposed instability is early (10-30 My), there is not\nenough time -- the collisional evolution of the Kuiper belt is almost negligible, especially for D\n> 100 km; the initial and final SFDs are almost the same (see their Fig. 1). Only for D < 1 km,\nthere is some evolution (Morbidelli et al. 2021). The only substantial difference between the two\nSFDs is a dynamical depletion of the massive inner disk, by two orders of magnitude (Fig. 1).\nConsequently, this is a more-or-less negative result.\nA question is, whether this result is suitable for the Nature Comm. journal? I'd say that even if the\nresult is negative, it should have important implications, it should address some\nproblem/contradiction. However, even authors themselves conclude with statements like \"as\ndescribed in recent literature\" and \"recalled in this research\".\nConsequently, I also do not think it is a 'break through'.\nAnother question is, whether this result is entirely new? The paper by Benavidez et al. (2022;\nMNRAS, 514, 4876) has a very similar contents. It is not explained, what is the difference\nbetween the two papers by the same author(s). Moreover, it is also conceptually similar to\nMorbidelli et al. (2021). It is not explained, why a new model is necessary and what are the\ndifferences between the two models.\nWe are afraid that we were not able to state the main point of our paper clearly, and we apologize\nfor that. We have deeply modified the argument layout in the paper now but let us clarify that\nhere as well. Using the last updated version of ALICANDEP-22 (see details below and in the\nmanuscript), we are able to model the collisional and dynamical evolution of both the cold and\ninner populations with the same functional form, up to an intermediate size of around 100-400\nkm, accordingly to what Lyra et al. (2023) state for planetesimal formation. In addition, we\ninclude a population of larger bodies in the inner population, following a power law distribution,\nas some simulations of the accretion process glimpse (Johansen et al., 2015). In this way, even if\nthe early onset of dynamical instability limits collisional evolution, we match current observables\nand make some predictions on the populations (SFD) of small objects, a matter of debate in the\ncommunity.\nIn our previous paper (Benavidez et al. 2022), we updated our collisional model to include the\nfeatures of the new Nice model, mainly including a new (inner) zone (22-30 au) to simulate the\nprimordial disc, which was depleted during Neptune outwards migration. In such a paper, we\nstudied four different initial SFDs, two adopted from Nesvorny ́ & Vokrouhlicky ́ (2016) and\nNesvorny ́ & Vokrouhlicky ́ (2019), and two power laws. Also, we experimented with different\nscaling laws. The outcomes of this study were compared to the most updated observables\nestimation at that moment, which was only the number of bodies with D>100 km in each\ndynamical population of the KB.\nHere we go a step further, triggered by the results on pebble accretion by Lyra et al. (2023). Based\non what we found in our previous paper and considering novel observational evidence (Petit et\nal. 2023) -suggesting a similar SFD for the hot and cold population in the range of 100 to 400\nkm- we used our model to test the following hypothesis. Is it possible that the primordial disc\nfrom which the hot population was implanted in the KB was characterized by a functional form\nsimilar to the cold population, consistent with the planetesimal formation process suggested by\nLyra et al. (2023) for bodies up to ~400 km?\nThus, we propose a new initial SFD informed by the functional form of the absolute magnitude\ndistribution of the cold KB population, which is thought to be primordial (Kavelaars et al, 2021).\nThis functional form must be scaled up in order to represent the primordial population between\n22-30 au. For large (>350 km) bodies, a mere scaling cannot ensure attaining the targeted mass\nof 30 M_E and the number of dwarf-planet-like bodies in the dynamically excited population as\nPlutinos and HKB. For this reason, it is necessary to combine the functional form common to the\ninner and cold populations, with a power law at the large-size end (D>100-400 km). Therefore, a\nsignificant part of the work involved finding the optimal combination of functions to match\ncurrent SFDs, assuming initial and evolutive conditions consistent with the Nice model.\nThe collisional model used in this work is not new, instead it has been widely tested in previous\nstudies.\nThe new finding is that a unique functional form of the initial SFD of the primordial disc (up to\nintermediate size, D<100-400 km), can accurately match the observed SFD of both the hot and\ncold populations under suitable conditions. This is substantially different from what was tried in\nthe past, when simply matching the number of objects with a diameter greater than 100 km was\nlooked for, starting from double-slope power law distributions adapted separately to each region.\nThis apparently ‘simple’ solution has two important implications. First, it contrasts with the power\nlaw approach used until now or other arbitrary distributions. Second (and more exciting), the\nproposed functional form links what is suggested by the streaming instability (SI) models and\npebble accretion (Lyra et al., 2023) for sizes up to ~ 400 km, with further accretion in the inner\nbelt, where mass density was high enough to start accretion of ices. Combining this function with\na power law for larger bodies indicates that the formation of dwarf planet-type bodies follows a\ndifferent process than the formation of bodies smaller than ~400 km.\nAction: We now explicitly comment in the paper on the difference between the current paper and\nBenavidez et al (2022). Also, we made several revisions in the manuscript, enhancing the\nrelevance of our findings.\n--\nMajor points can be summarized as follows.\nAbstract:\nUnfortunately, the abstract does not have a suitable structure for a Nature-family journal. In\nparticular, at the beginning, it is not clear, what is the problem to solve. The authors mention:\ni) hot vs. cold KB,\nii) pebble accretion + streaming instability,\niii) 30 vs. 40 au,\niv) dynamical instability; but not the problem itself.\nAbstract should also report the results, but the authors only mention one assumption. May be, the\nresults is \"to reconcile\"?\nIt is simply not clear...\nWe thank the reviewer for pointing this out.\nAction: we reworded the abstract following the Nature structure and clearly stated the problem\nand the main result.\nline no. 52\nI believe that the description of the streaming instability (SI) is incorrect. Pebbles are accumulated\nby aerodynamic drag into planetesimals.\nThanks for pointing this out.\nAction. We reviewed the definition.\n125\n\"We prefer this SFD to shallower ones...\"\n<- I strongly believe the small-size end is more related to the Q^*, not initial conditions, if there is some collisional evolution...\nWe mean that if the SFD starts with a shallower slope for the small-size end, it leads to a shallower\nfinal SFD than the steeper one. That is because the collisional evolution is inefficient, steeping\nthe 1 to 100 km slope to the expected slope inferred for the Jupiter trojans. We agree that the\nminimum of the scaling law induces a dip in the SFD around the same size range, and the slope\nat the small bodies-end is modified as a by-product of collisional evolution.\nAction: We have reworded the text about the effect of a shallower slope combined with modest\ncollisional evolution.\n133\n\"However, slight discrepancies appear at sizes larger than ~ 500 km,mainly due to the Poisson\nrandom statistics...\"\n<- Consequently, it is not possible to run just one model, but -let's say-10 (or even better 100) to\ndetermine the mean value and its uncertainty.\nCould the authors compute also the statistics, please?\nAction: We performed 10 twin runs and estimated the mean value to make the plots. We did not\ninclude error bars in the plot to avoid confusion since they represent four cases.\n147-158\nD < 2 km slope change is not explained at all. What is the reason for it in your model, please?\nOne also cannot say that Morbidelli et al. (2021) is better or worse, instead, one should recall\nwhat are the differences (in physics, not only outcomes) between the two models.\nWe thank the reviewer for their comment.\nAction: we now explain the comparison with Morbidelli et al. (2021) in the line suggested by\nthe referee. In addition, we also include a comparison with Singer et al. (2019). As we also were\nasked to modify the structure of the manuscript this text is now in Sec. “Considerations for the\nKuiper Belt Population”\nyour Kuiper belt model, please? In Tab. 1, you only summarize the values of a, e, i.\nAction: We included the requested data in Tables 1 and 2. We also included in Methods a new\nsection, Impact velocity and intrinsic probability, to describe in more detail the reviewer's\nconcerns on this matter. We also plotted the evolution of the relative impact velocity and intrinsic\nprobability for the inner belt, where most collisional evolution occurred.\nWe did not show the collisional rates because we had to fix a specific size to estimate it, and the\ncode does not save this data automatically. However, we provided all other requested information\nand included a plot of the mass evolution that helps understand how much mass is lost by\ncollisions and dynamical effects.\n197-212\nThe discussion (Sec. 3) is focussed on the streaming instability (SI) or pebble accretion (PA), but\nneither of these processes is included in the Alicandep, which only models late phases of\nevolution. Consequently, it seems there are no (direct) implications of this work for SI or PA\n(cf. the statement \"Our findings perfectly match this theory, ...\").\nThis is correct, our model doesn’t try to reproduce planetesimal formation. What we mean is that,\naccording to planetesimal formation theories, in particular by Lyra et al. (2023), ALICANDEP-\n22 is able to match current observables for both the inner and outer populations starting from\nsuitable initial SFDs, which are common to both up below 400 km. This is possible only if the\nprocess that created such a distribution is also common to both regions, as stated by Lyra et al.\n(2023).\nAction: we addressed this point in a new Sec. “Conclusions”.\n--\nMinor points can be summarized as follows.\nline no. 30\n\"the cold component ... a steeper slope than the hot component.\"\n<- Please, quantify the slopes and their uncertainties.\nAction: we have included the slope values for the hot and cold populations.\n33\n\"Explaining the distinct initial SFDs...\" <- It seems to be a premature statement. Here, one should\nstill talk about the observed SFDs?\nAction: We deleted this sentence in the new manuscript.\n36\n\"original disc surrounding the Sun\" 22-30 au\n<- it's not the \"original\" disk\n<- it's the planetesimal disk\n<- it's surrounding the giant planets, remaining after their formation\nAction: we reworded the sentence as the referee suggested.\nMention high mass (30 M_E) of the planetesimal disk somewhere here?\nAction: we added this data.\n40-42\nPlease, do not use abbreviations so often:\nHKB = hot KB\nCKB = cold KB\nAction: we modified following the suggestion.\nAction: This data is mentioned now.\n47\n\"within the same size range\"\n<- Which s. r. do you refer to, please?\nWe refer here to the size range with common Hr function corresponding to the size range of 100-\n400 km .\nAction: We included that info in the manuscript.\n49\nplanetary -> planetesimal?\nAction: thanks for pointing this out. We change to “planetesimal”.\n55\nexplain \"polydisperse\" somewhere?\nAction: the text was reworded to include the definition.\n76\n\"Here we contribute to constraining such a population by comparing our outcome to estimations\nextrapolated from observational restrictions.\"\n<- I'm sorry, but this sentence seems to me too complicated...\nAction: we reworded this paragraph.\n80\n\"comprehensively reassesses\" <- better wording?\nAction: Reworded.\n85\ntrans-planetary == inner == primordial == original?\n<- Please, use consistent wording throughout the paper!\nAction: We went through the text using consistent wording.\n94\n\"able to replicate\" <- better wording?\nAction: we replaced for “match”.\n95\n\"achievement is accomplished\" <- better wording?\nAction: Reworded.\n99\n\"Notably, such maximum size is determined by the locally available initial mass.\"\n<- Please, add a ref. (this is not from your model).\nThis phrase is not in the new manuscript.\n110\n\"focus our comparison on\" <- better wording?\nThis phrase is not in the new manuscript.\n123\n\"across all comparable size ranges\"\n<- What s. r. do you refer to?\nWe mean here the size range where there is observational data.\nAction: we reworded as: across the comparable size range (\\sim 70 <D< 400 km).\n174\nuncovered -> realized?\nAction: we changed to realized.\n175\nreplication -> match?\nAction: we changed to match.\n252\n\"Petit & Farinella and later updates.\" <- Which updates, please?\nWe included the correction detected by O'Brien and Greenberg, Icarus 178 (2005) 179-212, in\nthe estimation of the potential energy (W_h) of the fragments separated by a distance on the order\nof the Hill radius. Eq 13 in Petit and Farinella (1993) is only valid for the case when both bodies\nare shattered, while O'Brien and Greenberg (2005) developed specific equations for the case in\nwhich one body is shattered and the other one is cratered, and for the case in which both bodies\nare cratered.\nAction: We included a reference to O'Brien and Greenberg (2005).\nwith D_min = 20 m or not? The word \"corresponding\" is misleading, because in Fig. 4, it is not\ncorresponding to your D_min = 200 m.\nWe understand the confusion at this point. Marschall et al. use Q*_D scaling law, while our\nfragmentation code uses Q*_S.\nIt is important to note that while Q*_D and Q*_S are equal in the strength regime (small bodies),\nthere is a difference in the gravity regime (large bodies). This difference occurs because, for large\nbodies, it requires more energy to shatter a target and disperse fragments than simply shatter it.\nThus, there is a moderate gap in the gravity regime, moving the minimum of Q*_S to larger\ndiameters.\nAction: We included a more detailed explanation of the scaling law used in Methods.\n303\nunder ... to -> upon ... from\nAction: fixed\nFig. 1\nIt is useful, when captions contain not only a description, but a substantial statement, like \"Why\nthis Fig. is used here?\". What is demonstrated here -- possibly, a negligible collisional evolution\nfor D > 10 km bodies?\nAction: The figure caption was rewritten to provide a more detailed explanation of the plot's\nmeaning.\nFig. 2\nditto\nAction: The figure caption was rewritten to provide a more detailed explanation of the plot's\nmeaning.\nFig. 3\nditto\nWhich Q^* is preferred? Is Marschall et al. used? Etc.\nAction: The figure caption was rewritten to provide a more detailed explanation of the plot's\nmeaning.\nFig. 4\nThis figure is not very informative. Will be the evolved SFDs different for different initial\nconditions, or not?\nAction: This figure only shows the initial SFD. To make it more illustrative, the two functions\nthat combine to form the initial distribution of the inner disk have been included.\nBtw., is the Yarkovsky effect accounted for in your Kuiper belt model?\nAnswer: The Yarkovsky effect is not taken into account in our model for a number of reasons. 1)\nSuch an effect is interesting in the main asteroid belt (MAB) as it drives bodies into many mean\nmotion and/or secular resonances with Jupiter and Mars present in that region (between 2.1 and\n3.6 au). Instead, in the KB, only Neptune resonances are of interest, but they are separated by 5\nto 7.5 au each other in the current belt. 2) In the MAB, the maximum semi-major axis drift is\nbetween 0.01 and 0.1 au, depending on surface conductivity values. In the current KB, the solar\nradiation flux, F, is ~400 times smaller than in the MAB. The intrinsic impact probability is about\n4 orders of magnitude smaller in the current KB than in the MAB, allowing for Yarkovsky drift\nundisturbed action, so one may argue that it may act over Gyrs instead of Myrs, allowing\naltogether for semiaxes drift on the order of a few au, so only a small fraction of objects already\nclose to resonances may be disturbed by the affect. Moreover, at the time when most of the\ncollisional evolution took place, the drift over collisional time would be much smaller because of\nlarger collisional probability, and the only resonance in the primordial belt was a 4:3 mmr with\nNeptune. Again, only a small fraction of bodies might reach that resonance at the corresponding\npossible drift rates (<<1 au, compared to 8-10 au width of the belt). Such drifts shouldn’t make\nany significant change on overall populations and may just affect a fraction of objects that were\nalready very close to that resonance, with little or no impact on the number of surviving bodies\nafter the dynamical instability.\n--\nNormally (A&A, MNRAS, Icarus, ...), I'd suggest a substantial revisión of the paper. However,\nin the case of Nature Comm., referees should be even more 'strict', right? Here, I leave the decisión\nto the editor.\nWith kind regards, Miroslav Broz"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have done a very comprehensive job at revising the paper. I am very satisfied that all my comments were taken into account. I am therefore happy to recommend the paper for publication in Nature Communications once these three minor points are addressed:\n(1) It is confusing that the authors write the inventory of the Kuiper Belt is complete up to Hr=5.3 in the first sentence of the paper, and then later write that it is relatively complete up to Hr=8.3. The authors should make sure that these two sentences are consistent with each other.\n(2) The legend in Figure 1 overlaps with the lines. This should be avoided for clarity in reading the figure.\n(3) The error interval on the S19 data in Figure 3 is very useful, but could be marked better with a grey background or similar.\nAlso, why does the data correspond to the upper limit around D ~ 1 km? And what is the meaning of the inverted error limits beyond D ~ 2 km?", "reviewer_2": "The authors have done a decent job in altering the claims in the manuscript to reflect what they have found from their simulations. They have also improved the presentation of the modelling, and the assumptions made therein. In that sense, the paper presents a lot better and lot more complete than the original version.\nI still have numerous significant concerns however that have not been addressed:\n1. The authors have failed to include a robust discussion of uncertainties. This goes two ways. Observationally, the data they are matching their models against (the OSSOS results) have significant uncertainties regarding detection rates that the authors do not make any use of. I have checked with Petit and Kavelaars and their publications and the uncertainties are available their, or through an email (which it seems the authors must have done at some point already).\nSimilarly there is no robust analysis of numerical uncertainties, most importantly the effects of model assumptions on the outputs. They have done a better job of discussing this, but only present a few models, which is hardly a complete discussion.\nThe major outstanding question is how far can the nominal model be bent and still be consistent with the range of possible SFDs allowable by the OSSOS data?\nThe above question immediately reflects upon whether or not the author's simulation outputs are common to a large range of possible models, or require very particular assumptions to get it \"just right\". Either scenario would be interesting, but as it stands, the reader cannot currently know this.\n2. Similarly to point 1 above, FIgure 3 has data points with no uncertainties, and curves with no uncertainties, which makes the presentation disingenuous. Given the uncertainties on the crater counts, it is reasonable to expect that a number of different models could be made to fit the possible range of crater count values, not just the singular model they argue for.\n3. I didn't raise this issue with the first review because it was already too long, so I raise it here. There are numerous other sources of measurement the authors have not considered. In particular the Bernstein 2004 result which presents the smallest size objects for which a direct SFD can be measured. Similarly, the bright end can be further constrained with the MPC dataset which is mostly complete at sizes larger than the OSSOS upper size range. Both of these results must be matched for the model to stand. Even an approximate match in the way done in Kavelaars 2021 would be sufficient, and so this shouldn't take too much effort.\nOther more minor concerns:\n4. The authors should state how many 1000+ km objects they form in their model outputs. This is an interesting number.\n5. The first sentence of the second paragraph on page 3 reads \" The inventory of large bodies in the Kuiper belt seems nowadays relatively complete for (cid:1)(cid:2) < 8.3, \"\nReasonable estimates of the number of H<8.3 objects predict an oder of magnitude more objects in that size range than the entire known sample. We provide a few here:\nFrom the early works of Gladman, Jewitt, Fraser, Fuentes, and Jones, the apparent magnitude distribution has slope alpha~0.7. Using the normalizations of Fraser or Fuentes, the cumulative distribution is approximately N(<m_r)~10^(0.7*(m_r-23.4)) objects per square degree. H_r~8.3 corresponds to m_r~25 at classical KBO distances (44 au). That results in 13 objects per square degree. The region of highest density is +-6 degrees of the ecliptic, resulting in 13 x (6x2) x 360 = 56,000 KBOs.\nFrom Fraser et al. 2014, the ecliptic number density of objects with H_r<8.3 ~7 per square degree, or 7 x (6x2) x 360 = 30,000 objects near the ecliptic.\nFrom Kavelaars 2021, directly sampling off their figure 2, they infer 10,000 bodies with H_r<8.3 in the cold classicals alone.\nThe important first sentence clearly needs correction.\n6. In the observational constraints section, they quote assumed albedos with no references. Obvious choices could include the Fraser 2014 result which shows albedos against size and the Lacerda 2014 result which shows albedos against color, or the appropriate TNOs are Cool papers (and the website)", "reviewer_3": "Dear editor, dear authors,\nI am sorry for a delayed response.\nRegarding the revised manuscript, I have to thank the authors\nfor revising their manuscript and explaining of some the points.\nAt the same time, I have to say that some of them still remain\nvalid. Since the collisional evolution of the Kuiper belt is\nalmost negligible (in the early-instability scenario), their\nwork is adjusting the initial conditions (\"functional forms\")\nin a way which matches the observations. The major driver of\nthe evolution is dynamics, which is size-independent (as seen\nin their Fig. 1).\nConsequently, I am sorry for this statement, but I still think\nthat the novelty of this work is not enough to warrant a publication\nin the world-leading journal like Nature Communications.\nThe authors themselves admit that the differences with respect to\nBenavidez et al. (2022; MNRAS, 514, 4876) are not so substantial.\n(See \"Here we go a step further, ...\" or l. 118 from the manuscript\n\"We show that bodies larger than 10 km currently preserve the\nsignature of the formation process, while smaller bodies were\njust slightly affected by collisional evolution in the early phase\nwhen Neptune started its outward migration...\")\nThe fact that the SFDs of hot and cold KB are similar has been\nalready discussed from the observational point of view in\nPetit et al (2023; ApJL 947, 4).\nFinally, let me briefly comment on the fact that the manuscript\nhas been rewritten substantially (most of the text is now \\bf).\nThe discussion has been extended to include also constraints\nat small sizes (D ~ 1 km), where the collisional evolution is\nalready apparent on the SFD. Its shape is controlled by the\nspecific shattering energy Q^* . This is a nice addition, although\nit does not improve the concerns raised above.\nWith kind regards, Miroslav Broz"}, "author_response": "RESPONSE TO REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have done a very comprehensive job at revising the paper. I am very satisfied\nthat all my comments were taken into account. I am therefore happy to recommend the\npaper for publication in Nature Communications once these three minor points are\naddressed:\n(1) It is confusing that the authors write the inventory of the Kuiper Belt is complete up\nto Hr=5.3 in the first sentence of the paper, and then later write that it is relatively\ncomplete up to Hr=8.3. The authors should make sure that these two sentences are\nconsistent with each other.\nWe thank the reviewer for their positive comment and for bringing this point to our\nattention. We have revised the sentences as indicated in the new manuscript.\n(2) The legend in Figure 1 overlaps with the lines. This should be avoided for clarity in\nreading the figure.\nOk, done.\n(3) The error interval on the S19 data in Figure 3 is very useful, but could be marked\nbetter with a grey background or similar.\nWe modified Figure 3 according to the reviewer's suggestion, including a grey\nbackground. We also added error bars to the Singer et al.’s (2019) data indicating the\nrange of uncertainty.\nAlso, why does the data correspond to the upper limit around D ~ 1 km?\nThe upper limit of the shaded zone overlaps values close to 1 km due to the flattening of\nthe distribution at small sizes. For a better interpretation, we included now a dashed line\nin Fig 3 that fits the distribution of projectiles smaller than 1.5 km, corresponding to a\nslope of 0.7. We chose this size because it is the closest bin to the inflexion point (1-2\nkm) suggested by S19 and Robbins & Singer (2021). Then, it is easy to see that the\ndashed line is in the middle of S19+- 2 \\sigma (corresponding to 0.7+2*(0.2)=1.1 and\n0.7-2*(0.2)=0.3).\nAnd what is the meaning of the inverted error limits beyond D ~ 2 km?\nThanks for pointing this out. In the revised manuscript we plotted the error interval on\nthe S19 data up to 1.5 km, as extending the line to larger sizes has no meaning here.\nReviewer #2 (Remarks to the Author):\nThe authors have done a good job in addressing the various concerns I raised in the previous,\nand first reviews, rather than just punting those issues to \"later works\". The manuscript is in\nexcellent shape and I see no further need to hold up its publication. Commendable work.\nI recommend this paper be accepted.\nReviewer #2 (Remarks on code availability):\nThe code is not available for review.\nWe now provide an executable version of our code and the necessary input files, enabling\nthe reproduction of the results.\nReviewer #3 (Remarks to the Author):\nDear editor, dear authors\nI am sorry for a delayed response.\nRegarding the revised manuscript, I have to thank the authors for revising their\nmanuscript and explaining of some the points.\nAt the same time, I have to say that some of them still remain valid. Since the\ncollisional evolution of the Kuiper belt is almost negligible (in the early-instability\nscenario), their work is adjusting the initial conditions (\"functional forms\") in a way\nwhich matches the observations. The major driver of the evolution is dynamics, which\nis size-independent (as seen in their Fig. 1).\nConsequently, I am sorry for this statement, but I still think that the novelty of this work\nis not enough to warrant a publication in the world-leading journal like Nature\nCommunications.\nThe authors themselves admit that the differences with respect to\nBenavidez et al. (2022; MNRAS, 514, 4876) are not so substantial. (See \"Here we go a step further, ...\"\nor l. 118 from the manuscript \"We show that bodies larger than 10 km currently\npreserve the signature of the formation process, while smaller bodies were just slightly\naffected by collisional evolution in the early phase when Neptune started its outward\nmigration...\")\nThe fact that the SFDs of hot and cold KB are similar has been\nalready discussed from the observational point of view in\nPetit et al (2023; ApJL 947, 4).\nWe respect the opinion of the referee but let us present some arguments supporting\nthe novelty of our work. The recent evidence that both the hot and cold populations share\na similar size distribution within a specific size range suggests a common planetesimal\nformation process. This hot topic in planetary sciences contradicts what was previously\nbelieved in the last decade. This hypothesis, proposed by Petit et al. (2023), necessitates\nrigorous testing with a proven and reliable collisional and dynamical evolution model,\nwhich is the primary objective of our study. In addition, compelling new theoretical work\nhas recently become available (Lyra, 2022), suggesting a common formation pattern\nbased on polydisperse pebble accretion for the formation of planetesimal up to a few\nhundred km. As stated in the manuscript, the model presented reconciles, for the first\ntime in the literature, to our knowledge, observational evidence to theoretical prediction\nthrough accurate modelling of the overall evolution of the trans-Neptunian populations.\nWe use a collisional and dynamical evolution model, updated in our latest publications,\nto test Petit's hypothesis with the support of theoretical findings. Our findings demonstrate\nsuch a hypothesis and define constraints on the initial conditions for the evolution of the\nprimordial population, refining—in addition—the scaling law characterizing collisional\noutcome.\nOur submitted article represents the first comprehensive study on this topic, leveraging\nthe most recent and advanced literature. We believe that this characteristic makes our\nwork particularly relevant and significant as to be published in a journal such as Nature\nCommunications."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have addressed all my comments in the this revision. I am happy to recommend the paper for publication.", "reviewer_2": "The authors have done a good job in addressing the various concerns I raised in the previous, and first reviews, rather than just punting those issues to \"later works\". The manuscript is in excellent shape and I see no further need to hold up its publication. Commendable work.\nI recommend this paper be accepted.", "reviewer_3": "Yesterday, I had a fresh look at the manuscript by Benavidez et al.,\nits latest version (42676_2_art_file_9375056_sh52xq_convrt.pdf).\nLet me discuss separately the results for the large (>100 km)\nand small (~<10 km) Kuiper belt objects (KBOs).\nFor >100 km objects, the collisional evolution is negligible\n(see Figs. 1, 2, 7, 9). This is correctly admitted by the authors.\nOk, but then their results are not (cannot be) novel, because\ntheir collisional model does \"nothing\". The only exception is\nthe dynamical instability acting on the inner disk -> hot KBOs,\nwhich is size-independent and is modelled as a depletion factor\n(Fig. 1, top panel). Consequently, they essentially proceed as:\ninitial conditions\n|\nV\nno collisions\n|\nV\n\"good match\"\nof observed SFD\nSetting-up initial conditions in accord with observations is\n-to me- trivial. Even using a different functon(exponential\n+ power-law) to approximate the observations is trivial.\nConsequently, I do not think this result alone warrants\na publication in the Nature Communications journal.\nThe fact that collisional evolution is negligible is known and\npublished previously (Benavidez et al. 2022, Figs. 3, 4, 5, 7).\nCollisional probabilities are low (Tab. 2) and the instability\nis assumed to occur early (e.g., rev. by Nesvorny 2018).\n--\nFor ~<10 km objects, there is a serious problem with observations,\nnamely Singer et al. (2019) and Morbidelli et al. (2021) seem to\nbe contradicting each other. For reference, I put the original\nfigs. here:\nhttps://sirrah.troja.mff.cuni.cz/~mira/tmp/benavidez/Singer_2019_Fig2.png\nvs.\nhttps://sirrah.troja.mff.cuni.cz/~mira/tmp/benavidez/Morbidelli_2021_Fig5.png\nThe Singer's error bars are relatively small, because the smallest\ncraters are the most numerous. Different locations (count areas)\nshow some mutual differences, but they all agree on a systematic\ndecrease of the slope below ~5 km (i.e., crater, not projectile).\nBenavidez et al. claim that these observations are uncertain and\nthat their model is in agreement (\"2-sigma\"). No, this is not\npossible. None of their collisional models shows a break together\nwith a continuation of the very shallow slope of Singer et al.\nInstead, the model continues with a steeper slope (see Fig. 3).\nMorbidelli's SFD is smoother, likely because it fulfills more\nobservational constraints (surveys, Trojans, Pluto craters,\ndust abundance). Ok, but none of the four models by Benavidez\never comes close to it. Even the preferred one (with the weakest\nscaling law, \"SL1x0.1 (alpha=0.3)\") is offset (up) between 1-100 km\nand offset (down) between 0.01-0.1 km.\nGenerally, I agree with the methodology, i.e., selecting/modifying\nthe scaling law to fit the observations. I agree with the statement\nthat objects approx. <0.1 km are probably weak.\nHowever, the fit presented in the paper is by no means \"the best fit\".\nThe statistical analysis is missing, all observations should have\nindividual error bars, the model parameters (C_1', lambda_1', C_2',\nlambda_2', alpha_s, alpha_l, beta, H_0, H_B, N_i1, N_i2) should be\nvaried continuously and the statistical significance of results\nshould be estimated (e.g., using the standard bayesian MCMC analysis)\n-- especially if you aim for Nature Communications.\nThis is in agreement also with Referee #1.\n--\nMinor or language points:\nabstract, l. 7\nframe -> framework\nobservational -> observed\np. 2, l. 31\n\"broken power laws with different slopes\" is ambiguous\n2, 33\nWhy the slopes were so different? -7.5 vs. -4.5\n2, 58\n\"aimed at concentrating solid particles\" is not an explanation\nperhaps, mention a difference in speeds (gas vs. solids)?\n3, 67-71\nThere was nothing between 30-40 au, really?!\n3, 78\nadd \"by Singer et al.\"?\n3, 95\nThe model by Benavidez et al. was practically the same, because\nthere was no collisional evolution.\n5, 177-179\nCould bigger bodies in hot KBOs originate from planet-crossing orbits,\nsimilarly as in Nesvorny et al. (2024)? In other words, could they\nbe former protoplanets (left-overs)?\n5, Tab. 1\nThere was nothing between 30-42 au, really?!\n6, 224\n\"OSSOS uncertainties at larger sizes open a wide range of possible\ncombinations for the SFD\" is in contradiction with using the all-sky\nMPC data, cf. the completeness.\n8, 243\nDitto.\n9, 278\n\"... by suitably scaling it to the Kuiper belt distribution\".\nWell, this is not good, because you are discarding the absolute\nnumber of craters, i.e., an important observational constraint!\nInstead, your model should be used to estimate it and to compare\nto observations...\n10, 323-324\n\"This finding highlights the intricate connection between collisional\nprocesses and early dynamical instability phase.\"\nLet me repeat that collisional evolution at >100 km is negligible.\n11, Fig. 3\nWell, this is a third main-text figure which shows essentially\nthe same thing (SFDs). Figs. 1, 2, 3 could be merged in one fig.\n(albeit with multiple panels).\n14, Fig. 5\nSimilar as Figs. 1, 2, 3?\n18, Tab. 3\nIs it initial? Or final?\nAre these all runs that you have performed? Only 5?\n--\nSince the issues discussed above are a bit too serious, I do not\nsupport the publication in the Nature Communications journal;\nat least in the current form. Nevertheless, the decision is on\nthe side of its editors...\nWith kind regards, Miroslav Broz"}, "author_response": "RESPONSE TO REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe authors have addressed all my comments in the this revision. I am happy to\nrecommend the paper for publication.\nAnswer:\nThank you for your kind and positive feedback.\n–\nReviewer #2 (Remarks to the Author):\nThe authors have done a good job in addressing the various concerns I raised in the previous,\nand first reviews, rather than just punting those issues to \"later works\". The manuscript is in\nexcellent shape and I see no further need to hold up its publication. Commendable work.\nI recommend this paper be accepted.\nReviewer #2 (Remarks on code availability):\nI go to the above link and nothing loads. So it seems like there is no code to review?\nAnswer:\nWe truly appreciate your constructive feedback, which has helped improve the\nmanuscript. We are delighted that you find the revised version suitable for publication.\nWe apologize for the inconvenience regarding the code. We have uploaded the\nnecessary files again.\n--\nReviewer #3 (Remarks to the Author):\nYesterday, I had a fresh look at the manuscript by Benavidez et al., its latest version\n(42676_2_art_file_9375056_sh52xq_convrt.pdf).\nLet me discuss separately the results for the large (>100 km) and small (~<10 km)\nKuiper belt objects (KBOs).\nFor >100 km objects, the collisional evolution is negligible (see Figs. 1, 2, 7, 9). This is\ncorrectly admitted by the authors. Ok, but then their results are not (cannot be) novel,\nbecause their collisional model does \"nothing\". The only exception is the dynamical\ninstability acting on the inner disk -> hot KBOs, which is size-independent and is\nmodelled as a depletion factor (Fig. 1, top panel). Consequently, they essentially\nproceed as:\ninitial conditions\n|\nV\nno collisions\n|\nV\n\"good match\"\nof observed SFD\nSetting-up initial conditions in accord with observations is -to me- trivial. Even using a\ndifferent functon(exponential + power-law) to approximate the observations is trivial.\nConsequently, I do not think this result alone warrants a publication in the Nature\nCommunications journal.\nThe fact that collisional evolution is negligible is known and published previously\n(Benavidez et al. 2022, Figs. 3, 4, 5, 7). Collisional probabilities are low (Tab. 2) and the\ninstability is assumed to occur early (e.g., rev. by Nesvorny 2018).\nAnswer:\nWe appreciate the reviewer's feedback and understand that he holds a different\ninterpretation of the novelty of our study.\nWe have thoroughly explained the novelty of our study in previous responses and within\nthe manuscript itself. Specifically, the approach we present involves\nidentifying a suitable size distribution function and initial conditions that align with the\nrequirements of planetesimal formation models and most observational constraints. This\nrepresents a significant advancement compared to prior studies.\nWhile we recognise that some aspects of collisional evolution are negligible, this does\nnot diminish the originality of our contribution. The careful construction of initial\nconditions, encompassing the primordial populations and integrating them with a robust\nphysical framework, remains an innovative aspect of our work.\n--\nFor ~<10 km objects, there is a serious problem with observations,\nnamely Singer et al. (2019) and Morbidelli et al. (2021) seem to\nbe contradicting each other. For reference, I put the original\nfigs. here:\nhttps://sirrah.troja.mff.cuni.cz/~mira/tmp/benavidez/Singer_2019_Fig2.png\nvs.\nhttps://sirrah.troja.mff.cuni.cz/~mira/tmp/benavidez/Morbidelli_2021_Fig5.png\nThe Singer's error bars are relatively small, because the smallest\ncraters are the most numerous. Different locations (count areas)\nshow some mutual differences, but they all agree on a systematic\ndecrease of the slope below ~5 km (i.e., crater, not projectile).\nBenavidez et al. claim that these observations are uncertain and\nthat their model is in agreement (\"2-sigma\"). No, this is not\npossible. None of their collisional models shows a break together\nwith a continuation of the very shallow slope of Singer et al.\nInstead, the model continues with a steeper slope (see Fig. 3).\nMorbidelli's SFD is smoother, likely because it fulfills more\nobservational constraints (surveys, Trojans, Pluto craters,\ndust abundance). Ok, but none of the four models by Benavidez\never comes close to it. Even the preferred one (with the weakest\nscaling law, \"SL1x0.1 (alpha=0.3)\") is offset (up) between 1-100 km\nand offset (down) between 0.01-0.1 km.\nAnswer:\nSinger 19 clearly states in his abstract that there is a scarcity of small craters (<13 km),\nimplying a deficit of small KBOs <1-2 km (all in diameter). The later study by Robbins et al.\n(2021) places a transition range for the distribution of craters between 10-15 km, emphasizing\nthat the slopes corresponding to small craters are quite ambiguous, but that a change in slope\nis robustly confirmed. As it is stated in the manuscript, our model catches properly the slope\nSFD change around 1-2 km. However, both S19 and M21 models have large uncertainties that\ncannot be fitted at the same time.\nFigure 2 highlights how our initial distribution with alpha=0.4 combined with an increasingly\nweaker scaling law reaching that indicated as SL1x0.1 can reproduce the slope change\nobserved at 1-2 km. However, this slope change does not occur if bodies are strong, nor even\nin the case of the weakest bodies tested and an initial shallow distribution (alpha=0.3).\nIn particular, regarding the discrepancy with M21, we would like to recall a recommendation\nby Petit et al. (2023): “Morbidelli et al. (2021) concluded dominance of the hot population for\nArrokoth crater formation, but their numbers of 100 km and 2 km projectiles are not in line\nwith our current estimates and should be revisited. “ This opens the possibility that some\nupdate of Morbidelli’s work is needed.\nTo clarify these points, we have reworded a statement on page 8.\n--\nGenerally, I agree with the methodology, i.e., selecting/modifying\nthe scaling law to fit the observations. I agree with the statement\nthat objects approx. <0.1 km are probably weak.\nHowever, the fit presented in the paper is by no means \"the best fit\".\nThe statistical analysis is missing, all observations should have\nindividual error bars, the model parameters (C_1', lambda_1', C_2',\nlambda_2', alpha_s, alpha_l, beta, H_0, H_B, N_i1, N_i2) should be\nvaried continuously and the statistical significance of results\nshould be estimated (e.g., using the standard bayesian MCMC analysis)\n-- especially if you aim for Nature Communications.\nThis is in agreement also with Referee #1.\nAnswer:\nOur study includes a comprehensive analysis of the relevant parameters, as detailed in\nthe section Model flexibility to match observational uncertainties. In this section, we\nsystematically explored parameter values through extensive simulations, ensuring they\nremained within the observational error bars. This approach allowed us to determine a\nrange of variability for α and β, along with the other parameters, that preserves\nreasonable initial conditions—such as the number of dwarf planets and total mass—\nwhile still producing final distributions consistent with current observations.\nA more in-depth statistical analysis, such as a Bayesian MCMC approach, could be\nvaluable for refining population matching within specific size ranges, where more\nprecise estimates of body counts could be directly compared with modelling results.\nHowever, such an analysis falls beyond the scope of the present study.\nFinally, we have added a paragraph at the end of Section 4 of the manuscript, discussing\nthis potential future contribution and highlighting how it could be integrated into more\nspecialized studies.\n--\nMinor or language points:\nabstract, l. 7\nframe -> framework done\nobservational -> observed done\n--\np. 2, l. 31\n\"broken power laws with different slopes\" is ambiguous\nAnswer:\nWe have reworded this sentence slightly. However, the next sentence in the manuscript\ncompletes the message, where the slope values (at the ends of large bodies) for each\ncomponent (cold and hot) are given.\n--\n2, 33\nWhy the slopes were so different? -7.5 vs. -4.5\nAnswer:\nThis is a known issue in past literature, and indeed has nothing to do with our work. In\nmost surveys before Kavelaars et al. (2021), a broken power law was used to fit the\nSFD of the cold and hot KBO populations. This led to different slopes being estimated\nfor the bright/large body end (-7.5 and -4.5). Thus, the scientific community\ninterpreted this fact as the by-product of different formation processes acting at\ndifferent heliocentric distances.\nInterestingly, this issue is related to the novelty itself of our work.\nThe article introduction points out the relevance of the new fit, to the SFD of KBO\nproposed by Petit et al. (2023), based on an exponential taper function. This new fit is\nconnected to the formation mechanisms of planetesimals, contrasting with the\npreviously employed arbitrary broken power-law fits.\nUnderstanding the relevance of using a different function -an exponentially tapered\nfunction- to fit the observed SFD of both populations up to a certain size, is the key\nelement which supports the novelty of the initial SFD used in our article, as it links it\nwith a well studied formation process.\nWe added some words at the end of the paragraph linking different slopes at large\nsizes to different formation processes.\n--\n2, 58\n\"aimed at concentrating solid particles\" is not an explanation\nperhaps, mention a difference in speeds (gas vs. solids)?\nAnswer:\nReworded as:\nThe SI represents an aerodynamic process driven by the speed difference between\nsolid particles and gas and the back reaction of solids onto the gas (Morbidelli et al.,\n2020).\n--\n3, 67-71\nThere was nothing between 30-40 au, really?!\nAnswer:\nIn the manuscript, we discuss the results of recent theoretical studies. Thus, we just\nmention that Lyra et al., (2023) show that mass accretion to Pluto-size bodies beyond\n30 au is inefficient. We are not proposing any planetesimal formation scenario\nourselves.\n--\n3, 78\nadd \"by Singer et al.\"?\ndone\n--\n3, 95\nThe model by Benavidez et al. was practically the same, because\nthere was no collisional evolution.\nThe model used in our paper in 2022 was similar, but not the same: the scaling laws\ntested and the initial SFD are different. Moreover, in the submitted paper the\nobservational constraints are updated to the new published results of Petit et al.,\n(2023).\n--\n5, 177-179\nCould bigger bodies in hot KBOs originate from planet-crossing orbits,\nsimilarly as in Nesvorny et al. (2024)? In other words, could they\nbe former protoplanets (left-overs)?\nAnswer:\nOur model is based on the framework of various Nesvorny papers cited in the\nmanuscript from references 11 to 16, plus reference 47. According to such papers, the\nprimordial disc, located between 22 and 30 au (our inner region), may contain about\n2000 to 4000 Pluto-size bodies (or larger) that are scattered inward and outward when\nNeptune starts its migration. In a personal communication with Dr. D. Nesvorny, he\nconfirmed that grain migration still works even with 6000 Pluto-size bodies. So, under\nthis framework, bigger bodies in the hot KB could form in our modelled inner disc and\nbe implanted in both the plutinos and hot populations.\n--\n5, Tab. 1\nThere was nothing between 30-42 au, really?!\nAnswer:\nBeyond the arguments already explained at the previous identical question, in this\ncase, for simplicity, our model does not consider an initial size/mass distribution in that\nregion. Any added mass in that region should necessarily be very low to be consistent\nwith Neptune's migration models, so the effect on collisional evolution would be\nnegligible. However, not considering a distribution of bodies in the 30 to 40 AU region\nhas been our hypothesis from the beginning.\n--\n6, 224\n\"OSSOS uncertainties at larger sizes open a wide range of possible\ncombinations for the SFD\" is in contradiction with using the all-sky\nMPC data, cf. the completeness.\nAnswer:\nThanks for pointing this out. We believe that the methodology is clear in the text: our\nbaseline reference for comparison is the MPC at large sizes, but for completeness, we\nalso explored the uncertainty range in the OSSOS++ survey. This is now smoothed in\nthe text to avoid reader’s misunderstanding.\n--\n8, 243\nDitto.\nPlease, see the answer to the previous point.\n--\n9, 278\n\"... by suitably scaling it to the Kuiper belt distribution\".\nWell, this is not good, because you are discarding the absolute\nnumber of craters, i.e., an important observational constraint!\nInstead, your model should be used to estimate it and to compare\nto observations…\nAnswer:\nWe do not agree on this point. Our model follows the evolution of an entire\npopulation, whereas the crater SFD comes from observations of a single target\n(region). Consequently, the key issue is the ratio of projectiles to targets, represented\nby the slope of the SFD within a specific size range. Then, to address this, we need to\nscale, for instance, the crater SFD to match the modelled SFD and compare the slopes\nnear the elbow.\n--\n10, 323-324\n\"This finding highlights the intricate connection between collisional\nprocesses and early dynamical instability phase.\"\nLet me repeat that collisional evolution at >100 km is negligible.\nAnswer:\nIndeed. But collisional evolution takes place for sure at smaller sizes. It is well known\nthat collisional evolution is negligible at large sizes, we never assessed the opposite.\n--\n11, Fig. 3\nWell, this is a third main-text figure which shows essentially\nthe same thing (SFDs). Figs. 1, 2, 3 could be merged in one fig.\n(albeit with multiple panels).\nAnswer:\nWe consider that it is not convenient to put the 3 figures together in a single one\nfigure, since the details of the SFD would not be appreciated.\n--\n14, Fig. 5\nSimilar as Figs. 1, 2, 3?\nAnswer:\nFig. 5 represents the initial SFD in the inner and cold population, it is not similar to Figs\n1, 2 or 3.\n--\n18, Tab. 3\nIs it initial? Or final?\nAre these all runs that you have performed? Only 5?\nAnswer:\nTable 3 presents the parameters used to construct the initial SFD, as described in the\nfirst paragraph of the section Model flexibility to match observational uncertainties.\nWe conducted a systematic exploration of all relevant parameters, performing numerous\nsimulations with a wide range of values for each. However, only specific parameter\ncombinations successfully reproduced the observed distributions. The runs listed in\nTable 3 represent the best matches that satisfy the imposed threshold constraints.\n--\nSince the issues discussed above are a bit too serious, I do not\nsupport the publication in the Nature Communications journal;\nat least in the current form. Nevertheless, the decision is on\nthe side of its editors...\nWith kind regards, Miroslav Broz"}]

Decision letter