Formation of Comets

Questions regarding how primordial or pristine the comets of the solar system are have been an ongoing controversy. In this review, we describe comets' physical evolution from dust and ice grains in the solar nebula to the contemporary small bodies in the outer solar system. This includes the phases of dust agglomeration, the formation of planetesimals, their thermal evolution and the outcomes of collisional processes. We use empirical evidence about comets, in particular from the Rosetta Mission to comet 67P/Churyumov--Gerasimenko, to draw conclusions about the possible thermal and collisional evolution of comets.

with the Keplerian frequency Ω. Particles with the same Stokes number behave aerodynamically similar. For a given PPD model and a fixed heliocentric distance, the Stokes number corresponds to a specific value of R ρ p or R 2 ρ p , depending on the flow regime (see Equation (1)). Thus, for a fixed value of the mass density of the dust particle, the Stokes number represents the particle size. Hence, particles with different sizes (or Stokes numbers) obtain different speeds at the same heliocentric distance, which leads to relative velocities among the solid particles and, consequently, to collisions.
The frictional interaction of gas and dust and the subsequent sub-Keplerian speed of the dust particles result in an inward drift of the dust, with radial velocities depending on Stokes number (i.e., particle size). Very small and very large objects possess radial velocities that are very small; the maximum radial speed is reached for dust particles with St = 1. In general, dust drifts towards higher gas pressures, which also applies for vertical settling towards the disc midplane or local pressure bumps.
Turbulent gas motion influences the dust motion as well. Quantitatively, the turbulence of a disc can be described by the turbulent strength parameter α. Qualitatively, gas turbulence can lead to collisions among the dust particles as well as to particle transport. Thus, thermally processed particles from the inner disc might be transported in-or outward due to turbulence, which leads to radial mixing of the dust species formed at different temperatures.
For very small, i.e., sub-micrometre to micrometre-sized, grains, Brownian motion may also be responsible for verylow-speed collisions. Due to the Maxwell-Boltzmann distribution of the mean thermal kinetic energy of the grains, equal-sized particles can also collide due to Brownian motion [30].

A Growth Scenario from Grains to Pebbles
When dust or ice particles collide, the outcome depends on material properties, particle size and collision velocity [31,32,33,34,35,36,37,38]. The relative velocities between the grains, the number of particles of a given size and the collisional cross-sections determine the number of collisions per unit time. The collision velocities result from the processes that induce the relative motion between dust and gas (see Section 2.1). For radial, azimuthal and vertical drift, the velocity between dust and gas depends on particle size; hence, only particles of different sizes collide, because particles of equal size move with the same speed. However, equally sized particles can collide in case of stochastic motion, as induced by Brownian motion or gas turbulence.
When particles of comparable size collide, the outcome can be sticking, bouncing, abrasion or fragmentation. Sticking always occurs when the collision energy is small compared to the van der Waals binding energy of the formed contact [32]. In collisions at higher impact energies, the outcome is determined by the degree of inelasticity, which can be a function of temperature [39,40].
The hit-and-stick process can be accompanied by deformation and compaction [37]. However, in the ultra-low velocity regime, compaction does not occur, and the growing agglomerates develop a fractal structure [32,41,42,43,44], with a fractal dimension of D f 1.9 [32,42,45,35]. For higher collision speeds, the fractal dimension may increase to D f ≈ 2.5 [35] or even to the limit of D f ≈ 3.0. High-speed collisions can lead to fragmentation of the agglomerates [46,47,48,49,50]. The mass ratio of the largest fragment to the initial mass decreases with increasing collision speed [50]. In the transition regime between sticking and fragmentation of similar-sized aggregates, bouncing occurs [46,51,52,53,54,55]. Bouncing collisions are inelastic and result in compaction. Successive bouncing events result in a volume-filling factor of the packing of ∼0.36 [56]. We refer to these particles as porous agglomerates.
When a small aggregate (projectile) hits a larger one (target), along with the outcomes mentioned above, mass transfer, cratering and erosion are also possible. Mass can be transferred from the projectile to the target, while the projectile fragments [57,58,59,37,60,47,61,62,50]. The mass-transfer efficiency reaches from a few percent to roughly 50 percent of the projectile's mass, and the mass-transfer process leads to compaction of the target, with typical volume-filling factors of 0.3 to 0.4 [60,61]. For smaller size ratios between projectile and target, cratering occurs instead of net mass transfer. In this case, more mass is lost in the form of ejecta than transferred to the target [63,64]. The relative mass loss of the target agglomerate mainly depends on the impact energy [50] and material strength, and can be up to 35 times the projectile mass [63,64]. For very small projectile-to-target size ratios, erosion (mass loss without a visible crater) has been experimentally observed for silica dust [48,38] and water ice [40], and has also been numerically studied [65,66]. The efficiency of erosion depends on the impact velocity and the projectile mass [38].
In general, all mentioned processes depend on material and disc properties, e.g., monomer-grain material and size distribution, PPD model, turbulence strength and distance to the central star. Different materials can differ significantly in properties, such as sticking threshold or tensile strength, as shown, e.g., for organic materials [67]. This was also observed for collisions of "interstellar organic matter analogues", which showed enhanced stickiness around 250 K, but less sticking for lower or higher temperatures [68,69]. For CO 2 , the sticking threshold seems to be comparable to that of silica [70]. Recently, Arakawa et al. 2021 [71] reconciled the different experimental results on the stickingbouncing threshold for water ice, carbon-dioxide ice and silica particles. Besides the surface energy (closely related to the tensile strength) of the material, the visco-elastic dissipation of energy is the essential material property, and is much larger for water ice at low temperature than for the other two materials. Hence, water ice sticks at higher velocities than carbon dioxide or silica particles of the same size, in agreement with laboratory results [40].
Several barriers prevent the further collisional growth of dust agglomerates to kilometre-sized planetesimals. The most important ones will be briefly described below.
First, drift of particles into the central star, caused by the friction between dust and gas, is a limiting factor. The drift efficiency depends on particle and disc properties, with a drift time scale of Here, a, v d , St and v max are the heliocentric distance, radial drift velocity, Stokes number of the dust particle and maximum deviation of the orbital gas velocity from Keplerian velocity, respectively. Equation (3) is valid for St ≤ 1, and typical values for v max are 50 m/s [31]). This barrier provides a critical time constraint for the formation of planetesimals. The maximum radial drift velocity of 2 v max is reached for bodies with St = 1, which means such particles possess drift timescales on the order of t d ≈ 100 yr for a = 1 au. This timescale is usually much shorter than the collisional growth time beyond St = 1.
Second, the so-called bouncing barrier stops growth when the sticking-bouncing threshold is reached and the colliding aggregates bounce off. Due to the inelasticity of the bouncing process, compaction of the aggregates occurs [56,72,73]. The bouncing barrier is reached at aggregate sizes of millimetres to centimetres, depending on the particle and disc properties [72,73]. The resulting agglomerates have been termed "pebbles". The largest pebbles are achieved in the minimum mass solar nebula model [72]. For increasing heliocentric distance or increasing dust-to-ice ratio, the maximum pebble size decreases and the pebbles become more porous [73]. Millimetre-to centimetre-sized solid particles have been observed in PPDs [74,75,76,77,78,79,80,81,82,83].
Third, even if the bouncing barrier can be overcome, further growth is halted because collisions result in destruction of the aggregates. Fragmentation typically happens at collision speeds of 1 m/s [34].
Fourth, recently, the erosion barrier was experimentally discovered [38]. Erosion is present when dust grains, or small aggregates thereof, impinge larger dust aggregates at high speeds. In this case, impacts liberate other grains from the aggregate surfaces so that the mass of the aggregates is slightly reduced. Numerical simulations have shown that this erosion is a runaway effect even if all other other collisions are assumed to result in sticking. Due to the erosion barrier, maximum aggregate sizes are on the order of 0.1 m [38].
Particles in the weakly ionized disc can carry a nonzero negative charge, which also influences the collisional behaviour. Due to the asymmetric distribution of charges, a repulsive force acts between the particles, hindering sticking [84,85,86]. This effect is called the electrostatic barrier. It can already halt growth at sizes of fractal agglomerates [85]. When taking the photoelectric effect into account, a layer of uncharged dust can build up, which might overcome the electrostatic barrier [87]. However, in low turbulence regions of the disc, such as the dead zone, growth can be effectively suppressed [88]. Numerical models have also shown that the charge of the particles influences growth rate, size and compactness, as highly charged particles grow to larger sizes and get more compact [89], which influences the bouncing and fragmentation thresholds. This may lead to a bridge between bouncing barrier and pebble sizes required for streaming instability, as indicated by experimental and numerical work [90] and could therefore also be beneficial for growth. However, the electrical charging of particles and its influence on the growth process need to be studied further to understand the whole picture.
To summarise, the formation of pebbles, millimetre-to centimetre-sized agglomerates of microscopic dust and ice grains, can be understood with current empirical knowledge about the collision behaviour of protoplanetary dust. However, any further growth seems impeded by the presence of a number of barriers, which -at the current stage of research -makes the direct growth of planetesimals impossible. Thus, other processes that lead to the formation of planetesimals are required. In the next Section, we show how this might be possible.

Pre-concentration of Pebbles
Regions of enhanced concentration of pebbles can be formed in eddies, vortices and pressure bumps of turbulent PPDs [91]. Vortices and pressure bumps can be caused by several phenomena, e.g., magneto-rotational instability [92], 5 FORMATION OF COMETS baroclinic instabilities [93,94], tidal forces of early formed planets and gaps induced by them [95], the edges of the dead zone [96] and snow lines (e.g., [97]). For example, in the transitional disc IRS 48, large-scale vortex structures with a concentration of millimetre-sized pebbles were observed [98].
Turbulent rotating structures can be described by their spatial dimension l, their rotation velocity v e and their turnover time t e = l/v e . For example, particles with stopping times of τ f ≈ t e are efficiently concentrated in high-pressure regions between eddies, which results in a narrow size range for such concentrated particles. Coriolis force and shear are dominant on larger timescales. In such cases, particles with τ f ∼ (2Ω) −1 t e are efficiently trapped, independent of the eddy turnover time. When also taking Keplerian shear into account, the optimal stopping time is τ f = Ω −1 . The radial drift speed can then be deduced from the eddy velocity, and reads Here, the local value of ∆v is important. For example, if ∆v = 0, i.e., in pressure bumps, the particles are being concentrated [91].

Further Concentration of Pebbles by the Streaming Instability
Due to the concentration of particles as described above, the collective surface-to-mass ratio of a pebble cloud is reduced and consequently also the gas drag. If the local dust-to-gas mass ratio exceeds unity, streaming instability is triggered and concentrates the pebbles further [99]. The velocity of the pebble cloud reaches Keplerian velocity, and radial drift is reduced. As a back-reaction by the pebble ensemble, the gas is forced to move at Keplerian velocity. Individual pebbles or small pebble clusters on the same orbit, which still move at reduced speed due to gas friction, are overtaken by the pebble cloud. Inward-drifting pebbles are also incorporated into the pebble ensemble when crossing the respective orbit. Thus, the mass of the pebble cloud rapidly grows over a short timescale. An upper limit of the growth timescale was found to be ∼10 5 years [100]. The length scale of streaming instability is typically 5 percent of the gas scale height [101], and a dust concentration of up to several thousands can be reached [102,103]. The threshold when streaming instability is triggered can also be expressed by the metallicity Z = peb / g , i.e., the surface density ratio of pebbles and gas. For Stokes numbers St ≈ 0.1 (which corresponds to pebble sizes of ∼0.1 m at 1 AU and ∼1 mm at 100 AU), a minimum metallicity of Z min ≈ 0.015 is required, and Z min increases for larger and smaller Stokes numbers [99,104]. As the solar metallicity is Z = 0.0134 [105], enhancement in Z by partial dissipation of nebula gas might be required (see discussion in [106]). However, recent results show the possibility for particle clumping due to streaming instability for smaller metallicities [107]. Substructures of the concentrated regions have been observed in computer simulations [103]. The resilience of streaming instability was investigated by Carrera et al. 2021 [108], who found it to be very robust for centimetre-sized pebbles and in the case of only weak pressure bumps.
Most studies assume monodispersed pebble sizes. However, this is certainly an academic assumption. Recent works on streaming instability with polydispersed pebble sizes show that the size-frequency distribution has an important influence on, e.g., the parameters for instability and growth time scales [109,110,111,112,113]. Additionally, including turbulence reduces the growth rates compared to laminar discs [114] and increases the critical metallicity, depending on Stokes number, needed to trigger streaming instability [104].

Gravitational Collapse
If the mass concentration of the pebble cloud reaches the Roche density with gravitational constant G, gravitational collapse occurs, as described by Johansen et al. 2014 [91].
The cloud collapses into a gravitationally bound body. Whether the pebbles survive the collapse depends on the mass of the pebble cloud. Pebbles can be destroyed either during collisions in the free-fall phase or upon impact on the surface of the forming planetesimal [50,115]. Moreover, the hydrostatic pressures inside the planetesimal can also destroy the pebbles [116]. Thus, the integrity of the pebbles is only guaranteed for low-mass planetesimals. Estimates on the strength of pebbles have shown that this criterion is satisfied for planetesimals with radii 50 km. Size sorting during the collapse could lead to stratification of the pebble size inside a planetesimal [117]. It has been argued that during the collapse, the pebbles may capture smaller dust/ice particles, including those still in the fractal growth stage [118]. If the pebbles survive the collapse intact, the pore spaces between the pebbles may not be entirely empty.
The radii of planetesimals formed by streaming instability with gravitational collapse are predicted to be between 50 and 1000 km, depending on the place of birth in the PPD and the disc properties [119,120,121,103,122,123,124,125].

6
With increasing disc mass, the formed planetesimal sizes increase [103]. The resulting mass-frequency distribution function of planetesimals can be fitted with a power law with exponent −1.6 ± 0.1 [126,127,128], but it was also found that an exponentially truncated power law with a low-mass slope of −1.3 fits better [129]. Klahr and Schreiber 2020 [125] discuss that a Gaussian distribution of the initial planetesimal size could provide a reasonable fit. High variability of the resulting planetesimal size distribution for simulations with slightly different initial conditions was found by Rucska et al. 2020 [130].
However, due to restrictions of the numerical simulation resolution, the lower threshold of the planetesimal size is unknown. With increasing numerical resolution, ever smaller pebble clouds and their collapse can be resolved [103]. Therefore, it is an open question as to which sizes of planetesimals can be formed by streaming instability and the subsequent gravitational collapse. Klahr and Schreiber 2020 [125] suggest that large planetesimals form first and closer to the Sun, and the formation of smaller planetesimals is possible further out in the PPD. Consequently, these objects would form later in time.
It seems likely that pebble clouds can also collapse into binaries, which can merge into contact binaries via a low-velocity impact [131,132]. The influence of the water-ice snowline opens the possibility for a bifurcation in planetesimals with different properties inside and outside the snowline [133].
As stated above, planetesimals larger than 50 km in radius are mechanically altered under their own weight. Pebbles are compacted and destroyed through local hydrostatic pressure, depending on their radial position inside and the total mass of the planetesimal. The loss of the pebble structure is accompanied by an increase of the tensile strength from ∼1 Pa [134] to ∼1-10 kPa [44] due to considerable increase in grain-grain contacts when the pebbles are compressed. Thus, cometary activity by gas pressure cannot be explained without pebbles, because pressures in excess of a few Pa cannot be reached in the cometary subsurface regions [135]. Because we focus on comet formation in this review, we further exclude planetesimals larger than 50 km in radius.

Evolution of Planetesimals towards Comets
In this Section, we discuss the main evolutionary processes acting on the bulk of the planetesimal volume over a timespan of 4.6 Gyrs. This encompasses, on the one hand, internal heating from radioactive decay of short-lived and long-lived radionuclei, and, on the other hand, possible collisional encounters of the planetesimal with other bodies in the outer regions of the solar system.

Radiogenic Heating
Due to the decay of radioactive isotopes, planetesimals can be heated in their interior. This alters the body dramatically if critical temperatures are reached. In this case, sublimation of volatiles depletes ices in the central regions, where temperatures are typically higher than those closer to the surface. Consequently, the advecting volatiles recondense in cooler regions and thus accumulate in shells in the outer parts of the planetesimals. At higher temperatures and pressures, some materials can also melt. Both melting and vaporisation lead to differentiation of the objects. The efficiency of radiogenic heating mainly depends on the size of the planetesimal, the thermal conductivity of the material and on the amount and half-lives of the radiogenic isotopes incorporated into the body.
The radioactive materials are divided into two groups, short-lived radionuclides (SLNs) and long-lived radionuclides (LLNs). SLNs (see Section 4.1.1) possess half-lives typically on the order of one to a few million years, which means that their heating effectively occurs only early in the lifetime of the solar system in the PPD phase. Moreover, as SLNs decay quickly, the formation time of the planetesimals also matters. In contrast, LLNs decay over billions of years (see Section 4.1.2).
Along with thermal alteration due to radiogenic heating, solar illumination can also increase the temperature, depending on long-term orbital dynamics. As shown by modelling of Jupiter-family comets on their path inwards from the scattered disc or the Kuiper Belt [136], this external thermal alteration can penetrate deep into the body, so that hypervolatiles may be redistributed and lost completely. Further, temperatures above 80 or even 110 K can be reached; hence, super-volatiles and amorphous water ice can be affected, too.

Short-lived Radionuclides
The effect of SLNs is dominated by two nuclei, namely 26 Al and 60 Fe. It seems that the distribution of 26 Al was homogeneous in the solar nebula with a ratio of 26 Al(t = 0)/ 27 Al = 5 × 10 −5 [137], with t = 0 denoting the time of the formation of calcium-aluminium-rich inclusions (CAIs). Due to their short half-lives, SLNs can heat bodies to high temperatures. Even the melting point of water ice can be reached in planetesimal-sized objects [138]. In the last two decades, more and more complex studies have been performed to understand the influence of radiogenic heating on the evolution of planetesimals [139,140,141,142,143,144,145,146,147,133]. Most of these studies neglected LLN and did not consider pebble-structured objects. The sizes of bodies investigated in their thermal evolution vary from a few kilometres to several 100 km in radius. In addition, formation time has a significant influence on the evolution of the objects; for example, Lichtenberg et al. 2021 [133] focused on early formed bodies with formation times < 3 Myrs and found mean temperatures of ∼150 K, even for bodies as small as 1 km in radius. Mousis et al. 2017 [146] found for planetesimals with radii of 1.3 km and 35 km that they need to be formed later than 2.5 Myrs and 5.5 Myrs after CAI, respectively, to prevent amorphous water ice from crystallising. Golabek and Jutzi 2021 [147] include, along with radiogenic heating, heating by collisions, and provide estimates for the maximum radius and formation time of planetesimals to avoid heating above 40 K in the interior. The maximum radius can be up to 40 km for a formation later than 5 Myrs after CAI, and even 100 km, including catastrophic collisions for a critical temperature of 80 K (see Figure 9 in Golabek and Jutzi 2021 [147]).
However, the thermal conductivity of the body, which is linked to the internal structure and composition, is of major importance for modelling radiogenic heating. In cases of high thermal conductivity, the heat can be transported efficiently to the surface, where it can be radiated away, resulting in lower peak temperatures, whereas low thermal conductivities impede energy transport from the interior towards the outer regions. Consequently, the peak central temperature can be very high.
It is the pebble-pile structure that makes planetesimals formed by gentle gravitational collapse so interesting for the study of their thermal evolution, because (i) pebbles possess low thermal conductivity due to their porous structure, (ii) the thermal contact between pebbles is minimal, and (iii) radiative heat transport is very inefficient at low temperatures.
The heat-transport model of Gundlach et al. 2012 [148] shows that at low temperatures and for pebbles with radii of 5 mm, the thermal conductivity is less than 10 −3 W m −1 K −1 (see also Figure 1 in [149]) and, thus, is more than three orders of magnitude lower than that of solid (i.e., non-porous, non-hierarchical) materials with the same bulk composition.
The study by Malamud et al. 2022 [150] has for the first time included the pebble structure of small planetesimals (see Section 3.3) in their analysis of radiogenic heating by SLNs and LLNs. They found that even for planetesimals as small as 0.5 to 20 km in radius, the influence of radioactive heating can be significant. Only bodies up to 2 km in radius with a late formation time of >5 Myrs after CAI obtain a peak temperature below 70 K. Higher temperatures force the (super-)volatiles to diffuse outwards and to condense in regions where the temperatures are below the condensation temperature. Differentiation of the ices, with a depleted interior and an enriched outer layer, is the result of this process. However, the near-surface regions are unaltered because outward diffusion stops below the surface. This differentiation of the volatiles is visualised in Figure 1 with red (pebbles depleted in ices) and white (ice-enriched shell) colours. When water ice is redistributed within the planetesimal, no hyper-or super-volatiles can remain in the water-ice depleted interior. Reaching the melting point of water ice was observed in the model when planetesimal radii are 5 km or larger; these planetesimals were formed rapidly after CAI and have large mineral abundance.
The evolved planetesimals after the phase of radiogenic heating by SLNs can be divided into undifferentiated and differentiated objects. Small planetesimals that avoided excessive heating remain undifferentiated and are therefore considered pristine (Type A). The differentiated planetesimals (Type B) are split into subgroups. The first stage of differentiation is the redistribution of super-volatiles, as observed in the simulations of Malamud et al. 2022 [150]. The most crucial differentiation occurs when melting temperatures and pressures (e.g., of water) are reached. In this case, aqueous alteration starts and pebble structures are destroyed, which also decreases porosity. In general, radiogenic heating by SLNs seems to be an important evolutionary process for most planetesimals, except for the smallest ones, which were formed late after CAI. Within the timescale of the first radiogenic heating phase of up to 10 Myrs after CAI, the influence of collisions can often be neglected.

Long-lived Radionuclides
In contrast to the heating by SLNs, LLNs decay over a much longer period of time. Relevant LLNs for planetesimals are 235 U, 238 U, 40 K and 232 Th, and these materials are important for larger objects and for long time scales. Their influence was first analysed by Whipple 1965 [151] and Prialnik et al. 1987 [152]. Both publications found that LLNs can increase the central temperature by a few 10 K. The second phase of radiogenic heating is accompanied by collisional evolution (see Section 4.2).

Collisional Evolution
Besides internal heating, collisions among planetesimals can have a significant influence on their evolution. Collisions can occur over the whole time-range after formation until today. However, collision probability evolves over time and 8 FORMATION OF COMETS depends on the planetesimal size and orbit. Hence, bodies with different sizes and thermal evolutions may collide with one another. The outcome of these collisions depends on the physical properties of the planetesimals, the collision angle, the impact velocity and the size of the colliding bodies. In this review, we only consider collisions that have potentially global impact on the planetesimals, and define three possible collisional types, defined below.
1. First, a collision can be sub-catastrophic if the largest fragment possesses at least half of the initial mass of the larger bodies. We assume that in this case, the mechanical bulk properties of the planetesimals are essentially conserved. 2. Second, if the fragments are smaller, the collision is considered catastrophic. However, if the energy input into the colliding bodies was not sufficient to destroy the pebble structure, the fragments are still composed of pebbles, but the original bulk properties might have changed. 3. Third, the collisions can be super-catastrophic, in which case the pebble structure and the pebbles themselves are destroyed by the collision event.
Several researchers have performed simulations of the collisions of planetesimals for different sizes and impact velocities. Abedin et al. 2021 [153] investigated the collision probabilities in the trans-Neptunian region of today and found a wide range of possible radii and velocities of colliding bodies, but favouring low relative speeds (< 1 km/s) and the main classical belt (30 AU < R h < 50 AU) as the collision site. The scattered and detached populations of the trans-Neptunian region show the smallest collision probability. The physical parameters are connected to those introduced in Section 2 of the dust and ice particles, and the pebbles thereof, and supported by experimental work. However, for planetesimal-sized objects, experimental results need to be extrapolated from laboratory scales. In addition, the impact angle can have an important influence on the induced heating [154]. However, impact heating is most relevant for impact velocities of several kilometres per second and can cause melting of the material. It is interesting to note that porous bodies are heated more efficiently by collisions [155] than solid objects. The numerical study of Charnoz et al. 2007 [156] found that the formation of the Kuiper belt, the scattered disc and the Oort cloud can be explained by dynamical depletion of those regions or low-efficiency implantation of bodies, in both cases accompanied by little collisional activity.
In the following subsections, we briefly review the previous work of sub-catastrophic, catastrophic and super-catastrophic collisions in the framework of comet formation.

Sub-Catastrophic Collisions
We term the cases in which the largest fragment after a collision has at least half of the initial mass of the larger of the two colliding bodies as sub-catastrophic collisions. These include cratering, erosion and collisions resulting in mass growth, such as merging and mass transfer. Golabek and Jutzi 2021 [147] studied the combined effect of radiogenic heating and collisions and found that only the unbound material is heated significantly, whereas the material that remains bound is hardly effected. A study of Jutzi and Asphaug 2015 [157] showed the possibility of the formation of a bilobate shape by a low-velocity collision, assuming non or small cohesion (≤100 Pa) and a collision speed of 1.5 m/s. Furthermore, the formation of bilobed objects by sub-catastrophic collisions was shown by Jutzi and Benz 2017 [158], who derived that such an evolution is more likely to happen than collision-avoiding or a catastrophic breakup of planetesimals. In the following, we assume that sub-catastrophic collisions do not change the morphology, micro-and macro-porosity, mechanical properties and volatile abundances of the largest surviving fragment. Thus, this fragment has the exact same properties as the planetesimal before the collision, with the exception of mass and shape.

Catastrophic Collisions
If the largest fragment after a collision possesses less than half of the initial mass, we term the collision as catastrophic.
The difference from a super-catastrophic collision (see Section 4.2.3) is that catastrophic collisions are not energetic enough to compact or destroy pebbles or to significantly reduce the amount of volatiles. The formation of comet 67P from a catastrophic collision and the following re-accumulation of its current mass, with conservation of volatiles and low bulk density, was proposed by Schwartz et al. 2018 [159]. Jutzi and Michel 2020 [160] showed that the material ejected and not re-accreted is indeed the most processed part of the original object. However, the available simulations do not take a pebble structure of the colliding planetesimals into account. As seen for radiogenic heating (see Section 4.1), the microphysical properties of the material can have a huge impact on the collisional outcome (see, e.g., [161]). In addition, it should also be studied how differentiated bodies behave in collisions. With all these uncertainties in mind, we assume here that bodies that re-accumulated after a catastrophic impact still consist of pebbles and contain a large percentage of the volatiles contained in the planetesimal predecessors. In contrast to the objects stemming from sub-catastrophic collisions (see Section 4.2.1), bodies arising from catastrophic collisions are rubble piles and may consist of spatially separated parts stemming from different regions of the precursor planetesimals. We also assume that this rubble is so mechanically weak [161] that its outer shape does not survive the gravitational re-accretion process. Thus, re-accumulated bodies from sub-catastrophic collisions should not possess voids on length scales exceeding that of the pebble size.

Super-Catastrophic Collisions
The term super-catastrophic collisions describes impacts that provide enough energy to destroy the pebble structure and to sublimate the volatile constituents. In this case, no pebbles or ices will remain inside the object, which is indicated in Figure 1 by grey colour. Super-catastrophic collisions occur at high relative velocities. Objects that gravitationally re-accumulate after a super-catastrophic collision among planetesimals are rubble piles. In contrast to the bodies stemming from catastrophic collisions (see Section 4.2.2), the rubble from super-catastrophic collisions is highly compacted, and thus mechanically strong enough to survive the accretion process intact. Thus, rubble piles from super-catastrophic collisions possess macro-porosity on length scales characteristic for the collisional outcome, and reduced microporosity due to compaction of the dusty material upon impact.

Discussion -Which Planetesimal Can Become a Comet?
With the limited knowledge about actual collisional encounters and outcomes (see Section 4.2), we try in the following to physically characterise a planetesimal after 4.6 Gyrs of thermal and collisional evolution. The formation and evolution scenarios described above may end in a variety of bodies, which we divide into three categories (see Figure 1 for an overview): • Icy Pebble Piles-Type A1 and Type B1 • Icy Rubble/Pebble Piles-Type A2 and Type B2 • Non-Icy Rubble Piles-Type A/B3 One very important question for researching comets is: Which of these types could represent the comets currently observed in the solar system? In the following, the most important physical properties, which can be used for the validation of the formation/evolution model and the corresponding observables, will be discussed. The expectation from the formation and evolution scenario for each observable will be described and compared to real physical properties of comets measured by (spacecraft) observations.

Intrinsic Physical Properties
We will start this discussion with those physical properties that are directly provided by the formation and evolution process. These properties, such as internal morphology and total porosity, are referred to as intrinsic physical properties and can be used to directly link the formation and evolution of planetesimals with the observable properties of contemporary comets in the solar system. Table 1 summarises the different intrinsic properties of cometary nuclei.

Internal Morphology and Pebble Properties Expectation
Planetesimals formed by the gentle collapse of a pebble cloud are initially very homogeneous throughout the whole body, despite the possibility of size-sorting of the pebbles while the cloud collapses [117]. The latter would generate a gradient in the size of the void spaces between the pebbles, but not necessarily in density. In any case, the size of the void spaces between the pebbles is always on the order of the pebble size. As described above, compaction due to lithostatic compression of large planetesimals would decrease the void space in the inner parts of the objects. However, for simplicity we excluded these bodies (i.e., >50 km in size) due to their large sizes compared to typical comets. As pebble size is expected to be on the order of millimetres to decimetres, density contrast on length scales much larger than decimetres is not expected in the interior. For type A1, no alteration of the internal structure due to evolutionary processes is expected, and the internal morphology remains pristine, with pebbles and void spaces of pebble size that are filled with fractals (see Section 3.3). For type B1, the internal morphology is almost identical to A1, with the exception that in the deep interior, the porosity is increased due to the evaporation of ices, and in the ice shell(s), the pore size is reduced due to recondensation of the ice(s). If a catastrophic collision occurred, fractals may have been lost where the void spaces were opened. Therefore, fractals may be locally depleted for type B1 and B2, but the internal morphology of pebbles and void spaces remains unaltered because catastrophic collisions do not destroy the pebble structure.
From the modelling point of view, pebble sizes are determined by the bouncing barrier [37,72]. For pebbles containing microscopic water-ice particles, Lorek et al. 2018 [73] derived maximum pebble sizes on the order of 1 mm to 1 cm, With a super-catastrophic collision (type A/B3), pebbles and fractals are destroyed, and the resulting rubble piles are highly compacted, and thus relatively rigid. Within this rubble, large void spaces can remain. The size of the rubble may vary and is not further constrained in this work.

Comparison to Observations
The internal morphology of comets is only measurable by instruments with the ability to scan the interior of the objects. The Rosetta Mission was equipped with instrument packages capable of performing these investigations, namely Comet Nucleus Sounding Experiment by Radiowave Transmission (CONSERT [165]) and Rosetta Radio Science Investigation (RSI [166]) instruments. Both devices performed radar observations of the interior of the nucleus.
Their measurements showed that the observed volume of comet 67P is highly homogeneous on length-scales larger than 3.3 m (CONSERT [167]). However, CONSERT could not investigate the entire body, so the results are limited to a fraction of the volume of comet 67P. CONSERT found a denser surface layer for the uppermost < 25 m [168]. Such variation may be explained by the redistribution of volatiles due to radiogenic heating (see Section 4.1). RSI measurements showed an offset of the centre of gravity, which indicates density difference of the two lobes. The big lobe seemed to be slightly denser than the small lobe [169].
Direct observations of the pebbles have not been unambiguously performed. Imaging of the surface of comet 67P with the highest spatial resolution was performed by the CIVA instrument onboard the Philae lander [170] and showed substructures with sizes of a few millimetres to one centimetre, which might be identified as pebbles. Other, indirect methods have been applied by Blum et al. 2017 [3] to derive the pebble size and result in values for the pebble radii between 3 mm and 6 mm. Cometary trails have been investigated by Lisse et al. 1998 [171], who found that a high abundance of macroscopic particles with β values of β < 10 −3 , corresponding to dust-particle sizes of 1 mm, are present in the trails. Moreover, cometary meteor showers exhibit particle sizes up to ∼1 cm [172]. As stated above, direct observational evidence of pebble substructures is not available.

Bulk Porosity Expectation
The bulk porosity of planetesimals is mainly influenced by formation processes and can then be altered due to evolutionary processes. As a simple rule-of-thumb, the packing density of a planetesimal consisting of two hierarchy levels (dust/ice particles forming pebbles; pebbles forming planetesimals) is Φ bulk = Φ d Φ p ≈ 0.24, with Φ d ≈ 0.4 and Φ p ≈ 0.6 being the packing densities of the dust/ice particles inside a pebble and of the pebble packing inside the planetesimal, respectively [173]. For pebbles with substructures, the internal packing might be Φ d < 0.4, as this can be regarded as another hierarchy level [2]. Thus, the bulk porosity adds up to ∼76 %. For larger planetesimals, the porosity is expected to decrease with increasing depth due to lithostatic compression of the pebble structure [150]. However, as described in Section 4.1, radiogenic heating can affect even rather small bodies when they form early enough (type B1). The outward diffusion of hyper-and super-volatiles goes along with a change in porosity so that the porosity increases above the primordial value close to the centre and decreases at locations where the volatiles condense. However, when no volatiles are lost into space, the bulk porosity does not change and remains pristine. For larger objects, a change of porosity can also be induced due to the melting of water ice. When water ice melts, the pebbles collapse and immerse into the water so the initial porosity is lost, and consequently, the body shrinks. However, as this only happens close to the centre, the outer layers maintain their pebbles and their porous structure. For type B2 bodies, the pebbles from the various regions possess different dust-to-ice ratios, which implies that the observed density across the surface and for different depths changes. For the compacted structures of type A/B3, very little porosity remains inside the rubble, but macroporosity may exist within it. Thus, the bulk porosity of the entire body of a type A/B3 should be lower than the pristine value.

Comparison to Observations
The porosity of a whole body can be derived from its volume, mass and material density. Material composition and porosity can be estimated by measurement of the electric permittivity [174], but this procedure is model-dependent and relies on calibration experiments [175]. Measurements of total porosity with high resolution were only performed for comet 67P, resulting in values for the bulk packing density of Φ bulk = 0.15 to Φ bulk = 0.37. The RSI instrument measured values between Φ bulk = 0.25 and Φ bulk = 0.30 [176]. An estimate of Φ bulk = 0.21 to Φ bulk = 0.37 was retrieved by directly measuring the pebble density by Grain Impact Analyser and Dust Accumulator (GIADA [177,178]). Moreover, permittivity measurements done by CONSERT provided a range of Φ bulk = 0.15 to Φ bulk = 0.25 [167]. The lower value of Φ bulk = 0.15 leaves room for a third hierarchy level, i.e., for pebbles with substructures, whereas the upper value requires the pebbles to be homogeneous.
Due to the overall high porosity of comet 67P, evolutionary processes that considerably reduce the porosity, such as, e.g., super-catastrophic collisions, can be excluded.

Dust-to-ice Mass Ratio Expectation
The dust-to-ice mass ratio (DIMR) of a planetesimal is given by the material distribution inside the pebble cloud that undergoes the gravitational collapse. The DIMR depends on the location of formation as well as on the preceding mixing processes of the materials during pebble growth. We refer to the DIMR at cloud collapse as pristine, but it should be noted that this DIMR is not necessarily equal to the DIMR of the solar nebula or of the bulk protoplanetary disc outside the various snowlines. During the thermal and collisional evolution of the planetesimals, the DIMR can be altered.
For comet types A1 and A2, the pristine DIMR is globally and locally preserved. Throughout the entire volume, the distribution of ices is homogeneous down to the pebble scale. For type B1, the bulk abundances of all ices are preserved, but the local DIMR is changed. The interior of the planetesimal is depleted in ices (higher DIMR), and the ice layers are enriched in ices (lower DIMR). The near-surface layer remains mostly pristine in the DIMR, despite locations where material was ablated by collisions. As type B2 is the breakup product of a differentiated body, the distribution of ices is inhomogeneous and depends on the depth of the place inside the body where the rubble was originally stored. The manner in which such a body re-accumulates the collisional fragments determines the distribution of the ices [150]. Whether the DIMR of the rubble is changed compared to the pristine value depends on the details of the preceding catastrophic collision event. For types B1 and B2, super-volatiles could be lost entirely if the corresponding threshold temperature was reached at the moment of the collision. Type A/B3 does not contain any ices due to the hyper-velocity nature of the collision event.

Comparison to Observations
The dust-to-ice mass ratio F can also be addressed by porosity and density estimations, as shown by Lorek et al. 2016 [179] and Patzold et al. 2019 [180]. For cometary refractory material, the material density is not known and can only be constrained within a range of reasonable values. Measurements by Rosetta found an organic-to-mineral mass ratio of roughly 1:1 [181]. Typical densities of organic matter (1000-2000 kg/m 3 ), combined with typical densities of minerals (2600-5000 kg/m 3 ), result in an expected density range of 1800-3500 kg/m 3 for the refractory material of a comet. GIADA measured comparable densities with 1925 +2030 −560 kg/m 3 [182]. Some particles showed a bulk density >4000 kg/m 3 . For water and CO 2 ice, the densities are ρ H2O = 934 kg/m 3 and ρ CO2 = 1600 kg/m 3 , respectively. Porosity estimations for comet 67P are shown in Section 5.1.2. With this, one can calculate the dust density ρ d with: with f d + f H2O + f CO2 = 1 the mass fractions of dust, water and CO 2 , respectively, ρ n = 532 kg/m 3 the nucleus density [169] and bulk porosity Φ. The results for varying DIMR are shown in Figure 2. From this plot, one can constrain the global porosity and the ice content. Porosities larger than ∼81% cannot be explained with reasonable dust densities and should therefore be excluded. For a porosity of 76% (see Section 5.1.2), the DIMR has a range of 3 to 4 depending on the assumed CO 2 fraction. For the smallest porosity of 63% derived from Rosetta instruments, only a small range of DIMR values between 0.8 and 2.5 is reasonable.
These values are in agreement with the review written by Chourkoun et al. 2020 [183] about the refractory-to-ice(s) mass ratio of comet 67P. In this work, the authors showed, based on different Rosetta observations, that the DIMR is most probably higher than 3 (see Figure 5b in [183]). An upper limit of the DIMR, which is in agreement with most of the Rosetta measurements, can be between 6 to 8. Only two papers argue for much higher maximum values. The works claiming higher DIMR values are based on coma dust grain density estimation around perihelion [182] and on the analysis of the fallback of volatile-bearing chunks in the northern hemisphere [184]. Figure 2: This plot shows the resulting dust density for varying dust-to-ice mass ratios and porosities. The greencoloured area indicates reasonable dust densities. The given CO 2 percentage corresponds to the CO 2 content of the total ice mass.

Binarity, Flattening and Rotational Orientation Expectation
Numerical simulations have shown that gravitationally unstable pebble clouds can collapse into binary planetesimals. These objects rotate in a prograde manner in 80% of the simulations [132]. This finding is almost not affected by different parameters influencing streaming instability. Such a formation of binaries results in very similar compositions of both bodies because they form from the same reservoir of material. The radii of both bodies would also be very similar [185]. The spinning of the planetesimals induced by the collapse results in a flattened shape [186]. Additionally, the parts of such a binary align over time because this reduces the energy of the system. The alignment can be followed by a reduction of the distance between both bodies until they contact due to, e.g., Kozai-Lidov cycling, the YORP and BYORP effects, tides, collisions or gas drag [186]. However, the collapse of a pebble cloud can lead to a wide variety of binary distances, including contacting binaries [131]. The survival of binaries in collisions was addressed by Nesvorny et al. 2019 [187], and they found that more tightly bound binaries have a higher survival probability than wider ones. Regarding radiogenic heating, due to the small contact area compared to the total area, the influence is negligible and should be equal for each body individually. It should be mentioned that the binarity of low-density objects also has important consequences for the escape of dust particles when dust activity is present (see, e.g., [188,189]).

Comparison to Observations
Many of the observed cometary nuclei show a binary shape. The most detailed observed one, comet 67P, exhibits this in great detail. However, the formation of binary shapes is still discussed and researched. As shown before, the formation of binaries is an intrinsic characteristic of formation by a collapsing pebble cloud. Nesvorny et al. 2019 [132] showed that for binaries that have formed by capturing another object, the distribution of inclinations would favour retrograde over prograde rotations. Comparison to the known binaries in the Kuiper belt shows that the distribution of binary inclinations can be best explained by direct formation of binaries from a single source [132]. Arrokoth shows alignment of the two lobes, which is in agreement with the predicted scenario [186]. However, for comet 67P, the observation of layered surface structures surrounding each lobe [190,191,192] suggests that the lobes formed individually and contacted later.

Cometary Activity Expectation
If we define cometary dust and gas activity as a state in which insolation produces subsurface pressures large enough to overcome the cohesion of the dust, then types A1, A2, B1 and B2 can in principle become active due to the presence of volatile ices close to the surface. However, these objects show significant differences, which are noteworthy and could become decisive to answer the question of how comets formed and evolved. The distribution of volatiles over the entire surface of the bodies is different. Due to the homogeneity of types A1 and A2, these objects should become active at each surface location given sufficient insolation. When the surface of a type-B1 object gets eroded locally, the activity can be different at this location because an ice shell or the depleted interior could be excavated. We could observe this effect by local variations of the outgassing rates of different volatile species. For type B2, the different ice-containing rubble is distributed inhomogeneously over the whole body. Hence, also for type B2 we expect locally varying outgassing rates. Type A/B3 cannot show cometary activity driven by ice sublimation as it no longer contains any ices.

Comparison to Observations
The first-ever observed cometary nucleus was the nucleus of comet 1P/Halley. Halley was expected to show global activity, but when the Giotto spacecraft arrived, only ∼1/3 of the surface was contributing to global activity [193]. A similar observation was made for comet 19P/Borrelly -one distinct jet originated from the cometary nucleus, which actually consisted of three smaller outbursts [194]. In contrast to other comets, comet 81P/Wild 2 possessed a diverse and complex variety of surface structures [195]. Longer-exposure camera images revealed the presence of many jets around the entire nucleus. More than twenty highly collimated jets were observed, which indicated a heterogeneous surface. In addition, analysis of captured particles indicated that comet Wild 2 may have experienced aqueous alteration in the past [196]. Comet 9P/Tempel 1 showed rather variable dust activity, with frequent natural outbursts all over the surface [197]. In comparison to the above-mentioned comets, 103P/Hartley 2 had very distinct activity patterns. During the flyby of the EPOXI spacecraft [198], water vapour emission was detected at the waist of nucleus. Localised CO 2 -active areas also showed ejection of solid water ice chunks of up to cm size. In addition, solid ice particles were also observed in an area where only water was active.
Comet 67P is by far the best-studied comet, and this is why we are aware of several local activity hotspots (such as sunset jets; see, e.g., [199]), although the nucleus seems to be globally active [200]. The activity pattern changed with time because of seasonal variations of solar illumination and because of material that fell back onto the nucleus [184]. During perihelion, the southern hemisphere showed an increase of CO 2 activity that led to the ejection of larger ice-containing chunks (up to several cm in size [201]) and to a bluing of the surface because more water-ice-containing pebbles were excavated by CO 2 erosion (for details see [164]). These chunks either escaped into space or fell back into the northern hemisphere of the nucleus. Although we understand the principle of how gas activity leads to the ejection of pebbles or chunks from the surface, it is still not known how to explain the particle-size distribution in the coma [3]. The smallest (sub-)micrometre-sized particles cannot be lifted directly from the surface because of their strong adhesion forces. This means that the chunks or pebbles must disintegrate in the coma after lift-off to produce the observed fraction of small particles in the coma. Fulle et al. 2020 [202] proposed that pebbles possess a substructure called agglomerates, typically ranging from micrometre to millimetre in size. Probably charging effects, fast rotation of the pebbles or extreme illumination conditions lead to break up of the pebbles into unbound agglomerates with a size distribution similar to the observed one.

Thermal Conductivity
Expectation Heat can be transported by conduction, radiation, gas diffusion and convection in the case of fluid material. The latter mechanism can be neglected in contemporary comets because they have certainly cooled down, even if they had at some point reached the critical temperature for water-ice melting [150]. Further, heat transport by gas diffusion can be neglected in most cases due to the low gas pressures involved. For the conduction and radiation of thermal energy, material and structural properties determine the efficiency of heat transport. In a granular medium, conduction is less effective than in solid material due to the reduced contact areas between neighbouring particles, which is described by the so-called Hertz factor. For pebbles, heat conduction is reduced further because of the granularity inside the pebbles. Therefore, and due to the large void spaces between the pebbles, thermal radiation contributes to the energy transport and can even dominate for large pebbles and high temperatures, as shown in Bischoff et al. 2021 [149] (their Figure 1). Only for low temperatures (<100 K for pebbles with 5 mm radius) does heat conduction dominate, but then the thermal conductivity is extremely small (<10 −3 W/(K m)). For the five types of bodies proposed, this means the following: Types A1 and A2 with their pebbles show low network conduction, and radiation dominates for higher temperatures. For type B1, the shell structure influences the conductivity such that the interior with reduced ice content shows decreased network conductivity compared to the pristine material due to the higher porosity and the low internal temperatures. As the ice shells are enriched with material, their network conductivity is locally increased. For the outer pristine material, the network conduction is the same as for types A1 or A2. For type B2, which possess regions of different dust-to-ice ratios, thermal conductivity will vary among the different locations. Ice-depleted areas, which possess higher porosities, are subject to strong day-night thermal conductivity variations. Areas of higher ice content generally possess higher thermal conductivity because of the network part, but temperature changes do not have the same importance as for the ice-depleted parts. In case of type A/B3, the rubble is expected to be compacted, and therefore, network conduction is the relevant transport process.

Comparison to Observations
Instead of the thermal conductivity λ, often the thermal inertia Γ = √ λρc is measured for small solar system objects, with ρ and c being the mass density and heat capacity of the material, respectively. However, these properties are often not known with high accuracy, and additionally, the temperature dependency of all these values is often not taken into account. As described above, the temperature dependency can be crucial for highly porous pebble-structured materials. Thus, measured values of the thermal conductivity or thermal inertia should be considered with caution and in the corresponding temperature range. For comet 67P, several instruments onboard Rosetta and Philae were able to measure temperatures. On the orbiter, Visible InfraRed and Thermal Imaging Spectrometer (VIRTIS) [203] measured the surface temperature at depths of tens of micrometres. Self-heating and shadowing effects were observed by VIRTIS [204]. However, VIRTIS had a lower limit for temperature detection of ∼156 K due to instrument noise. The Microwave Instrument for the Rosetta Orbiter (MIRO) [205] used wavelengths of 0.5 mm and 1.6 mm, resulting in penetration depths on the order of a few millimetres to centimetres [3]. The MIRO temperature data are difficult to interpret because the measurement depth is comparable to the diurnal skin depth [206,207], and therefore a wide temperature range contributes to the measured brightness temperature. Marshall et al. 2018 [208] found a best-fitting thermal inertia of 80 J/(K m 2 s 0.5 ) under the assumption of a complex roughness distribution when comparing MIRO and VIRTIS data. Onboard the lander Philae, Multipurpose Sensors for Surface and Subsurface Science (MUPUS) [209] was capable of probing the surface with the MUPUS PEN (thermal probe) and the MUPUS TM (16 resistance temperature detectors and an infrared radiometer). They found a best-fitting thermal inertia of 85 J/(K m 2 s 0.5 ) [210]. Assuming a heat capacity of 560 J/(kg K) (as used in Bischoff et al. 2021 [149]) and a density of 532 kg/m 3 (for the bulk nucleus, [169]), these thermal inertia values translate into thermal conductivities of ∼0.02 W/(K m). As described before, the uncertainty in heat capacity and density of the measurement spots is forwarded into uncertainties of the thermal conductivity. Thus, such a value alone cannot be used to conclude the kind of structure of the cometary surface. The method proposed by Bischoff et al. 2021 [149], which allows the distinction between micro-and macro-porous subsurface structures by relating the sunrise temperatures to the insolation at noontime, could not be used with measured data of comets due to the lower limit of available temperature data. However, in future missions to comets, such as ESA's Comet Interceptor [211], the application of this method should be kept in mind when measuring surface temperatures.

Tensile Strength Expectation
The cohesion of granular matter is in general influenced by the material composition, its structure and the observed length scale. The tensile stress a body can withstand before it breaks is called the tensile strength. For non-porous solid material, the tensile strength is highest due to the direct contact of the atoms and is typically found in the MPa range. Tensile strength is reduced by any kind of porosity, granularity or roughness because of the reduction of contact area among the atoms. For objects composed of micrometre-sized grains, i.e., for the inside of pebbles, the tensile strength is on the order of kPa [212,67]. For pebble-piles, the tensile strength is further reduced to a few Pascals or even less, depending on pebble size [173,134,54]. For planetesimals formed from pebbles, such a low tensile strength is expected on the metre-scale or above. Compaction due to collisions or redistribution of volatiles leads to higher tensile strength values. The relation between volume-filling factor and tensile strength has been measured for several materials, such as silica, graphite and organics, and depends strongly on the compaction [212,213,214,67]. Applying this knowledge to the five types of bodies discussed at the bottom of Figure 1 leads to the conclusion that types A1 and A2 possess low tensile strength on the order of 1 Pa down to the length scale of the pebble size. As shown in Section 4.2, sub-catastrophic and catastrophic collisions can occur without causing significant damage to the pebble structure. Hence, the tensile strength is not affected by collisional evolution in these two cases. However, differentiation causes significant changes of the tensile strength. Locations of enhanced ice content will possess higher tensile strength values, whereas depleted volatile material consequently leads to a decrease in tensile strength. This is why the redistribution of material in the rubble in type B2 leads to variation of the tensile strength over the surface and volume of these bodies. The material of type A/B3 bodies is highly compacted, and this is why the rubble should have a high intra-rubble tensile  Figure 3: Response function of a granular material to compressive stress. The green curve shows the reaction of a homogeneous, i.e., non-hierarchical, granular material, whereas the blue curve represents the volume-filling factor of a pebble-pile structure, i.e., a hierarchical material. Both functions can be described by S-shaped curves as shown by [216]. The homogeneous material can only be compressed in a single stage (transition regime I), while for the hierarchical material, first the pebble packing will be compressed (transition regime II), then at higher pressures the pebbles will be deformed and destroyed (transition regime III; dashed curve).
strength. However, the tensile strength between the rubble particles is almost zero, which means that the rubble is only bound by gravity.

Comparison to Observations
As described above, the tensile strength of a body varies for different length scales. For large scales of several metres to 100 m, the tensile strength of comet 67P was estimated by the observation of cliff collapses [215] and found to be on the order of 1 Pa, as expected for a pebble-pile body. For length scales of a kilometre or above, the tensile strength is difficult to measure because of the increased importance of gravitational binding. The central hydrostatic pressure of comet 67P, for example, is on the order of 10 Pa, and thus larger than the expected tensile strength of a loose pebble pile. However, for bodies with radii of several kilometres or above, a memory effect can lead to increased tensile strength as shown by Blum et al. 2014 [134]. These results may indicate why the lower limits of the tensile strengths for a number of comets are slightly above the canonical 1 Pa value [44]. Much smaller length scales can be probed through the breakup of cometary meteors in Earth's atmosphere. These typically millimetre-sized bodies possess tensile strengths on the order of 1-10 kPa [134] and may thus represent intact pebbles or fragments thereof.

Compressive Strength Expectation
Compressive strength describes the ability of a material to withstand pressures tending to reduce size. For granular materials, this material property depends on the volume-filling factor of the material as well as on the grain properties, such as shape and size distribution. This is why the compressive strength of a granular medium cannot be seen as a single value, but should rather be considered a volume-filling-factor-dependent value that changes with respect to the applied compressive stress (see Figure 3). The reason is that fluffier materials can be compacted much easier in comparison to more compact structures. The functional behaviour of the compression curve of granular materials was derived using theoretical considerations and laboratory experiments [216]. The result is that the response of the volume-filling factor to compressive stress can be described by an S-type function (see Figure 3), with a transition regime that is characterised by the so-called turnover point, p m . Homogeneous granular materials are compressed in the transition regime I. Larger pressures cannot lead to any further compaction because the grains are already positioned in the densest possible packing. For pebble piles, however, the compression function can be separated into two different regimes, namely transition regime II, in which the pebbles are brought into a denser packing without destruction [150]; and transition regime III, in which the pebbles cannot withstand the applied pressure and are, hence, destroyed [217].
For type A1 and A2 bodies, low-stress compressive strength should be low and comparable to the tensile strength.
For types B1 and B2, the general behaviour should be similar, but due to the enhanced volume-filling factor at the ice shell, the compressive strength should be increased, and inhomogeneously distributed in type B2. The compact rubble in bodies of type A/B3 possess a high compressive strength on length scales smaller than the rubble size.

Comparison to Observations
The only available estimations of the compressive strength come from the landing attempts of Philae on the nucleus of comet 67P. Heinisch et al. 2019 [218] determined an upper limit of 800 Pa, while O'Rourke et al. 2020 [217] derived an upper limit of only 12 Pa for the compressive strength. Based on the compression model for granular matter under reduced gravity conditions by Schräpler et al. 2015 [216], these low values can be interpreted as being due to loosely packed pebble-sized bodies on the surface of comet 67P [217].

Gas Permeability Expectation
The gas permeability of a small solar system body is determined by its porosity and the characteristics of the pore spaces, such as their size distribution and percolation properties. Tortuosity describes the actual length of the path the gas takes divided by the closest distance between the start and end points. Hence, high tortuosity hinders the gas flow, whereas low values of tortuosity provide favourable gas-flow conditions. For planetesimals, the gas permeability therefore depends on the porosity distribution inside the body. For larger pebbles, the gas permeability increases. Therefore, bodies of types A1 and A2 possess relatively high gas permeability. However, if the redistribution of volatiles due to radiogenic heating results in filling the voids in the ice shell (types B1 and B2), the gas permeability can be greatly reduced, and the flow could even be stopped if the pores are closed completely. In the deep interior, the permeability could be increased where icy pebbles are extinct, but this depends on whether there are pebbles consisting only of ices, leaving large voids behind, or pebbles that are mixtures of ice and dust, with the voids created by ice evaporation being small. For type B2 bodies, the permeability is inhomogeneous and follows the distribution of rubble from the interior, the ice shell or the pristine mantle. As type A/B3 bodies do not contain any ices which could sublimate, gas permeability is not of interest. In general, the gas permeability directly influences the possible pressure build-up inside a body with evaporating volatiles, which can lead to ejection of material if the cohesion is overcome.

Comparison to Observations
Measurements of the gas permeability of the cometary subsurface material are not directly possible. However, using thermophysical models for the interpretation of the gas production rates of comets may shed light into the vapour transport inside the cometary nucleus. This was done by, e.g., Gundlach et al. 2020 [219], who found that a pebble pile can explain the observations of the dust and gas emission of comet 67P around perihelion. Skorov et al. 2021 [220] investigated how variations of porosity and inhomogeneities of the surface material affect the outgassing behaviour. They found high influence due to porosity changes, but only minor contributions due to cavities and cracks. Skorov et al. 2022 [221] showed that heterogeneous microstructures of dust layers can be neglected, so that effectively a homogeneous dust layer is a valid assumption for modelling outgassing. Further modelling of the transport characteristics was performed by Reshetnyk et al. 2021 [222], which is useful for the application to comet observations.

Permittivity Expectation
Permittivity is a measure of the electric polarisability of a material in response to an electric field and is dependent on the frequency of the electromagnetic wave, porosity, temperature and composition of the body. It can be divided into a real and an imaginary part. Brouet et al. 2014 [223] investigated the permittivity of porous granular matter in different grain size ranges. Besides refractory material, water ice can have a large influence due to the strong dependency of the permittivity on frequency. For frequencies below 10 4 Hz, the permittivity is roughly three orders of magnitudes higher than above 10 4 Hz [224]. This effect could be used to investigate the water-ice abundance. Due to the homogeneity of objects of types A1 and A2, the permittivity is homogeneous throughout the whole body. The differentiation in type B1 objects results in a radial gradient of the permittivity. As describe above, water ice and silica can have similar permittivities in some frequency ranges. Therefore, distinction between the two is only possible in frequency regimes in which the materials possess different permittivities. For objects of type B2, the permittivity is distributed inhomogeneously, following the rubble material composition. As A/B3 bodies are by definition ice-free, there is no influence due to water-ice permittivity, but the permittivity is expected to be homogeneous inside single rubble.

Comparison to Observations
With the Comet Nucleus Sounding Experiment by Radiowave Transmission (CON-SERT) [165], the permittivity of comet 67P was investigated with a frequency of 90 MHz [225,174,168]. The relative permittivity ranged from 1.7 to 1.95 in the subsurface (<25 m) and from 1.2 to 1.32 in the interior. A denser surface layer or a difference in composition could result in this dichotomy [168], which could have evolved due to activity and redistribution of volatiles. The interior seems to be more porous. This characteristic is similar to that of a type B1 body, but the length scales could be different, and more detailed simulations are needed to address this point. Further, the Surface Electric Sounding and Acoustic Monitoring Experiment (SESAME) [226] instrument package, including the permittivity probe onboard Philae, found hints of a more compacted surface layer compared to the interior [227]; however, the instrument only probed the first metres. Additionally, MIRO data can be used to address the material properties, as performed by Bürger et al. (subm.) [228] for optical constants in the millimetre-wavelength range, which are linked to the permittivity. They found that a pebble-structured surface model can fit the measured MIRO data. In general, the Rosetta results of permittivity suggest high homogeneity, similar to bodies of type A1 and A2.

Conclusions and Open Questions
In this review, we presented our current understanding of planetesimal formation and evolution in the context of comets being possibly the most primitive survivors of these eras. We illustrated the most crucial processes and pathways from protoplanetary dust to contemporary bodies of the solar system in Figure 1. As ices are an important ingredient of comets, comet formation as macroscopic bodies must start beyond the snowlines of the volatiles and super-volatiles found in comets, e.g., CO 2 and CO. Through three phases of coagulation (see Figure 1 and Section 2), pebbles consisting of microscopic ice and dust grains form in the size range of millimetres to decimetres and cannot grow further due to growth barriers. These pebbles get concentrated by hydrodynamic effects, leading to streaming instability, which enhances the concentration of pebbles to the gravitational-instability limit, followed by a gentle collapse of the pebble cloud into a many-kilometre-sized object. These planetesimals follow a broad size-frequency distribution, but to evolve into comets, their size cannot exceed 50 km due to the lithostatic pressure, which destroys the pebbles for larger planetesimal sizes. Present numerical models of the collapse stage do not have the resolution to predict frequency at the small end of the expected size distribution of planetesimals, particularly for sizes in the 1-10 km size range applicable to comets. Here, more numerical work is needed to investigate the expected properties of the smallest planetesimals.
The first phase of planetesimal evolution is dominated by radiogenic heating by short-lived radioactive nuclei, which can redistribute volatiles inside the planetesimals. This phase is followed by collisional evolution, whose outcome depends on the specific impact energy. In this review, we differentiate between sub-catastrophic, catastrophic and super-catastrophic collisions. All these evolutionary processes result in five distinct types of evolved planetesimals, designated by us as types A1, A2, B1, BS and A/B3 (see the bottom of Figure 1). Additionally, it should be kept in mind that thermal alteration by solar illumination on long time scales regarding the path from the Kuiper Belt or the scattered disc inwards can influence the abundances of volatiles [136].
The Rosetta Mission to comet 67P Churyumov-Gerasimenko provided a huge number of scientific results as a baseline for the study of comet formation and evolution. However, there are still many open questions, and it should be kept in mind that these detailed measurements are only available for one comet whose representation of the class of all comets in general is unknown. Only objects of type A/B3 as non-icy rubble piles can be excluded to represent a contemporary comet. Icy pebble piles (objects of types A1 and B1) preserve most of the pristine planetesimal properties and would be most interesting for the understanding of planetesimal formation. However, the possibilities to distinguish icy pebble piles from icy rubble/pebble piles (types A2 and B2) are restricted. As shown in Tables 1 and 2, the only distinction between types A1 and A2 seems to be the remaining amount of fractal dust aggregates and the binarity. Additionally, the expected differences between the various types need to be modelled numerically, as here we have only made a relatively generic description of their physical properties. Hence, simulations of collisional outcomes of colliding pebble piles are urgently required, including features such as the survival of fractals in the void spaces between the pebbles and differentiation of the bodies following heating episodes.
Measurements of the permittivity of ice and dust pebbles would be helpful to investigate the possibility to differentiate between volatiles and super-volatiles within comets and to compare observations with thermophysical modelling of the radiogenic heating. However, this is a complex endeavour because deep-penetrating radar requires the usage of large wavelengths. DC geo-electric methods with high penetration depths require spacecraft missions with at least two landers.
Further measurements of temperatures at the surface as well as in the upper layers of cometary nuclei are needed to constrain the makeup of the cometary subsurface through comparison with thermophysical models. However, as there are many factors affecting temperature, such as surface roughness, self-heating and shadowing, the models need to be applied with care. In general, the cooling phase in the cometary night seems to be most promising because the above influences are reduced, and models have shown that much information can be extracted.
Future missions to and astronomical observations of comets may be able to test the following most crucial predictions of the pebble-cloud-collapse model of planetesimals and subsequent evolutionary processes into comets: i Strong positive correlation between the surface temperature at sunrise and the insolation at local noon for a subsurface made of pebbles, in contrast to no such dependency for a makeup without large void spaces [149], measurable by thermal IR mapping. Such measurements would deliver, from remote observations only, invaluable information about the presence and size of pebbles in a shallow subsurface layer.
ii Proof of the absence or presence of internal volatile differentiation of comets, as predicted by Malamud et al. 2022 [150] for comets consisting of pebbles, measurable by long-wavelength radar. Such measurements would deliver information about the formation time and/or abundance of radiogenic nuclei and about the thermal conductivity in the deep interior of the comet nucleus, and thus would confirm or disprove the gravitational-collapse theory.
iii Search for traces of the collisional evolution of cometary nuclei, measurable by radar through the amplitude and length scale of internal inhomogeneities. Due to the obvious distinction of the internal makeup of bodies of type A1, B1, A2, B2 and A/B3 (see Figure 1), deep-penetrating radar measurements covering the entire body would provide information about the internal distribution of water ice, the refractory component and the void spaces inside the body.
iv Determination of the physical properties and orbital parameters of extinct comets. Due to memory effect (see [134]) or differentiation (see [150]), the depth to which dust activity is possible strongly depends on the original size of the comet-precursor body, so larger comets should be going extinct on a different timescale than smaller ones (see [229]).
v Determine whether positive relief features [106] are local remnants of impacts, measurable by high-frequency radar through local permittivity enhancement. If moderate-velocity impacts on the surface of a planetesimal lead to local compaction of the material and thus to loss of the pebble structure, further dust activity is impossible due to considerable enhancement of the tensile strength. Such measurements would deliver information about the collision history of the planetesimals and the abundance of small-scale (metre-to decimetre-sized) objects in the region in which the planetesimals resided throughout most of their lifetime before becoming comets.
It should be noted that the picture of formation and evolution of planetesimals drawn in this review cannot be complete because it ignores evolutionary processes that are usually summarised under space weathering. We again emphasise that we only provide one plausible, but not the only possible scenario of the formation of planetesimals and evolution into comets. Advantages of the pebble-collapse model of comet formation are that it is entirely based on a large body of empirical evidence from the laboratory as well as on established numerical simulations of hydrodynamic processes.
The pebble-collapse model explains a number of observations of comets, particularly of comet 67P, such as the low bulk density, ultra-low mechanical strength, low thermal conductivity, presence of ultra-low-density dust particles and internal homogeneity. However, there are observations of comet 67P that the model cannot explain. Among these are the observed layering, the presence and variety of macroscopic ( 1 m) geologic features and the measured increase of permittivity in a shallow ( 25 m) subsurface region, for which alternative formation models may have explanations (unless they are evolutionary processes). These models are discussed in Weissman et al. 2020 [4].