Erosion of volatiles by micro-meteoroid bombardment on Ceres, and comparison to the Moon and Mercury

(1) Ceres, the largest reservoir of water in the main-belt, was recently visited by the Dawn spacecraft that revealed several areas bearing H$_2$O-ice features. Independent telescopic observations showed a water exosphere of currently unknown origin. We explore the effects of meteoroid impacts on Ceres considering the topography obtained from the Dawn mission using a widely-used micro-meteoroid model and ray-tracing techniques. Meteoroid populations with $0.01-2$ mm diameters are considered. We analyze the short-term effects Ceres experiences during its current orbit as well as long-term effects over the entire precession cycle. We find the entire surface is subject to meteoroid bombardment leaving no areas in permanent shadow with respect to meteoroid influx. The equatorial parts of Ceres produce $80\%$ more ejecta than the polar regions due to the large impact velocity of long-period comets. Mass flux, energy flux, and ejecta production vary seasonally by a factor of $3-7$ due to the inclined eccentric orbit. Compared to Mercury and the Moon, Ceres experiences significantly smaller effects of micro-meteoroid bombardment, with a total mass flux of $4.5\pm1.2\times10^{-17}$ kg m$^{-2}$ s$^{-1}$. On average Mercury is subjected to a $50\times$ larger mass flux and generates $700\times$ more ejecta than Ceres, while the lunar mass flux is $10\times$ larger, and the ejecta generation is $30\times$ larger than on Ceres. For these reasons, we find that meteoroid impacts are an unlikely candidate for the production of a water exosphere or significant excavation of surface features. The surface turnover rate from the micro-meteoroid populations considered is estimated to be $1.25$ Myr on Ceres.


Apr 2021
Ceres is the only dwarf planet with a global high-resolution shape model thanks to the Dawn mission and the Dawn spacecraft's Framing Camera 1 . The dwarf planet itself is surrounded by the zodiacal cloud, a cloud of dust and meteoroids enveloping the entire solar system, that is mainly sourced from main-belt asteroids and short/long-period comets (Nesvorný et al. 2010). Since the zodiacal cloud is increasingly denser with decreasing heliocentric distance (Leinert et al. 1981), Ceres is expected to experience a smaller meteoroid flux and impact velocities than those seen by the Moon and Mercury.
Two sets of telescopic observations suggest that Ceres has at least temporarily a water exosphere (A'Hearn & Feldman 1992; Küppers et al. 2014). Despite negative results from follow-up observations with different observational facilities (Rousselot et al. 2011;Roth et al. 2016;McKay et al. 2017;Roth 2018;Rousselot et al. 2019), the production of the water exosphere on Ceres was the subject of many works (Tu et al. 2014;Schorghofer et al. 2016;Formisano et al. 2016;Landis et al. 2017;Villarreal et al. 2017;Schorghofer et al. 2017). Landis et al. (2019) has analyzed the surface evolution based on the flux of impactors larger than 100m in diameter, however the effects of smaller meteoroids are yet to be assessed. Therefore, it is natural to investigate this gap and address how different meteoroid populations imprint their activity on the surface of Ceres, how they are affected by the complex topography, and how they influence the regions with exposed water-ice.
Until recently, the availability of precise shape models and detailed dynamical meteoroid models was scarce; to say nothing of the existence of models to study the effects of meteoroid bombardment using ray-tracing techniques. Pokorný et al. (2020) analyzed the effects of meteoroid bombardment on the lunar surface, specifically on both lunar poles using detailed topography derived from Lunar Orbiter Laser Altimeter (LOLA) observations (Smith et al. 2010). Advanced solar illumination models (Mazarico et al. 2018) show that the lunar poles harbor many cold traps and permanently shadowed regions (PSRs), consistent with spacecraft data (e.g., Diviner instrument Paige et al. 2010). On the other hand, Pokorný et al. (2020) showed that meteoroids, due to their broad range of impact directions, reach even the deepest craters on both lunar poles, including those that are permanently shadowed from the Sun. Pokorný et al. (2020) also showed that the crater walls, that were more inclined toward the high energetic meteoroid flux concentrated close to the ecliptic, experienced a higher rate of impact gardening compared to the crater floors. The floors were partially shielded from the meteoroid flux and experienced the high energetic impacts at grazing angles.
Compared to the Moon, Ceres is on an inclined, eccentric orbit, embedded inside the source region of a significant portion of the zodiacal dust cloud (main-belt asteroids; e.g. Nesvorný et al. 2006). The obliquity of Mercury and the Moon is stable in the long term (Goldreich 1966), which allows the existence of long-lived permanently shadowed environments, whereas the obliquity of Ceres oscillates between = 2 • and = 20 • with a period of 24.5 kyr and complicates the existence of permanently shadowed regions (Ermakov et al. 2017). However, as shown in Ermakov et al. (2017), PSRs indeed exist on Ceres and can potentially retain water-ice. In addition to the ice deposits in the PSRs, Combe et al. (2019) showed nine areas with H 2 O absorption features detected in Dawn VIR spectra (2.0µm line). These areas are not currently accessible by direct solar radiation, however, as shown in Pokorný et al. (2020) they might be accessible to meteoroid impacts. These meteoroid impacts bring exogenous energy to PSRs that can remove volatiles from PSRs. It has been suggested that destabilized volatiles might produce a tenuous exosphere at airless bodies such as Mercury and the Moon (Cintala 1992), which motivates us to quantify the effect of meteoroids on these interesting regions at Ceres. Landis et al. (2019) discussed various explanations for the temporary existence of the water exosphere on Ceres: (a) sublimation from sub-surface water-ice tables (Fanale & Salvail 1989;Prettyman et al. 2017); (b) sublimation from transient surface exposures of water-ice ); (c) sputtering by solar energetic particle events from surface ice (Villarreal et al. 2017); (d) seasonal, optically thin water ice polar deposit (Schorghofer et al. 2017). The telescopic evidence from A' Hearn & Feldman (1992); Küppers et al. (2014) suggests 3 -6 kg s −1 of water vapor produced to sustain the observed exosphere. None of these mechanisms or their combination provides such a high rate (Landis et al. 2019). The residence time of the water exosphere was shown to be around 7 hr (Schorghofer et al. 2016), thus a continuous active source is needed to sustain the longer lasting exosphere observed by Küppers et al. (2014). Since meteoroids provide a quasi-continuous source of energy to surfaces of airless bodies, we aim to estimate the meteoroid impact driven production rates and provide a missing piece to this puzzle. Moreover, meteoroids and dust may also erode water ice deposits.

Ceres topography model
For the representation of Ceres, we use an object (OBJ) file containing 1,579,014 vectors and 3,145,728 faces (triangles) derived from the digital terrain model (DTM) constructed from Framing Camera 2 (FC2) images taken during Dawn High Altitude Mapping Orbit (HAMO). HAMO DTM covers approximately 98% of the cererean surface with a lateral spacing of ≈ 136.7 m/pixel 2 . We converted the original ICQ file to OBJ file using an example code described here https://sbnarchive.psi. edu/pds3/dawn/fc/DWNCSPC 4 01/DOCUMENT/ICQMODEL.ASC. We also provide the working version of this code at the project's GitHub page (see Software section).

Ray-tracing code
We improved our own procedure from Pokorný et al. (2020) in terms of computational speed and combined it with two external libraries written in C++: (1) tinyobjloader for loading up to 10 million polygon models and (2) fastbvh -a Bounding Volume Hierarchy algorithm that uses axisaligned bounding box (AABB) trees to efficiently calculate line-triangle intersections. Our program loads the OBJ file using tinyobjloader into memory and then constructs using fastbvh, an AABB tree that allows us to quickly evaluate ray-face intersections for the entire triangular mesh. For each combination of longitude, latitude, and velocity in the meteoroid model our code determines whether each surface triangle is reachable or obstructed/shadowed by any other triangle. We also calculate the incidence angle ϕ for each surface triangle and meteoroid direction/velocity combination independently. Each surface triangle represents a plane for which we calculate a normal vector pointing outside of the object, − − → n sur . Then cos ϕ = − − − → n sur · − − → e imp , where − − → e imp is the velocity vector of the impacting meteoroid normalized to unity. For impacts perpendicular to the surface cos ϕ = 1, whereas for meteoroids at grazing angles cos ϕ → 0, since ϕ → 90 • .
2.3. Model for meteoroid environment at Ceres and its variations over current Ceres' orbit 2 The mission data product can be found at https://sbnarchive.psi.edu/pds3/dawn/fc/DWNCSPC 4 01/DATA/ICQ/ CERES SPC181019 0512.ICQ The meteoroid model in this work uses the same constraints and configuration as the model used in Pokorný et al. (2018Pokorný et al. ( , 2019Pokorný et al. ( , 2020 for studies of Mercury's and the Moon's meteoroid environments. This means we do not change the configuration of the zodiacal cloud/meteoroid model with respect to those previous studies, ensuring they are comparable with each other. Our model combines contributions of four meteoroid populations originating from main-belt asteroids (MBAs), Jupiter-Family comets (JFCs), Halley-type comets (HTCs), and Oort Cloud comets (OCCs). These four populations dominate the meteoroid mass and number density flux in the inner solar system (Nesvorný et al. 2010(Nesvorný et al. , 2011a, while the outer solar system is mostly dominated by Edgeworth-Kuiper Belt (EKB) meteoroids (Poppe et al. 2019). The details of the model used here are summarized in Table 1 and references therein. We calculate the distribution of directions and velocities of impacting meteoroids for one orbit from January 1 st , 2015 to August 18 th , 2019 in 10-day intervals resulting in 169 individual snapshots of the meteoroid environment at Ceres. Each of these snapshots provides the meteoroid mass flux distributed in sun-centered longitudes λ − λ and latitudes β with 2 degree resolution (i.e., directions) and impact velocities with 2 km s −1 resolution, thus providing the full three-dimensional map of velocity vectors for meteoroids in each meteoroid population separately. An example of the meteoroid environment map and the impact velocity distribution is shown in Figure 1. The directionality of meteoroids on Ceres is similar to that seen on other airless bodies like Mercury (Pokorný et al. 2018) or the Moon (Pokorný et al. 2019). The mass flux is dominated by meteoroids impacting Ceres from directions close to the orbital plane (β ≈ 0 • ) and within 90 • of the ram/apex direction (−180 . HTC and OCC meteoroids preferentially originate from the apex direction (a circle around [λ − λ , β] = [−90 • , 0 • ] with 45 • radius). MBA meteoroids have the smallest relative impact velocities V imp < 10 km s −1 and impact Ceres preferentially from higher Sun-centered latitudes |β| > 60 • and/or from the sun/anti-sun directions (λ − λ ≈ −180 • and λ − λ ≈ 0 • ). The impact velocity distribution ( Figure 1B) shows that MBA meteoroids have a median V imp 50% = 5.5 km s −1 , JFC meteoroids are slightly faster with V imp 50% = 9.2 km s −1 , whereas the long-period comet sources provide much more energetic impactors with V imp 50% = 20.5 km s −1 for HTC meteoroids and V imp 50% = 25.3 km s −1 for OCCs meteoroids.
Due to the non-zero eccentricity and inclination of Ceres, the meteoroid environment undergoes significant changes during one orbit (Fig. 2). During one of its orbital cycles, Ceres crosses the ecliptic twice (gray regions in Fig. 2); this is where we record the global and local maxima of the meteoroid mass flux M for MBA and JFC meteoroids. This is not unexpected because MBA and JFC meteoroid models show that they are concentrated close to the ecliptic (Nesvorný et al. 2010(Nesvorný et al. , 2011a. On the other hand, HTC and OCC meteoroids are unaffected by Ceres' departure from the ecliptic plane due to the broad range of inclinations of their parent bodies (Nesvorný et al. 2011b;Pokorný et al. 2014). Overall, JFC meteoroid mass flux dominates the total mass influx at Ceres throughout the entire orbit due to the dominance of JFC meteoroid in the Zodiacal cloud. This stems from the population mixing ratios obtained from Carrillo-Sánchez et al. (2016) and Pokorný et al. (2019), where JFCs dominate the terrestrial flux by a factor of 10 compared to other meteoroid populations. The departure from the ecliptic is not the only factor shaping the meteoroid mass flux at Ceres. From spacecraft observations and modelin, we know that the Zodiacal Cloud density inside 1 au increases with heliocentric distance, r (e.g, Leinert et al. 1981;Pokorný et al. 2019), while the outer portions of Table 1. Description of meteoroid dynamical models used in this work. The meteoroid model used here has six free parameters: the collisional lifetime multiplier F coll (Pokorný et al. 2014), differential size-frequency index α, and the average daily mass influx at Earth M for each of the four populations in metric tons per day (1000 kg per day or 11.57 g s −1 ) (Carrillo-Sánchez et al. 2016;Pokorný et al. 2019). For more detailed information refer to references in the table or Pokorný et al. (2018). the Zodical Cloud show constant density (Poppe et al. 2019). With Ceres being inside the main-belt, we can expect significant differences from what is observed inside 1 au for MBA and JFC meteoroids. For the outer solar system sources (HTC and OCC), the meteoroid mass flux scaling should be similar to that observed/modeled inside 1 au. We calculate the proportional changes of M for all four sources assuming the single power-law scaling M ∝ r α and positions close to the ecliptic |z| < 0.01 au (labels above arrows in Fig. 2). While not directly comparable to the Zodiacal Cloud density, due to the velocity changes that Ceres experiences during its orbit, the proportionality of HTC (∝ r −2.65 ) and OCC (∝ r −2.23 ) meteoroids is quite similar to those modeled in Pokorný et al. (2019). On the other hand, JFC (∝ r −3.08 ) and most significantly MBA (∝ r −4.71 ) meteoroids show much steeper scaling with heliocentric distance than what was observed inside 1 au. The two main factors that drive the proportionality of different populations are: 1) the location of the source region with respect to the target (Ceres) where MBA and JFC meteoroids are ejected on semimajor axes close to that of Ceres while long-period comet meteoroids can be considered distance source regions. This causes larger variations for the close source regions due to the fact that some meteoroids are released at semimajor axes smaller than that of Ceres; and 2) the impact velocity of different populations at Ceres which is modulated by the populations' semimajor axis/eccentricity/inclination distributions. MBA and JFC meteoroids can acquire very small impact velocities at Ceres due to their orbital similarity, while for the long-period meteoroids this is insignificant. We calculate that, averaged over one of its current orbits around the Sun, the mean meteoroid flux at Ceres is M = 4.49±1.18×10 −17 kg m −2 s −1 , which is approximately 9.5 times smaller than that at

Quantities induced by the meteoroid bombardment
In this paper we discuss four different quantities that result from the bombardment of the Cererean surface by interplanetary meteoroids: (1) the meteoroid mass flux M, (2) the meteoroid energy flux E, (3) the ejecta mass produced by impacting meteoroids P + , and (4) the area of craters produced by meteoroid bombardment A.
In Section 2.2 we describe how we determine whether each surface element is reachable by the flux of meteoroids incoming for a selected direction (i.e., it is not shadowed by any feature on the surface of the dwarf planet). This analysis results in a list of cos ϕ values (cosines of incident angles) for each triangular element on Ceres, and the value S(λ−λ , β) which represents the percentage of shadowing of a particular element (S = 0 is a completely shadowed element, S = 1 is unobstructed). Here, λ − λ and β are the sun-centered longitude and latitude, i.e., the angles representing the meteoroid directionality. The meteoroid model we use here gives us the five-dimensional information for each time snapshot we analyze in this manuscript: for each combination of directions λ − λ , impact velocity v imp , and meteoroid diameter D we obtain the meteoroid number flux N met . Assuming that all meteoroids are spheres, we get the meteoroid mass flux M met and cross-sectional area of impacting meteoroids per unit time A met as: where we adopt the meteoroid bulk density ρ met = 2, 000 kg m −3 used in meteoroid models we show in Table 1. Equipped with these quantities, we then calculate for each triangular element on Ceres M, E, P + , and A, defined as follows: where C = 7.358 km −2 s 2 is a scaling constant determined from the laboratory experiments reported by Koschny & Grün (2001).We assume the impactor-to-crater area ratio F cr = 63, which is based on fitting results from Table 1 in Koschny & Grün (2001). The scaling constant C describes the amount of ejecta the surface produces. The value of C used here characterizes impacts of glass projectiles into ice-silicate surfaces. The ejecta mass production rate, P + , is likely orders of magnitude smaller for regolith dominated surface than those covered by water-ice surfaces, as suggested by the discrepancy of the meteoroid modeling and Lunar Dust Experiment discussed in Pokorný et al. (2019). See Section 6.2 for a discussion on the scaling of impact processes for different surface types with respect to other methods and lab experiments. While the area of craters produced by meteoroid bombardment A is a useful quantity, we also define the e-folding lifetime T A that defines the time that craters produced by meteoroids will take to cover 1 − e −N of the surface, where N is the number of lifetimes the surface experienced. Since A represents the area of craters per unit surface area pre time, it has units of s −1 , and then T A = A −1 . This means that if T A is 10 Myr, then in 10 Myr approximately 63% of the surface will be covered by craters, while in 30 Myr 95% of the surface will be covered by craters, and so on.
In order to decrease the computation time we investigate only the centroids of each triangular surface element, thus S(λ, β) is either one or zero. To evaluate this simplification we run our raytracing procedure using 10, 50, and 100 points per surface triangle. When we average these results over all directions of the meteoroid environment, we saw only negligible changes in all the quantities analyzed here. This is in contrast to solar radiation simulations, whose source is a small disc (< 0.1 • in radius as viewed from Ceres) in the sky and requires finer resolution (Mazarico et al. 2018).

Ceres' rotation, pole precession, obliquity oscillation
The effect of Ceres' rotation period of P rot = 9.074 hr is implemented by rotating the meteoroid map (i.e., the celestial sphere) along Ceres' rotational axis and averaging the contribution of different position into one mean map. Since our meteoroid model time stamps are recorded in 10 day intervals, our simplified rotation method introduces ∼ 2% uncertainty/error to the daily flux, which is much smaller than the intrinsic uncertainties in the meteoroid model; for the model uncertainty discussion see Pokorný et al. (2018Pokorný et al. ( , 2019 and Section 6.1. The precession of the pole is implemented by rotating the Ceres rotational axis around the axis perpendicular to the orbital plane and adopting a pole precession rate of 210 kyr (Ermakov et al. 2017). The default rotational axis pointing in our model is set to right ascension α = 291.42751 • and declination δ = 66.76043 • , which translates to ecliptic longitude λ = 11.18622 • and latitude β = 81.55038 • , or longitude and latitude with respect to the orbital plane λ O = 328.24905 • and β O = 85.96854 • . The rotation from ecliptic coordinates to orbital coordinates is where R a (y) is the rotation matrix with respect to axis a which represents a clock-wise rotation by an angle y. In a more specific way: The effects of obliquity/axial tilt are calculated for each obliquity value separately and then the contribution of each obliquity regime is weighted by the time Ceres spends in that obliquity regime, based on Fig. 1 in Ermakov et al. (2017).

RESULTS -PRECESSION CYCLE AVERAGE
First, we look at results of the meteoroid bombardment of Ceres, averaged over the entire precession cycle, i.e., approximately 210 kyr (Ermakov et al. 2017). Here, we assume the proper orbital elements of Ceres 3 , with a = 2.7671 au, e = 0.116198, and sin(i) = 0.167585 (i = 9.64744 • ). Figure 3A shows in detail the entire surface of Ceres color coded by the value of the meteoroid mass flux M. Despite the seemingly significant difference between the equatorial and polar regions, the average M is almost constant over the entire surface, where the difference between the maximum and minimum values is 7.8%. The average value over the entire surface is M = 4.54±0.04×10 −17 kg m −2 s −1 , where M peaks at mid-latitudes, around 60 • . The local topography of Ceres plays a negligible role when the meteoroid environment is averaged over a longer timescale. This is also apparent on both Cererean poles (right side in Fig. 3A), where we see < 2% differences between maximum and minimum values.
We also highlight nine areas of interest from Combe et al. (2019) that showed detections of exposed H 2 O from Visible and InfraRed (VIR, the mapping spectrometer of the Dawn mission) remote sensing observations (white labeled rectangles in Fig. 3). The coordinates of these nine areas are shown in Table 2. Unlike solar irradiation which is almost a point source at the distance of Ceres, the meteoroid environment is able to influence any surface element on the Cererean surface and even the most significant features do not provide permanent shadowing from the meteoroid bombardment. Ermakov et al. (2017) identified seven bright crater floor deposits (BCFDs) that are correlated with the most persistend PSRs. Since all of these seven areas have high latitudes (|β| > 69.7 • ), they are all comparable to area (C) -Messor crater from Combe et al. (2019).
The global map of the ejecta mass production rate P + is shown in Figure 3B. For P + the equatorial areas experience the maximum exposure, while the polar regions experience 43% less exposure (i.e., a factor 1.76 difference). A similar situation holds for the meteoroid energy flux, where the pole receives ∼ 19% less E than the equatorial regions. The reason for the difference between the global behavior of M and P + shown Fig. 3A is due to the different impact velocities of meteoroid populations. As shown in Fig. 2 the HTC and OCC meteoroids have the highest velocities with median velocities 2.1-2.6 times larger than the JFC population that dominates the meteoroid mass flux. High velocity meteoroids from both HTC and OCC populations are impacting Ceres preferentially from the invariable plane which is close to the ecliptic. Due to their high impact velocities, these meteoroids are not sensitive to Ceres' own orbital velocity variations. This results in a continuous high energy bombardment of the equatorial regions, while polar regions experience high energy impacts on grazing angles. The mean values for the meteoroid energy flux and ejecta production rate are E = 5.16 ± 0.31 × 10 −16 MJ s −1 m −2 and P + = 1.56 ± 0.22 × 10 −16 kg s −1 m −2 , respectively.
The last quantity we discuss here is the surface e-folding lifetime T A (Fig. 3C). The shortest efolding times are in the equatorial area, while T A increases toward the both poles. However, the difference between the maximum and minimum values of A is small: < 6%. Unlike P + and E the e-folding time is not modulated by the impactor velocity, since we assume a constant impactor-to-  Costello et al. 2018Costello et al. , 2020. The typical depth of meteoroid-induced craters is a factor 2-4 larger than the particle radius (Koschny & Grün 2001). When we average over the entire size range simulated here, we find a mean crater depth around 100 − 300 µm. The mean crater depth depends on the location of Ceres as well as the meteoroid population causing the impacts. This is due to the different size-frequency distributions that meteoroid populations simulated here have when impacting Ceres. Note, that the Koschny & Grün (2001) explored only impact velocities up to 10 km s −1 , so the crater depth has high uncertainty that we cannot quantify at the moment. A closer look at the nine areas showing exposed H 2 O signatures, the crater morphology and the general topography have a rather small effect even for the Juling crater that shows the largest ratio between the maximum and minimum values for all three quantities: R M = M max /M min = 1.03, R E = 1.24, R + P = 1.75 and R T A = 1.06. The only exception is the ejecta mass production rate P + , which for the Juling crater shows 75% difference between the crater rim P + = 1856 × 10 −16 kg m −2 s −1 and partially shadowed crater floor P + = 1063 × 10 −16 . This is due to its cubic dependence on the cosine of the incident angle and the v 2.46 imp velocity scaling, which emphasize the effect of fast meteoroids close to the ecliptic. The surface features facing the ecliptic produce significantly more ejecta than those that are effectively shadowed by the crater rim. This effect is most efficient for mid-latitudes 30 • < |β| < 50 • , where it provides ∼ 70% difference between the exposed and shadowed portion of the crater. This efficiency of shadowing drops toward the equator and both poles, where for the Messor crater we see a drop to 49% difference between the exposed/shadowed areas.

RESULTS -ONE ORBIT OF CERES
The image of the meteoroid environment imprint on Cererean surface averaged over one precession cycle of Ceres greatly simplifies the picture. The averaged picture is important for understanding the long-term surface evolution, but Ceres experiences significant meteoroid flux changes during the  Cererean year (Fig. 2). In this Section we analyze variations of several quantities over one orbit of Ceres, 4.61 years. Such a time segment is longer than the duration of the scientific stay of the Dawn spacecraft at Ceres (RC3 orbit phase started on April 23, 2015 and on October 31, 2018 the spacecraft ran out of its propellant).
In Fig. 2 we show that the meteoroid mass flux M impacting Ceres is most significantly modulated by the asteroid's distance from the ecliptic. To put this in a different perspective, we show the orbit of Ceres in Cartesian coordinates in Fig. 4. This Figure shows 12 different time records of M on Cererean surface averaged over one rotation period (9.1 hours) in SI units from January 1st, 2015 to March 21st, 2019. In Fig. 4 we see the effects of the orbital motion of Ceres on the global shape of the meteoroid mass flux. On January 1st, 2015 Ceres just passed the descending node and its negative zaxis velocity is close to its maximum, v z = −3.32 km s −1 (note that the median impact velocity of JFC meteoroids is V 50% (JFC) = 9.32±0.46 km s −1 ). This increases the number of impacts to the southern hemisphere, since Ceres is plunging downwards and the relative velocity of meteoroids impacting the southern hemisphere is increased, while the northern hemisphere experiences attenuated meteoroid mass flux. The opposite effect can be seen on April 20th, 2017, when Ceres is close to the ascending node and its positive z-axis velocity is close to its maximum v z = 3.40 km s −1 . When the z-axis velocity is close to zero, i.e., the sum of the argument of perihelion and true anomaly is sin(ω +f ) = 0, the meteoroid mass flux is symmetric around the Cererean equator as seen on July 14th, 2016 and June 14th, 2018.
The variations of M, P + , and T A for the Combe et al. (2019) areas (Table 2) are shown in Fig.  5. The mass flux during one orbit of Ceres shows a double peaked structure, where depending on the latitude of the feature, the global maximum occurs at true anomaly angle TAA = 106 • (for southern hemisphere features D -Juling and I -Baltay) or at TAA = 286 • (for nothern hemisphere features) coinciding with the ecliptic crossing (Fig. 5A). These mass flux spikes are correlated with the density of the zodiacal cloud that is densest close to the ecliptic and gets more tenuous further from the ecliptic. The north/south hemisphere seasonality is caused by the orbital velocity in the z-axis described in Fig. 4. As shown in Fig. 2, the TAA = 286 • peak is stronger than the TAA = 106 • peak due to Ceres' smaller heliocentric distance; the zodiacal cloud gets denser with decreasing heliocentric distance. The north/south asymmetry is changing as Ceres undergoes nodal precession and when averaged over the entire precession cycle the difference between the areas in the south and north are negligible (see Table 2).
The ejecta mass production rate P + true anomaly profile differs from that of M mainly due to two factors. First, P + is strongly dependent on the impact velocity of individual meteoroids (v 2.46 imp , see Eq. 5). The highest impact velocities come from OCC and HTC meteoroids (see median velocities for each population in Sec. 2.3), which shifts the dominance over P + from JFCs and MBAs to long-period comet particles. As we show in Fig. 2, the long-period comets are not very sensitive to the distance from the ecliptic which attenuates the variations of P + during its orbit to a factor of 2-3 (compared to a factor of 4-7 variations of M). Second, the 5-95 percentile (gray shaded area in Fig 5) is wider than for M due to the cos 3 of incidence angle dependence. This accentuates the effect of surface features, where surface patches experiencing almost perpendicular impacts produce significantly more ejecta than those subjected to grazing impacts.
The surface e-folding timescales T A show a similar trend in the north/south asymmetry as the two previous values, but their maxima are inverse with respect to the true anomaly angle (Fig. 5C). This is understandable, since the higher impactor flux produces shorter e-folding timescales. T A is not sensitive to the impact velocity, thus is similar to the mass flux variations with slightly smaller minimum/maximum ratios of 3-4.
Analysis of the Combe et al. (2019) areas showed significant variations of all quantities analyzed here during one Ceres' orbit. From Fig. 5 we can infer that areas on similar latitudes are undergoing very similar variations within one orbit. Since the Combe et al. (2019) areas are preferentially closer to the Cererean poles, many potentially interesting areas are left out. Furthermore, due to the short rotation period of Ceres, that longitudinally averages meteoroid bombardment effects at a given latitude, we can simply divide Ceres into latitudinal strips and analyze the entire dwarf planet as a whole. In Figure 6 we show variations of M, P + , and T A with true anomaly angle for 20 • -wide latitudinal strips, revealing the variations experienced by the entire surface. Lines in Fig. 6 can be used to determine M, P + , and T A for any point on the surface. The uncertainty of this approximation is quite small, because for each latitudinal strip the difference between the median value (color coded lines in Fig. 6) and the 5% or 95% percentile is < 35%. Figure 6 shows that the equatorial latitudes experience smaller variations in all quantities analyzed here, where the ratio between the maximum and minimum value during one orbit increases with latitudes closer to the Cererean poles. For example, the M max-min ratio for −10 • < β < 10 • is 2.43, while the same quantity for the north pole areas 70 • < β < 90 • is 6.88. As mentioned before,  Table 2 for description). The two peaks in M correspond to the ecliptic crossings and the maximum flux from MBA and JFC meteoroid populations. Middle: The same as the top panel but now for the ejecta mass production rate P + . Bottom: The same as the top panel but now for the surface exposure e-folding time T A . the north pole peak values are smaller than those on the south pole because Ceres passes through the ascending node closer to the Sun, hence, through the denser portion of the Zodiacal cloud.
Telescopic observations of the water exosphere on Ceres by A'Hearn & Feldman (1992); Küppers et al. (2014) suggest production rates of 3 − 6 kg s −1 . Let us assume an extreme case where all ejecta produced by meteoroids convert to the detectable water exosphere. We can obtain the water exosphere production rate via multiplying P + by the surface area covered with surface water-ice. Using our values of P + = 0.5 − 2.5 × 10 −13 kg m −2 s −1 , this would require a surface water-ice area of 1.2 − 12 × 10 13 m 2 to obtain production rates 3 − 6 kg s −1 , which is equivalent to a sphere with a radius 1000 − 3000 km. This would mean that even if the entire surface of Ceres was covered by water-ice, the meteoroid impacts would not be able to produce enough ejecta required to sustain the water exosphere observed by A'Hearn & Feldman (1992); Küppers et al. (2014). Suppose that the entire Occator crater with diameter of 92 km is covered in water-ice and continuously bombarded by meteoroids. The Occator crater area is about 6.65 × 10 9 m 2 and the ejecta mass produced per second is then 0.3 − 1.7 × 10 −3 kg s −1 . Such a hypothetical value is much smaller than the sublimation from known water-ice patches (0.16 kg s −1 Landis et al. 2019).
In the next Section we will compare these values with the identical quantities on Mercury and Moon and draw conclusions for the meteoroid bombardment effects on these three bodies.

COMPARISON OF CERES TO MERCURY AND MOON
Similarly to Sec. 3, we calculate the effects of meteoroid bombardment using our meteoroid model on surfaces of Mercury and the Moon (Pokorný et al. 2018(Pokorný et al. , 2019(Pokorný et al. , 2020. For Mercury we use the following orbital elements and escape velocity: a = 0.3871 au, e = 0.2056, i = 7.0056 • , and V esc = 4.25 km s −1 , whereas for the Moon we use a = 1.0000 au, e = 0.0167, i = 0.0005 • , and V esc = 2.38 km s −1 . The remaining orbital elements are assumed to be randomly distributed between 0 and 360 • . In order to simplify our calculation, we use the shape model of Ceres as a substitute for the shape models of Mercury and Mars. Since in this Section we aim to analyze the global models of three airless bodies, we assume that simplification is acceptable. To further check our assumptions, we compare our results to those calculated on a smooth sphere representing Mercury and the Moon. The values shown in this Section agree within 1%. Note that all three objects rotate much faster than their nodal precession periods and they are not in 1:1 spin-orbital resonance with the Sun and we can approximate that the meteoroid bombardment is longitudinally uniform. This ensures that the surface of each of these three bodies is uniformly affected with respect to the surface longitude. The 3:2 spin-orbit resonance of Mercury with respect to the Sun can possibly introduce some secondary non-uniform variations in effects of meteoroid bombardment. Mercury is still rotating with respect to the meteoroid environment, thus we expect that the effects of the meteoroid bombardment average out. We reserve the quantification of the 3:2 spin-orbit resonance for future work due to the necessity of simulating the entire Hermean two-year cycle in fine detail. For each of the three airless bodies, we analyze 60 segments that are uniformly distributed in surface latitude β (i.e, 3 • wide segments). The results for Ceres are in Fig. 7, where we show variations of M, E, P + , and A with respect to surface latitude. To make our results more readable and transferable, we fit simple four-parameter functions to each quantity. For M, E, P + we achieve a best fit using the Gaussian function G(β): where a 0 is the maximum amplitude of the function occurring at β = µ 0 , b 0 is the offset of the function in the y direction, µ 0 is the offset from the center in the x-axis (i.e., in the latitude), and σ 0 is the standard deviation. From Figure 7, it is evident that M is not a Gaussian distribution, but rather resembles the sum of two symmetric Gaussians. For fitting such a profile, we use a four-parameter function G sym (β): where the only difference between the two Gaussians is the sign at µ 0 . Furthermore, A on Mercury responds the best to fitting a cosine function C(β): where a 0 is the amplitude, b 0 is the offset of the function in the y direction, µ 0 is the offset from the center in the x-axis (i.e., in the latitude), and σ 0 describes the period of the cosine. These three functions are chosen for their simplicity and ease of interpretation. We tested more than 80 other probability density distributions available in the SciPy framework, but we find no distribution that would provide significantly better fits than our three functions. Table 3 shows the results of our function fitting to the latitudinal distributions of M, E, P + , and A for Mercury, Moon, and Ceres. For each planetary body, we show the four fit parameters, the difference between the maximum and the minimum value, and the ratios, R, between the three objects for the average values of quantities M, E, P + , and A. Mercury and Moon both have µ 0 = 0.0 except for cases when we fit G sym . This is because we enforce the latitudinal symmetry in our fits since both bodies have negligible flatness (i.e., almost perfectly spherical shape).
Mercury experiences the highest values of all quantities analyzed here, which are an order of magnitude higher than those at the Moon, and 2-3 orders of magnitude higher than those at Ceres. The meteoroid energy flux E and the ejecta mass production rate P + at Mercury are amplified with respect to those at Ceres due to a combination of higher impact velocities and meteoroid fluxes closer to the Sun. Mercury, assuming that both Mercury and Ceres have the same surface material, produces ∼ 700 times more ejecta than Ceres via meteoroid impacts. Similarly, the E is ∼ 400 times higher on Mercury with respect to Ceres. These quantities have been shown to have strong correlation with the existence of a tenuous dust cloud around the Moon Pokorný et al. 2019;, sustaining an exosphere of several metals around Mercury (Killen & Hahn 2015;Merkel et al. 2017;Pokorný et al. 2018), and potential stability of water-ice in permanently shadowed regions at Mercury and the Moon (Pokorný et al. 2020;Hayne et al. 2015;Deutsch et al. 2019).
The lunar surface water-ice stability was shown to suffer similar erosion rates from meteoroid bombardment (about 13 × 10 −8 m per year; Pokorný et al. 2020) as those from H Ly-α radiation (about 7 × 10 −11 m per year; Morgan & Shemansky 1991). This means that the surface water-ice on Ceres should be primarily excavated by extrasolar radiation because the meteoroid energy flux is ∼ 25 times smaller on Ceres as compared to that on the Moon. Meanwhile, the meteoroid bombardment should dominate the shadowed surface water-ice excavation on Mercury, due to an order of magnitude higher value of E as compared to effects of H Ly-α. The quantification of the surface water-ice loss through meteoroid impacts at Ceres is a complex process that involves the combination of impact vaporization (e.g., Eq. 10 in Cintala 1992) and loss of ejecta produced in impacts. Since both of these effects depend linearly on the meteoroid mass flux, scale with some power of impact velocity, and depend on material characteristics, it is beyond the scope of this article to obtain more concrete values. We can only conclude that the loss of surface water-ice through meteoroid impacts at Ceres is an order of magnitude smaller compared to that on the Moon, and 2-3 orders of magnitude smaller with respect to that at Mercury.
The area of craters produced by meteoroid bombardment A shows that fresh deposits on the Hermean surface are covered by meteoroid-induced craters 74× faster than those on Ceres, whereas the lunar surface shows 12× faster rates. This means that while at Ceres the mean surface e-folding time is T A = 1.25 Myr, the mean value for the Moon is 107 kyr and for Mercury 17.1 kyr. These values ignore the effect of secondary impactors, which might enhance the gardening rates on the surface (for the effect on the Moon see Costello et al. 2018). As such, 95% of the surface of Mercury's surface Table 3. Comparison of four different quantities: rows (1,5,9) meteoroid mass flux M, rows (2,6,10) meteoroid energy flux E, rows (3,7,11) ejecta mass production P + , and rows (4,8,12) meteoroid-induced crater cross-section A at Mercury, Moon, and Ceres. Column (1) represents the selected quantity Q, columns (2-5) show the function fit parameters, column (6) shows the percent difference between the maximum and minimum values, and columns (7-9) show the ratios R for all quantities in this table with respect to other airless bodies analyzed here. For instance, R Moon =5.45 in the first row indicates M Mercu /M Moon , the ratio between the average meteoroid mass flux on Mercury with respect to that on the Moon. The last column (10) shows the function F used to fit the quantity. Quantities M, E, and P + are fitted using Gaussian distributions described in Eq. 8, whereas A is fitted using the cosine function described in Eq. 10. Parameters a 0 and b 0 are in SI units, while σ 0 and µ 0 are in degrees. is covered by meteoroid-induced craters within three e-folding lifetimes, i.e., in 50 kyr. However, the crater depth resulting from meteoroid bombardment in the size range in our model is < 1 mm according to laboratory experiments in Koschny & Grün (2001). One caveat of such an experiment is the absence of V imp > 10 km s −1 impacts that are dominating the impacts on lunar and hermean surfaces, thus the penetration depths could significantly change with new laboratory experiments.

Meteoroid model
There are several uncertainties that stem from the meteoroid model. Firstly, the meteoroid mass flux at Earth has an intrinsic uncertainty of about 50-60% based on the latest estimates (Carrillo-Sánchez et al. 2016. This uncertainty linearly scales all quantities in this paper for all three airless bodies discussed here. Secondly, the collisional lifetimes used in our meteoroid model are also subject to uncertainty as discussed for example in Pokorný et al. (2018) and Pokorný et al. (2019) for Mercury and Moon, respectively. We analyze different values of the collisional lifetime multiplier F coll ∈ [10, 50] to test the model sensitivity and compare it to the settings used in Pokorný et al. (2018Pokorný et al. ( , 2019. Due to Ceres' proximity to MBA and JFC source populations, the variations of the collisional lifetime has a negligible effect on the model results. The effect on the long-period comet meteoroids (HTCs and OCCs) is similar for Mercury, Moon, and Ceres and does not significantly alter the results. Thirdly, the size-frequency distributions (SFDs) of our model meteoroid populations are poorly constrained due to the fact that JFC meteoroids dominate the inner solar system budget and the direct measurements of different components of the meteoroid complex are extremely rare (see e.g, Section 2.4 in Pokorný et al. 2019). We ran our meteoroid model calculation using a range of mass indices [3.4, 4.6] for each meteoroid population separately, to test their sensitivity to different SFDs. Since our model is constrained by the meteoroid mass flux at Earth, the changes in the size-frequency distribution produce < 10% changes in the meteoroid mass flux M, meteoroid energy flux E, and ejecta production rate P + at Mercury and the Moon, similar to Pokorný et al. (2018Pokorný et al. ( , 2019. However, the SFD has a higher impact on Ceres, where we record up to 46% higher mass flux for a shallower SFD α = −3.4 as compared to our nominal model α = −4.0. On the other hand, steeper SFDs result in smaller mass flux by up to 25% for α = −4.6. The highest sensitivity to SFD stems from the main-belt meteoroids, due to their extreme proximity to Ceres and their high intrinsic collision probability with Ceres without the need for any dynamical evolution. Unlike M, E, and P + , the area produced by meteoroid bombardment A and consequently the efolding lifetime T A do not scale linearly with the meteoroid mass; they scale as M 2/3 . The meteoroid model we use here is scaled such that the mass flux at Earth is held constant (see Table 1). For this reason, A and T A are more sensitive to the the SFD setting than our other variables , which we see for all three airless bodies studied here. The values of A are approximately two-times larger for α = −4.6 compared to our nominal SFD α = −4.0, while for the shallow SFD α = −3.4 the values are two-times smaller compared to α = −4.0. This means, that for steeper SFDs the smaller particles dominate the total impactor cross-section area, while this value is attenuated for shallower SFDs. This impacts also the e-folding lifetime T A , which is ∼ 2× shorter for steeper SFD α = −4.6 and ∼ 2× longer for shallower SFD α = −3.4 as compared to values for α = −4.0. The intermediate values of SFD indices fall between the two extremes presented here.
Our meteoroid model represents a range of meteoroids with diameters between 10 µm and 2, 000 µm. Particles smaller than D = 10 µm exist in the solar system and add mostly to the number flux experienced by various bodies in the solar system. Their mass flux is small compared to our modeled sample due to their shallower SFD driven by the Poynting-Robertson drag (for SFD at Earth see Love & Brownlee 1993). Meteoroids with D ≤ 1 µm are effectively blown out of the solar system via radiation pressure (Burns et al. 1979) and do not significantly contribute to the quantities analyzed here. By expanding our model to smaller sizes, the results in this article would not change significantly, because our meteoroid model is scaled to provide a certain mass flux at Earth (Carrillo-Sánchez et al. 2016). The influence of meteoroids larger than those in our population D > 2, 000 µm is more difficult to estimate. We expect meteoroids of D = 2, 000 µm to dynamically resemble larger meteoroids because the Poynting-Robertson drag magnitude decreases with increasing particle diameter. This consequently increases the dynamical timescales of larger meteoroids making their dynamical evolution similar to their parent bodies. Asteroidal impacts at Ceres are a focus of Marchi et al. (2016). We are not aware of any study that would deal with impacts of comets at Ceres.

Scaling of impact processes
In Section 2.4, we establish that the crater-to-projectile cross-section ratio is F cr = 63 using the Koschny & Grün (2001) experimental results. This experiment showed that the crater diameters were not correlated with impactor velocity, thus the velocity part is missing in Eq. 6. On the other hand, a literature overview of impact processes at larger sizes and impact velocities < 8 km s −1 shows there is a strong correlation between the impactor velocity and crater size (Holsapple 1993).
In order to quantify differences between our estimates for the crater area production rate and scaling laws shown for larger and slower impactors than those we analyze in this article, we use the cratering volume V cr estimates from Table 1 in Holsapple (1993). We use the strength regime for our impact estimates because our impactors are smaller than 1 mm in radius. Assuming that all meteoroid-induced craters are spherical caps, we get for the volume of the crater where a is the crater radius, h is the crater depth. Assuming a crater radius-to-depth ratio of δ = a/h = 2.55, V cr , we get a simple form that allows us to estimate the crater radius a: Holsapple (1993) provides a general relation for the crater volume and the impactor characteristics: where C 1 is a material constant, m met is the impactor/meteoroid mass, d met is the impactor/meteoroid diameter, ρ met is the impactor/meteoroid bulk density, and is the velocity power index. Finally, the crater cross-section A is: For simplicity we assume ρ met = 2000 kg m −3 , the value we use in all meteoroid models in this article. From Eq. 14 we see that F cr , the crater-to-impactor cross-section ratio, is a function of impact velocity. We recalculate the values of the meteoroid-induced crater cross-section A for Mercury, the Moon, and Ceres that we show in Section 5 for two different surface compositions: dry soil and soft rock( Table 1 in Holsapple 1993), and summarize our results in Table 4. The influence of the impact velocity is the most important for Mercury, where the meteoroids that follow the Holsapple (1993) formulas are expected to produce on average 23.49 times larger area of craters on soft rock and 10.35 times larger area on dry soil. Then for the soft rock surface, we would expect that the mean surface e-folding time on Mercury is T A = 730 years. Note that formulas from Holsapple (1993) are not supported by experiments for the size and velocity regimes that we analyze in this article. All three airless bodies that we analyze here are commonly bombarded by meteoroids with much higher impact velocities and smaller sizes than those used in laboratory experiments. Our original decision to use Koschny & Grün (2001) is based on the fact that this experiment is the closest in impactor size and speed to the meteoroids that we model. Kato et al. (1995) present the results of an experiment with 15 × 10 mm cylinder impactors made of water-ice, polycarbonate, aluminium, and basalt. They test two types of targets: ice block and Table 4. Estimates of the average enhancement factor of meteoroidinduced crater cross-sections A using the Holsapple (1993) formulas (column 3) with respect to our default value F cr = 63 for two different surface compositions. Introduction of the impact velocity factor increases the average crater area produced by meteoroids by a factor of 23.49 assuming soft rock surface on Mercury as well as by a factor of > 4.74 for any combination of surface type and airless body. This would lead to a significant decrease of e-folding lifetimes T A for all three airless bodies that we analyze here. ground snow powder, with a maximum impact velocity of 1 km s −1 . This experiment uses an order of magnitude larger impactors than the largest meteoroid we model, and velocities much lower than those we record in our model. The Shrine et al. (2002) experiments showed that impact cratering of polycrystaline ice by 1-mm aluminium spheres depends heavily on impactor velocity/kinetic energy. Their maximum impact velocity was 7.34 km s −1 . The entire sample consisted of 16 shots with velocities between 1.07 and 7.34 km s −1 . Sommer et al. (2013) showed a summary of previous laboratory experiments and added impacts of iron meteorite and steel projectiles with velocities between 2.5 and 5.3 km s −1 onto dry and wet sandstone. None of the lab experiments except for Koschny & Grün (2001) recodred impact velocities larger than 7 km s −1 and particle diameters below 800 µm. These three works show that the impact cratering regime that the meteoroids we model in this article experience, is very poorly constrained by laboratory experiments.
A possible solution to the lack of experimental records is extensive numerical modeling using stateof-the-art Hydrocodes (Elbeshausen et al. 2009;Kraus et al. 2011;Stickle et al. 2020). However, in this case, we are not aware of any work that analyzes impacts and cratering of micron-sized meteoroids onto regolith or water-ice.

Axial tilt of Ceres
The effects of the axial tilt are negligible when assumed over time frames longer than one orbit. This is because the meteoroids impacting the surface of Ceres come from a broad range of ecliptic longitudes and latitudes, as opposed to the Sun, which is represented by a singular point on the celestial sphere. We test the effects of the axial tilt for the range of [0 • , 20 • ], similar to the range of axial tilts shown in Ermakov et al. (2017), and record < 5% differences between the extreme values of all quantities analyzed in this work.

CONCLUSIONS
We present the first model for micro-meteoroid bombardment effects on the dwarf planet Ceres. Using a detailed shape model, efficient ray-tracing code, and a widely-accepted meteoroid population model for the diameter range of 0.01-2 mm, we estimate the effects of the meteoroid bombardment on the entire Cererean surface analyzed over one precession cycle (Section 3), and over one current orbital period (Section 4). Finally, the effects meteoroid bombardment experienced by Ceres are compared to those on Mercury and the Moon (Section 5).
Here, we summarize the most important findings: • There are no permanently shadowed regions with respect to meteoroid bombardment. The local topography creates up to 80% difference between occluded and exposed regions (floor vs. rim) • The equatorial regions are producing on average 80% more ejecta than the polar regions, whereas the mass flux is more concentrated at the Cererean poles. However, the mass flux is almost uniform over the entire surface with only 2% variations. The surface e-folding lifetimes are ∼ 8% shorter at the Cererean poles as compared to equator.
• All areas showing detections of exposed H 2 O from Combe et al. (2019) and all areas of interest from Ermakov et al. (2017) experience similar rates of meteoroid bombardment (Table 2). Meteoroids smaller than 2 mm generate orders of magnitude less ejecta from patches of exposed ice than would be required to sustain the water exosphere observed by Küppers et al. (2014) and are much lower than estimated surface water-ice sublimation rates. The surface turnover rate is expected to be 1.25 Myr, much longer than the period of obliquity cycles.
• Ceres currently experiences a factor of 3-7 seasonal variations in mass flux, energy flux, and ejecta production along its current orbit. This is due to its inclined orbit that takes the dwarf planet far from the invariable plane, i.e., away from the densest parts of the Zodiacal cloud.
• Ceres experiences a ∼ 10× and ∼ 50× smaller mass flux than the Moon and Mercury, respectively. The meteoroid energy flux and ejecta production rate differences are significantly enhanced by higher impact velocities on the Moon and Mercury resulting in ∼ 30× and ∼ 600× larger effects on the Moon and Mercury, respectively (Table 3). The same area on Mercury is covered by primarily meteoroid-induced craters 74 times faster than the equivalent area on Ceres, while the lunar surface is covered on timescales 12 times shorter than at Ceres.