The brittle boulders of dwarf planet Ceres

We mapped all boulders larger than 105 m on the surface of dwarf planet Ceres using images of the Dawn framing camera acquired in the Low Altitude Mapping Orbit (LAMO). We find that boulders on Ceres are more numerous towards high latitudes and have a maximum lifetime of $150 \pm 50$ Ma, based on crater counts. These characteristics are distinctly different from those of boulders on asteroid (4) Vesta, an earlier target of Dawn, which implies that Ceres boulders are mechanically weaker. Clues to their properties can be found in the composition of Ceres' complex crust, which is rich in phyllosilicates and salts. As water ice is though to be present only meters below the surface, we suggest that boulders also harbor ice. Furthermore, the boulder size-frequency distribution is best fit by a Weibull distribution rather than the customary power law, just like for Vesta boulders. This finding is robust in light of possible types of size measurement error.


INTRODUCTION
Boulders on planetary bodies bear information on past and present surface processes. In particular, the boulder properties and spatial distribution are related to the bulk properties of the parent body and the surface environmental conditions. On terrestrial planets, processes like impact cratering, volcanism, and mass wasting are typically responsible for the boulder formation. Degradation of boulders may result from processes like comminution by impacts and weathering, which on bodies with water and/or an atmosphere can include chemical weathering. Over the last decades, observations by spacecraft have revealed the existence of boulder populations on small airless Solar System bodies such as comets (Pajola et al. 2015(Pajola et al. , 2016, asteroids (Lee et al. 1996; Thomas et al. 2001;Michikami et al. 2008;Küppers et al. 2012;Jiang et al. 2015;Michikami et al. 2019;Dellagiustina et al. 2019), icy satellites (Pajola et al. 2021), and the protoplanet (4) Vesta (Schröder et al. 2020). In the absence of an atmosphere and volatiles like water, there are only a few processes that can produce and destroy boulders. The most important formation mechanisms are the destruction of a parent body (Michel et al. 2020) and spallation during large impacts (Krishna & Kumar 2016). The former is thought to be responsible for the boulderdominated surfaces of rubble-pile asteroids (Fujiwara et al. 2006;Michikami et al. 2019;Dellagiustina et al. 2019), whereas the latter process dominates on asteroids suspected to be more monolithic (Lee et al. 1996;Thomas et al. 2001;Küppers et al. 2012). Destruction by small impacts ) and thermal stress weathering (Delbo et al. 2014;Molaro et al. 2017) are the most important degradational processes.
Dwarf planet (1) Ceres maintains a position somewhat in between small bodies and the terrestrial planets, in the sense that it is a large, volatile-rich world, yet without atmosphere (Russell et al. 2016). As such, boulders on its surface may be affected by more processes than on small bodies, but by fewer than on the larger, more complex terrestrial planets. Here, we investigate the boulder population of Ceres and compare it with that of Vesta (Schröder et al. 2020). Both bodies were imaged by the same camera aboard the Dawn spacecraft (Russell & Raymond 2011). Such a comparison benefits from the fact that Vesta and Ceres have comparable distances to the Sun and very similar surface gravities (Basilevsky et al. 2013). So, any differences between the respective boulder populations may relate to compositional differences, with Ceres' crust harboring water ice, phyllosilicates, and salts De Sanctis et al. 2016;Prettyman et al. 2017), and Vesta's crust being basaltic (De Sanctis et al. 2012). The limited spatial resolution of the global Dawn image data set restricts our study to clasts larger than 100 m, for which Bruno & Ruban (2017) suggested the term megaclasts. But "boulder" has typically been used for clasts on small airless bodies irrespective of their size, and for consistency with the Vesta study we retain the term boulder.
Here, we study the global boulder population of Ceres using similar methods as those used for the Vesta boulder population (Schröder et al. 2020), as described in Sec. 2. The results of our analysis are reported in Sec. 3. We searched all Dawn images acquired in the Low Altitude Mapping Orbit (lamo) for boulders and determine general statistics of the global population related to boulder sizes and numbers (Sec. 3.1). We determine the size-frequency distribution (SFD) of boulder populations of individual craters and that of the global population (Sec. 3.2). The boulder SFD is traditionally fit with a power law, but the Vesta boulders rather follow a Weibull distribution. We evaluate whether the same holds true for Ceres boulders. We investigate the spatial distribution of boulders in and around individual craters, as well as the distribution of craters with boulders across the globe (Sec. 3.3). Furthermore, we estimate the average boulder lifetime by comparing the boulder density around craters for which an age estimate is available, and assess the Basilevsky et al. (2015) prediction that meter-sized boulders on Ceres have the same lifetime as on Vesta (Sec. 3.4). In Sec. 4 we discuss our results and the implications of the observed differences with the Vesta boulder population.

Boulder mapping
Boulders on Ceres can only be distinguished in framing camera images acquired in the Low Altitude Mapping Orbit (lamo) at an altitude of around 400 km and lower orbits of the extended mission (Russell et al. 2007). The framing camera is a narrow-angle camera with a field-of-view of 5.5 • × 5.5 • (Sierks et al. 2011). lamo coverage of the illuminated surface was near-complete for the camera's clear filter, but color imaging was sparse. The clear filter (F1) is a polychromatic filter with 98% transmission in the 450 to 920 nm wavelength range (Sierks et al. 2011). lamo images were acquired between 16 December 2015 and 27 August 2016 and have a typical scale of 35 m per pixel 1 (Roatsch et al. 2017). The average scale of 59 lamo images that we used in our analysis (at least one for each crater with boulders) is 35.8 ± 1.3 meter per pixel. The boulder-finding procedure was identical to that followed for Vesta boulders (Schröder et al. 2020). In summary, the second author browsed the entire data set of lamo clear filter images and identified, measured, and mapped all boulders using the J-Ceres GIS program, which is a version of jmars 2 (Christensen et al. 2009), after which the first author reviewed the results for accuracy and completeness. Boulders were identified as positive relief features in projected images at a zoom level of 1024 pixels per degree. The lamo resolution is about 230 pixels per degree at the equator, so this represents a zoom factor of about 4. Boulder size was determined using the J-Ceres crater measuring tool, which draws a circle around a boulder fitted to 3 points that are selected by the user on the visible boulder outline. The measurement uncertainty is about a single pixel. The limited accuracy of pointing information for lamo images leads to mismatches between projected images. We used small craters inside and outside the crater as tie points to align the projected images to a Ceres background mosaic and relative to each other. All this leads to uncertainty in the location of boulders on the order of 500 m. We are confident that we could reliably identify boulders with a size of at least 3 pixels (105 m), although a criterion of 4 pixels (140 m) is more likely to ensure that mapping is complete (Schröder et al. 2020;Pajola et al. 2021). We did not distinguish between boulders located either inside or outside the crater rim, a choice that we justify in Sec. 3.3.
The illumination conditions at the time of imaging affect the visibility of a boulder, mainly by the strength of its shadow. The photometric angles at the center of the lamo images, calculated for an ellipsoid Ceres, are plotted as a function of latitude L in Fig. 1. We see that the illumination conditions for imaging, and thereby boulder visibility, systematically changed according to latitude. The spacecraft looked at nadir most of the time (low emission angle), but the incidence and phase angle increased with L. If we distinguish three latitudinal zones as "low" (|L| < 30 • ), "mid" (30 • < |L| < 60 • ), and "high" (|L| > 60 • ), the average incidence angle at the image center is ι = 48 • ± 4 • for low latitudes, ι = 62 • ± 5 • for mid-latitudes, and ι = 75 • ± 5 • for high latitudes. For a spherical boulder that is half-buried in a plane surface, the maximum length of the shadow is l = r(cos −1 ι − cos ι), with boulder radius r and incidence angle ι. A boulder at low latitudes will cast the shortest shadow, with l = 0.83r. A boulder with a diameter of 3 and 4 pixels will cast a shadow of 1.2 and 1.6 pixels, respectively. Especially for the larger diameter, the shadow is long enough to be well visible. Thus, although boulders will be more easily recognized in high-latitude images, we are confident that all boulders with a diameter of 4 pixels can be recognized in low-latitude images. Nevertheless, it is likely that we missed boulders with a 3-pixel diameter in low-latitude images due to their short shadows. In Fig. 2 we investigate how the increase of ι with latitude affects our mapping. The figure shows lamo images of two craters with abundant boulders: high-latitude crater Jacheongbi (69 • S, ι = 78 • ) and low-latitude crater Unnamed17 (10 • S, ι = 42 • ). Both craters appear fresh, with boulders that are easily recognized by the shadows they cast on their ejecta blankets. At lamo resolution, the blankets appear equally smooth at either incidence angle. The rightmost panels show the distribution of boulders as mapped using J-Ceres. The stronger shadows of the Jacheongbi boulders did not lead us to recognize them in higher numbers compared to Unnamed17, which confirms that differences in visibility may only be consequential at the most extreme incidence angles (Wilcox et al. 2005). The figure also shows that half of Jacheongbi's interior is in shadow. Boulders are abundant on the sunlit half of the crater floor, which suggests that boulder numbers in high-latitude craters are severely underestimated, with potentially important consequences for the SFD.

Power law SFD
Various practices for displaying boulder SFDs seen in the literature were discussed by Schröder et al. (2020). In this paper, we display the SFD in both cumulative and differential format, using the incremental (binned, histogram) version for the latter (Colwell 1993). The cumulative distribution of boulders on Solar System bodies is often assumed to follow a power law, for which the number of boulders with a size larger than d is: with α < 0 the power law exponent and N tot the total number of boulders larger than d min . The exponent of a cumulative distribution of a quantity that follows a power law is identical to that of the associated incremental differential distribution with a constant bin size on a logarithmic scale, if the logarithmic bins are chosen wide enough (Hartmann 1969;Colwell 1993). The power law exponent is best estimated from the SFD by means of the maximum likelihood (ML) method (Newman 2005;Clauset et al. 2009). The ML power law exponent (α < 0) is estimated directly from the boulder size measurements aŝ with d i the size of boulder i and N the total number of boulders with a size larger than d min . The standard error ofα is plus higher order terms, which we ignore. The estimator in Eq. 2 is unbiased only for sufficiently large sample size. Clauset et al. (2009) also provide details to a statistical test that evaluates whether a power law is an appropriate model for the data 3 . The test randomly generates a large number of synthetic data sets according to the best-fit power law model (specified byα and d min ), and calculates for each the Kolmogorov-Smirnov statistic, which is a measure of how well the synthetic data agree with the model. A p-value, defined as the fraction of synthetic data sets that have a larger statistic than the real data set, quantifies how well the power law performs. The authors adopt p < 0.1 to mean rejection of the power law model. We fitted power laws to the global boulder population, but also to the populations associated with individual craters to investigate possible variations of the exponent over the surface. To evaluate whether any such variations found are meaningful, we simulate the Ceres boulder population on the basis of the power law, using the observed population sizes of individual craters as input. For all craters we adopt the same exponent, namely that of the power law that best fits the global population. The continuous power law probability distribution is also known as the Pareto distribution (Newman 2005). To simulate a size distribution of boulders associated with impact craters we draw a random variate U from a uniform distribution on (0, 1) using the randomu routine in IDL with an undefined seed. Then the boulder diameter follows a Pareto distribution, with α the power law index (associated with the cumulative distribution function, with α < 0). We adopted a minimum boulder diameter of d min = 140 m (4 pixels). We simulated populations for all craters and estimated the power law exponent for each using the ML method. We then compared the resulting distribution of exponents with the observed one. The power law exponent estimated according to Eq. 2 is biased for low boulder numbers (Clauset et al. 2009), which was illustrated for simulated boulder populations by Schröder et al. (2020). Another potential source of bias is measurement error of boulder size. While the boulders described in this paper are large in absolute terms, they typically measure only a few pixels across, and measurement errors on the order of a pixel can be expected. Here we investigate the consequences of measurement error using boulder populations simulated according to Eq. 4. First, we consider a group of "craters", each with a simulated boulder population of different size. The SFD of all boulder populations follows a power law with exponent −4, with a minimum boulder size of 40 m. We modified the boulder sizes according to three different definitions of measurement error, where we distinguish between systematic and random errors. Measurement errors are sized one pixel of 35 m, equal to the spatial resolution of the Ceres lamo images. We estimated the power law exponent for each crater, including only boulders larger than 4 pixels (140 m) in the fit. Note that we can only assess how many boulders have met this requirement after performing the simulation. Figure 3 shows the results. In (a), the boulder sizes are all measured correctly, and the power law exponent is retrieved reliably for craters with large boulder numbers. The negative bias at small boulder numbers inherent in Eq. 2 can be recognized clearly. In (b), boulder sizes are affected by measurement errors with a random character: Sizes are either decreased by 1 pixel (size overestimated), increased by 1 pixel (size underestimated), or left unchanged, each with equal probability of 1/3. Figures (c) and (d) explore the consequences of systematic measurement errors, by either underor overestimating all boulder sizes by 1 pixel. Systematic errors can be introduced by the method of measuring. Figure 3 shows that all types of measurement error lead to biased power law exponents. Underestimating the boulder sizes increases the exponent by a little less than unity (shallower power law), while overestimating the sizes decreases the exponent by unity (steeper power law). The bias is stronger for overestimating the sizes than for underestimating, which cause the random measurement errors in (b) decrease the exponent by a little less than unity (steeper power law). Another consequence of random measurement errors is that the exponent converges only at larger boulder numbers than without errors, which is not reflected in the formal uncertainty of the exponent (Eq. 3).
Another aspect of this problem is the shape of the cumulative SFD. Figure 4 investigates how the shape changes for the same three cases of random and systematic measurement errors. For each case we generated four populations of the same size, again with a power law exponent of −4. In (a), the boulder sizes are all measured correctly, and the SFD follows the straight line of the power law up to a diameter of about 200 m. Beyond this size, the simulated curves diverge considerably due to chance. In (b), the boulder sizes are affected by measurement errors with a random character, which steepens the SFD slightly. In (c), boulder sizes are systematically underestimated, which makes the SFD shallower and introduces a slightly convex curvature. In (d), boulder sizes are systematically overestimated, which steepens the SFD considerably and introduces a slightly concave curvature. Our simulations demonstrate that bias due to measurement error is unavoidable when typical boulder sizes are on the order of a few image pixels. The results in Figs. 3 and 4 may allow us to identify or predict such bias for the Ceres boulder population.

Weibull SFD
The SFD of the Vesta boulder population is better described by a Weibull distribution than a power law (Schröder et al. 2020). We investigate whether the same holds true for the Ceres boulder population. The Weibull distribution was initially derived empirically, and is often used to describe the particle distribution resulting from grinding experiments (Rosin & Rammler 1933). Where the power law follows naturally from a single-event fragmentation that leads to a branching tree of cracks that have a fractal character, the Weibull distribution results from sequential fragmentation (Brown & Wohletz 1995). Because we only include boulders larger than a certain size in the fit, we employ a left-truncated Weibull distribution with the cumulative form (Wingo 1989): where N is the number of boulders larger than d min . We estimate the Weibull parameters α and β = 3(γ + 1) from the boulder sizes d i > d min using the ML method.
To maximize the log-likelihood function, these two equations must be satisfied: We findβ from a simple grid search, andα by insertingβ.

General statistics
We identified a total of 4423 boulders on the surface of Ceres with a diameter larger than 3 image pixels (105 m), of which 1092 were larger than 4 pixels (140 m). All boulders are associated with impact craters. The details of all craters with at least one boulder larger than 4 pixels (n = 58) are listed in Table 1. First, we summarize some general statistics of the global boulder population. Figure 5a shows a (weak) correlation between the number of boulders per crater and crater size. The largest craters in the sample have a number of boulders that is much smaller than expected. For example, Occator is the largest crater with a diameter of 90 km, and it has 26 boulders larger than 4 pixels. From the general trend in the figure, we can expect to find a number at least an order of magnitude larger for a crater of its size. But many of its boulders may have been destroyed or hidden from view by the largescale flows that are present in and around the crater (Schenk et al. 2019). In fact, three out of the largest four craters in our sample show evidence of such flows: Azacca, Ikapati, and Occator (Buczkowski et al. 2018a;Krohn et al. 2018;Schenk et al. 2019). Boulders around these craters are also difficult to distinguish from features resulting from partly submerged topography. The fourth crater in this group, Gaue, is very old Pasckert et al. 2018). Age must affect the number of boulders, as they are destroyed over time. The figure also shows the number of boulders associated with craters on Vesta (Schröder et al. 2020), where we adopted the same minimum size (140 m) as for the Ceres boulders. Clearly, the number of boulders produced on average by impacts on Ceres is larger than on Vesta for craters of the same size.
We investigate the relation between the size of the largest boulder (L) and size of the crater (D) in Fig. 5b. Being a single measurement, the largest boulder size is a poor statistic, but is nevertheless often used to characterize boulder populations (Thomas et al. 2001;Küppers et al. 2012;Jiang et al. 2015;Schulzeck et al. 2018a). We compare the Ceres distribution with the relation provided by Lee et al. (1996) for craters formed in rocky targets (L = 0.25D 0.7 with L and D in m) and with the empirical range established by Moore (1971) for a selection of lunar and terrestrial craters (L = 0.01 1/3 KD 2/3 with K ranging from 0.5 to 1.5). The former relation represents more or less the upper limit of the latter range. The largest boulders on Ceres do not agree well with either relation from the literature, with sizes that are almost independent of the size of the craters. Again, it is mostly the largest craters that break the trend, probably because the aforementioned flow features and old age. The largest boulder we found on Ceres is a 500 m large block on the rim of Jacheongbi (figure inset), a relatively large crater (27 km). The figure also shows the largest boulders of Vesta craters. While the Vesta data also do not perfectly agree with either relation from the literature, most fall within the range of Moore (1971). On average, the largest boulders on Vesta are somewhat smaller than those on Ceres.

Size-frequency distribution
We aggregate all boulders counted on the surface of Ceres to find the cumulative power law exponent of the global boulder population. We note that the resulting global SFD is biased, as the boulder populations of the largest craters were almost certainly decimated by large scale flows (see Sec. 3.1). Figure 6 shows the SFD, both in cumulative and differential representation. At the top of the differential plot we show the uncertainty in size resulting from a 1 pixel measurement error. We chose a logarithmic bin size of 0.07 with the boulder size in meters, ensuring that the size is on the order of the measurement error at the larger end of the scale. As the error is larger than the bin size at the smaller end of the scale, we can expect boulders to end up in adjacent bins merely by chance. We recognize the characteristic roll-over of the distributions towards smaller diameters, caused by the limited spatial resolution and the measurement error. We fit two power laws to the data with the ML method, one with the minimum boulder size (d min ) fixed, and the other with d min estimated by the ML algorithm. When fixing d min to 4 pixels (140 m), we find a power law exponent of α = −5.8 ± 0.2 (n = 1092, black dashed lines in Fig. 6). By extrapolating the power law to smaller diameters we find that the number of boulders with a diameter around 3 pixels may be severely underestimated; the observed number in the bin closest to the 3 pixel limit is 2328, but the extrapolated, expected number is about 4000. The counts for boulders larger than 4 pixels are probably close to complete. We note that the counts at the largest diameters do not match well with the power law, both in the cumulative and differential representation. The statistical test provided by Clauset et al. (2009) confirms that this power law is not a good model for the data (p = 0). When we let the ML algorithm itself choose the minimum boulder size, we find a larger d min = 169 m and a steeper power law with α = −6.7 ± 0.3 (n = 400, red dashed lines in Fig. 6). The statistical test indicates that this power law is a good model for the data (p = 0.37).
We also estimated the power law exponent for individual craters that have at least 6 boulders larger than 4 pixels (Table 1), and plot these as a function of number of boulders in the population in Fig. 7. The figure also includes three simulations of the power law exponent distribution of individual craters. The simulation uses the observed population sizes and adopt the best-fit power law exponent for the global boulder population (α = −5.8). The observations and simulations agree in showing a large scatter in the exponents for smaller population sizes and their expected negative bias (Clauset et al. 2009;Schröder et al. 2020). The degree of scatter is similar for both simulated and observed data, indicating that the observed variety in power law exponents is merely due to differences in population size and not some physical property of boulders or craters. However, the observed exponents of craters with small boulder populations are typically more negative than in the simulations. Additionally, the power law exponent of the crater with the largest number of boulders, Jacheongbi, is further from −5.8 than that of the simulated "Jacheongbi's". It is almost as if the observed distribution of exponents is skewed with respect to the simulated distribution. This suggests that the power law model does not correctly describe the boulder SFD. Schulzeck et al. (2018a) independently counted boulders around several Ceres craters and fitted power laws to the SFD, using ML to estimate both the exponent and the minimum diameter. We compare their exponents with ours in Fig. 8. The exponents agree within the error bars for the craters with more than 70 boulders (Jacheongbi, Nunghui, and Unnamed11). The exponents for the other three craters (Juling, Ratumaibulu, and Unnamed17) agree less well, which is not surprising given the small size of their boulder populations. The exponents for the crater with the largest number of boulders (Jacheongbi with 160 boulders) match closely (−4.5 ± 0.4 versus −4.4 ± 0.7). This suggests that our counts are consistent with those of Schulzeck et al. (2018a). In Sec. 2.2 we uncovered evidence for bias resulting from measurement errors in the boulder sizes, which were probably similar in magnitude for Schulzeck et al. (2018a). Can this bias be responsible for the unusual steepness and the convex shape of the SFD of the global boulder population? Given that the boulder size measurements were almost certainly subject to errors of at least one image pixel, the "true" SFD is probably less steep (see Figs. 3 and 4). The power law exponent may be smaller by about unity, but with a value somewhere between −4.8 and −5.7, the SFD would still be unusually steep. Measurement errors also affect the shape of the SFD (Fig. 4). Random errors would actually lead to a slightly more concave shape of the SFD. Only systematically underestimating the boulder sizes would lead to a convex SFD, but this would also tend to make it shallower. Thus, we cannot attribute the downturn of the SFD towards large sizes to measurement error, and it must be an intrinsic property of the Ceres boulder population.
We conclude that the power law is not a good model for the SFD of boulders larger than 4 pixels. The ML algorithm was able to find an acceptable power law for a larger minimum size (d min = 169 m), but there is no reason to exclude the wellresolved boulders larger than 4 pixels but smaller than 169 m (more than half of the total) from the fit. This is the same situation as for the Vesta boulder SFD, for which the Weibull distribution proved to be a better model than the power law (Schröder et al. 2020). The best-fit Weibull distribution for the Ceres global boulder population has N = 1092, α = 1.32, and β = 0.45 (Fig. 9). The fractal dimension D f = 3 − β for the cracks in the rock is 2.5. The Weibull distribution fits the SFD better than the power law, and, contrary to the latter, it does not predict that the number of boulders with a size of 3 pixels is massively underestimated. Figure 9 also shows Weibull distributions for the Vesta global boulder population (Schröder et al. 2020). There are two best-fit curves, one including and the other excluding boulders of Marcia, the largest crater on Vesta. Just like the largest Ceres craters, Marcia shows evidence of flows, which may have destroyed many of its boulders. As Vesta is smaller than Ceres, its Weibull distributions plot below the Ceres distribution. Given the uncertainty surrounding the boulder populations of the largest craters on both worlds, a detailed comparison of the Weibull parameters provides little insight.

Spatial distribution
In Fig. 10, we plot the distribution of all boulders with a size of at least 3 pixels (105 m) on a color composite map of Ceres. In the same figure and on the same scale, we also show the distribution of boulders larger than 105 m on Vesta using counts from Schröder et al. (2020). With a surface area 3.2 times that of Vesta, the Ceres map is much larger. On both bodies, boulders are confined to craters. On Vesta, boulders are mostly absent from a large area of low albedo that appears to be enriched with carbonaceous chondrite material (Denevi et al. 2016;Schröder et al. 2020). The situation is different on Ceres. Several distinctly blue craters (Occator, Haulani, Kupalo) have boulders, which is consistent with blue being a marker of youth (Schmedemann et al. 2014). Large areas are devoid of boulders, especially at lower latitudes, but these do not stand out in terms of color or albedo like on Vesta. Even though the poles were partly in shadow during lamo, we find many boulders there.
We quantify the boulder abundance in three latitude zones (low, mid, and high) on both Ceres and Vesta in Fig. 11. We calculated both the density of boulders and the density of craters with at least one boulder larger than 105 m by dividing the total number of boulders/craters in a latitudinal zone by the total surface area of the zone (including shadowed terrain), calculated under the assumption that Ceres and Vesta are spheres with radii of 469 and 263 km, respectively. Adopting Poisson error bars allows to assess whether any differences are the result of chance. The boulder density graph (Fig. 11a) confirms our visual impression that the density is higher at the high latitudes of Ceres, despite the fact that polar terrain was partly in shadow. The small (Poisson) error bars indicate that this is very likely not due to chance. The Ceres boulder density at mid-latitudes is also a little higher than that at low latitudes. The boulder density at low latitudes is more uncertain than the error bar indicates. First, boulder numbers may be underestimated because of the limited visibility of the shadows cast by smaller boulders (see Sec. 2.1). Second, boulder counts for Occator and Haulani, with their large scale flows, are uncertain (see Sec. 3.1). In contrast to Ceres, the Vesta boulder density does not vary with latitude, within the Poisson error margins. It is a little lower than the Ceres boulder density at low and mid-latitudes and much lower at high latitudes. The density of craters-with-boulders (Fig. 11b) also increases on Ceres toward high latitudes, although the correlation is not as strong due to the larger error bars. The density of craters-with-boulders on Vesta does not vary with latitude, within the error margins. The crater density may be similar at high latitudes on both bodies, although the error bars are large. Interestingly, the density of craters-with-boulders on Vesta is significantly higher than that on Ceres for low and mid-latitudes, despite the boulder density being lower. This implies that, on average, craters on Ceres have more boulders than on Vesta, consistent with Fig. 5a.
Examples 4 of the spatial distribution of boulders around individual craters are shown in Fig. 12. Boulders are located both inside the crater and outside the rim, typically, all within one crater radius. The figure distinguishes between boulders in two different size classes, but sorting of boulders according to their size is not evident, consistent with the findings of Schulzeck et al. (2018a). There are no accumulations of boulders with a size range considered in our global study at the foot of steep slopes, which argues against a formation of such boulders by post-impact weathering. High-resolution images (scales of 3-10 m/px) acquired in the Dawn extended mission show evidence for boulder transport on crater walls, such as bounce marks on unconsolidated talus and boulders at the downslope end of tracks. Figure 13a shows an example of boulders collected at the foot of a crater wall. Such boulders are consistently smaller than those we consider here. As very high-resolution images are only available for a very small fraction of the surface, we do not consider them for our global study. We therefore decided to group boulders inside and outside the craters together and treat them as a single population. Another image from the extended mission, in Fig. 13b, shows clusters of boulders that appear to derive from former, larger boulders. Such fields of debris may originate from either impact of the larger boulder on the surface or weathering and/or erosion in place, demonstrating that boulders disintegrate by fracturing.

Boulder lifetime
Boulders degrade over time and eventually disappear from the surface. Basilevsky et al. (2015) predicted that the survival time of meter-sized boulders on Ceres is very similar to that on Vesta, based on estimates of the potential impactor flux and the expected impact velocities. Schröder et al. (2020) determined a survival time of about 300 Ma for Vesta boulders, much larger than the ∼ 10 Ma predicted by Basilevsky et al. (2015). The authors attributed this apparent discrepancy to the fact that the boulders in their sample were one to two orders of magnitude larger than the meter-scale associated with the prediction. The boulders in our sample are even larger than the Vesta boulders studied by Schröder et al. (2020) because of the lower lamo image resolution at Ceres, but typically by a factor of two rather than an order of magnitude. It should therefore be possible to test the Basilevsky et al. (2015) prediction of similar survival times on Vesta and Ceres, if accurate age estimates are available for our Ceres craters.
Estimating the age of a crater is typically done by counting smaller craters on a selected area on the ejecta blanket and modeling the resulting SFD. Just as for Vesta, two alternative chronologies have been used to model crater SFDs for Ceres: The lunar-derived model (LDM) adapts the lunar production and chronology functions to impact conditions on Ceres, whereas the asteroid-derived model (ADM) derives a production function by scaling the observed SFD from the main asteroid belt to the SFD of Ceres craters (Hiesinger et al. 2016). Most papers on the topic of Ceres dating employ both chronologies and provide two age estimates for a particular crater. Table 2 lists craters for which age estimates are available. The uncertainty associated with the age is typically large, as the two chronologies can yield widely different values.
Additional sources of uncertainty are the choice of counting area and the assumed strength of the target surface. The tabulated ages were estimated assuming an impact into hard rock. Williams et al. (2018) found that the ADM (and, presumably, the LDM) ages are much larger for a rubble surface (e.g. Cacaguat: 14.4 instead of 3.3 Ma, Rao: 133 instead of 30.4 Ma).
We define "areal density" as the total number of boulders identified in and around a crater divided by the crater equivalent area, calculated as the area of a circle with the diameter for that crater (Table 1). We determined the areal density of boulders larger than 105 m in and around the craters in Table 2. Some areal densities are unreliable. Boulder numbers are underestimated for craters at high latitudes, which were partly in the shadow during lamo (Shennong, Unnamed4/26/28/36/44). Boulder numbers are uncertain for craters that show post-impact modifications in the form of flows (Haulani, Ikapati, Occator). The distribution of boulders around craters with a reliable density was shown in Fig. 12. We relate the boulder density to crater age in Fig. 14a, where the age ranges are spanned by the LDM and ADM estimates. The two variables are anti-correlated. There is a large degree of scatter in the data, which may be due to the fact that the boulder density also depends on latitude (see Sec. 3.3). The data suggest that the maximum boulder survival time is around 150 Ma, where we note that the age of the oldest crater in the figure (Gaue) is very uncertain. Support for this maximum age comes from craters of this age for which we did not find boulders: One such, unnamed, crater at (162 • E, +78 • ) has an estimated age of 89-252 Ma . Two others, Messor at (234 • E, +50 • ) and an unnamed crater at (186 • E, +23 • ), have estimated ages of 96-192 Ma and 88-205 Ma, respectively (Scully et al. 2018). All craters estimated to be younger than 150 Ma in the papers referenced in Table 2 have boulders. Given the uncertainty in the crater age estimates due to the different models, we adopt an uncertainty of 50 Ma for the boulder survival time. Whereas the correlation with crater age is expected, we also find that the boulder density is correlated with crater size (Fig. 14b). This may at least partly be explained by the fact that large craters are, on average, older than small craters. Another explanation might be that larger craters sample different, deeper crustal layers in a stratified crust, a concept that we discuss in the next section. Basilevsky et al. (2015) estimated the survival time of meter-sized boulders on Vesta, Ceres, and the Moon by predicting the impactor velocity and density distributions. Boulders on Vesta and Ceres should live equally long, but 30 times shorter than on the Moon. Schröder et al. (2020) found that the maximum boulder survival time on Vesta is 300 Ma, the same as that of lunar boulders (Basilevsky et al. 2013). The authors attributed this apparent contradiction to the fact that the Vesta boulders in their study are larger than 60 m, rather than meter-sized. In other words, large boulders live longer than small ones. At 150 ± 50 Ma, the survival time of the Ceres boulders in our sample is only half that of Vesta boulders, despite being larger (> 105 m).
Thus, the lifetime of Ceres boulders may be less than half that of Vesta boulders of the same size.

DISCUSSION
In the previous section we described the properties of the population of large boulders on Ceres and compared them to those of the Vesta population. Our major findings are: (1) Ceres craters have, on average, more boulders than Vesta craters of the same size, (2) the largest boulders are, on average, somewhat larger for Ceres craters than Vesta craters of the same size, (3) the SFD of the global boulder population is better described by a Weibull distribution than a power law for both Ceres and Vesta, (4) boulders on Ceres are more numerous at high latitudes than at midand low latitudes, in contrast to Vesta, and (5) boulders live shorter on Ceres than on Vesta. How can we reconcile these findings? Let us start with finding (3), which supports the idea that the SFD of particles on the surface of small bodies follows a Weibull distribution rather than a power law (Schröder et al. 2020). Unfortunately, uncertainties regarding the boulder populations of the largest craters on both Ceres and Vesta prevent us from a meaningful comparison of the Weibull parameters.
To address the other findings, we need to consider how boulders degrade. The dominant mechanism responsible for degradation are weathering due to thermal stress (Delbo et al. 2014;Molaro et al. 2017;El Mir et al. 2019) and meteorite impacts . The efficacy of weathering due to diurnal thermal cycling or thermal shock correlates with the rate of surface temperature change (dT /dt). Rapid temperature changes occur at sunrise and sunset and during daytime shadowing. On quickly rotating bodies such as Vesta and Ceres, sunrise and sunset are the main drivers of dT /dt (Molaro & Byrne 2012). Vesta and Ceres have rotation rates of 5.34 h and 9.07 h, so dT /dt at the terminator should be larger on Vesta. Moreover, the thermal cycling rate is higher on Vesta, and a boulder of the same age will have experienced more thermal cycles than on Ceres. Therefore, boulders of identical lithology would degrade faster on Vesta through thermal stress weathering. This is inconsistent with the shorter boulder lifetime on Ceres (finding 5), but consistent with the fact that boulders for a given crater diameter are larger on Ceres (finding 2). It is a different story for boulder degradation due to meteorite impacts. Because of their similar locations in the main asteroid belt, Ceres and Vesta experience similar impact regimes, in terms of impactor size distribution, flux, and velocity (Basilevsky et al. 2015). Therefore, boulders of identical lithology would degrade equally fast on both worlds through meteorite impacts.
The crusts of Vesta and Ceres are of different composition, leading to differences in the compressive and tensile strengths that control the resistance to stress. Vesta's crust is mostly an assemblage of eucritic basalts and pyroxene cumulates ( (Prettyman et al. 2017;Schmidt et al. 2017). On average, the materials in Ceres' crust are mechanically weaker than those in Vesta's crust, so its boulders are less resistant to stress. This applies to both thermal stress and stress caused by meteorite impacts. Given the considerable uncertainties in quantifying the effects of thermal stresses and determining thermal strain thresholds (e.g., Boelhouwers & Jonsson 2013), it is impossible to predict which of these competing effects dominates (higher thermal stress experienced on Vesta or lower resistance to stress on Ceres). Nevertheless, our finding (5) that boulders on Ceres have a shorter lifetime suggests that the lower mechanical strength of the Ceres crust is primarily responsible. The lower crustal strength would also lead to the formation of a larger crater on Ceres than Vesta for an identical impactor. So a crater of the same diameter is, on average, younger on Ceres than on Vesta. Boulders disappear over time, and younger craters have, on average, more boulders than older craters of the same size. Therefore, craters of the same size are expected to have more boulders on Ceres than Vesta, explaining finding (1).
The prevalence of large boulders at high latitudes may be explained by a higher rate of physical weathering and boulder breakdown at lower latitudes as compared to higher latitudes. Ceres has a obliquity of only 4 • at present 5 (Ermakov et al. 2017), and therefore the diurnal temperature waves are expected to be larger at lower latitudes (Hayne & Aharonson 2015). The duration of sunrise and sunset would also be shorter at lower latitudes, increasing dT /dt. Both effects would lead to relatively faster boulder breakdown by thermal stresses at equatorial latitudes, consistent with finding (4). As water ice is likely abundant in the subsurface, Ceres boulders may harbor a significant fraction of ice. Ice would be relatively stable inside these large boulders, just as it is stable just meters below the surface (Fanale & Salvail 1989), yet ice-rich boulders could be more prone to degradation by thermal stress, as fractures may be widened by sublimation, further weakening the boulder structure. The hypothesis that Ceres' boulders are rich in water ice is consistent with Rivkin et al. (2014), who considered the question why Ceres does not have a dynamical family. Large impacts on Ceres would produce escaping fragment asteroids (essentially liberated boulders), but the absence of a dynamical family led the authors to suggest that such escaped fragments are ice-rich and prematurely destroyed by sublimation.
To identify other factors that may contribute to latitudinal differences in the areal densities of boulders and craters with boulders, we consider the independent global data set of floor-fractured craters. Craters with fractured floors are indicators for several possible processes, including updoming due to cryomagmatic activity, as for instance beneath Occator crater (Buczkowski et al. 2018a). Cryomagmatism may indicate the presence of a crustal column beneath floor-fractured craters that is enriched in volatiles and therefore mechanically weaker. A volatile-rich and weak target material is expected to eject boulders with these properties, which would be more susceptible to degradation. Twenty-one floor-fractured craters have been identified on Ceres (Buczkowski et al. 2018b), seven of which are marked in Table 2. Floorfractured craters tend to have low boulder densities, consistent with the expectation that boulders ejected from floor-fractured craters degrade faster, resulting, on average, in a lower boulder density. However, the floor-fracturing is likely related to post-impact processes, and it is not clear whether the target substrate was already weaker before the impact. Support for the hypothesis of a crust with a mechanical strength that depends on latitude comes from the observation that five of the seven floor-fractured craters in Table 2 (Azacca, Haulani, Ikapati, Occator, and Tupo) display concentric fracturing beyond the crater rim, suggestive of creep of a low-viscosity, and therefore mechanically weak, subsurface layer (Otto et al. 2019). Such craters with concentric fractures are mostly located between latitudes 46 • S and 34 • N. This is consistent with landslide morphology that suggests the presence of a relatively weak layer at low-to mid-latitudes (Chilton et al. 2019). This layer thins towards the poles and overlies a stronger layer, in agreement with lower temperatures at high latitudes increasing crustal viscosity and strength (Bland et al. 2016). Boulders excavated from the weaker layer would tend to degrade faster, and more boulders would be retained in the polar regions. Moreover, the subsurface ice content is lower in low-to mid-latitudes regions (Schmidt et al. 2017), hence boulders there would harbor more impurities like phyllosilicates, which tend to lower the albedo, raise temperatures, and enhance sublimation (Rivkin et al. 2014). This effect would lead to faster degradation of lowand mid-latitude, less ice-rich boulders as compared with polar boulders with possibly a higher ice content. Therefore, there may be a causal relationship between the boulder density and the pre-impact properties of the crust. So, while boulders may harbor water ice, the complexity of Ceres' crust, with its laterally and vertically varying properties (Bland et al. 2016;Otto et al. 2019;Park et al. 2020), precludes any definite conclusion on the observed distribution of boulder and boulder crater densities across the surface.

DATA AVAILABILITY
Dawn framing camera images are available from NASA's Planetary Data System at https://pds.nasa.gov/. Our Ceres boulder data, including maps of the boulder distribution around all craters, are available for download at https://doi.org/10.5281/ zenodo.4715154.

ACKNOWLEDGMENTS
We are grateful for technical support provided by J-Ceres developer Dale Noss and his team at ASU. We thank Maurizio Pajola and an anonymous reviewer for helpful suggestions to improve the manuscript. Table 1. All craters on Ceres with at least one boulder larger than 4 pixels (d > 140 m). Crater and boulder diameters are D and d, respectively, and α is the power law exponent of the (cumulative) boulder SFD as derived with the ML method (only for craters with n d>4px > 5).     . Investigating the effect of measurement error on the shape of the cumulative distribution when the true power law exponent is −4. The pixel size is 35 m, and the adopted power law is only shown for sizes "measured" larger than 4 pixels (140 m, dashed line). The four curves of different color represent repeated simulations of the same population. a. Boulder sizes measured correctly. b. Sizes either correctly measured, overestimated by 1 pixel, or underestimated by 1 pixel with equal probability. c. Sizes underestimated by 1 pixel. d. Sizes overestimated by 1 pixel. The population sizes for the four different scenarios were adjusted to achieve a cumulative number of boulders of around 1000 at 140 m. Figure 5. General statistics of boulders associated with craters on Ceres (this paper) and Vesta (Schröder et al. 2020). a. Number of boulders larger than 140 m as a function of crater diameter. b. Diameter of the largest boulder as a function of crater diameter. The uncertainty of the Ceres boulder size derives from a measurement error of 1 pixel (Vesta boulder error bars are omitted for clarity). The empirical range given by Moore (1971) for selected lunar and terrestrial craters is shown in gray. The dashed line is the relation given by Lee et al. (1996). The inset shows the largest boulder identified on Ceres, a 500 m sized block on the rim of Jacheongbi crater.  . Power law exponents for all craters with a population of at least 6 boulders larger than 4 pixels (n = 39). The observed exponents were derived by fitting a power law to the data of each crater. The best fit power law index for the observed global boulder distribution is α = −5.8 ± 0.2 (dashed line with gray confidence interval). The crater with the largest number of boulders (160) is Jacheongbi. We compare the observations to three simulations. The simulated exponents were derived by fitting randomly generated boulder distributions, assuming a Pareto distribution with α = −5.8 (dashed line), using the number of boulders in the population of each crater as input. Figure 8. Comparison of the power law exponents for the craters Juling, Jacheongbi, Nunghui, Ratumaibulu, Unnamed11, and Unnamed17 as determined in this paper and by Schulzeck et al. (2018a). Filled symbols refer to exponents that are more reliable, being associated with populations of more than 70 boulders.  , and the Vesta map has filters centered at 650 nm, 555 nm, and 438 nm in the RGB channels (Schröder et al. 2013). The center of the maps is at (lat, lon) = (0 • , 0 • ). Figure 11. The distribution of boulders larger than 105 m on Vesta and Ceres along latitude L. We aggregate the boulders in three latitude ranges: "low" (|L| < 30 • ), "mid" (30 • < |L| < 60 • ), and "high" (|L| > 60 • ). a. Number density of boulders. b. Number density of craters with at least one boulder larger than 105 m. The error bars derive from Poisson statistics. Figure 12. Spatial boulder distribution for several craters for which an age estimate is available. Green, small dots represent boulders with a size between 3 and 4 pixels (105 m < d < 140 m). Red, large dots represent boulders larger than 4 pixels (d > 140 m). The FC2 image number is indicated in the top right.   (Table 2). a. Density versus crater age. b. Density versus crater diameter. The error bars on the density were calculated assuming the number of boulders follows a Poisson distribution. The open symbols represent craters whose boulder density is unreliable, either because of uncertain boulder identifications or because the associated craters were partly in the shadow (underestimate).