Macroporosity and Grain Density of Rubble Pile Asteroid (162173) Ryugu

Rubble pile asteroids such as (162173) Ryugu have large bulk porosities, which are believed to result from void spaces in between the constituent boulders (macroporosity) as well as void spaces within the boulders themselves (microporosity). In general, both macroporosity and microporosity are estimated based on comparisons between the asteroid bulk density and both the bulk and grain density of meteorite analogs, and relatively large macroporosities are usually obtained. Here we use semiempirical models for the macroporosity of multicomponent mixtures to determine Ryugu's macroporosity based on the observed size‐frequency distribution (SFD) of boulders on the surface. We find that Ryugu's macroporosity can be significantly smaller than usually assumed, as the observed SFD allows for an efficient packing of boulders, resulting in a macroporosity of 16% ± 3%. Therefore, we confirm that Ryugu's high bulk porosity is a direct consequence of a very large boulder microporosity. Furthermore, using estimates of boulder microporosity of around 50% as derived from in situ measurements, the average grain density in boulders is 2,848 ± 152 kg m–3, similar to values obtained for CM and the Tagish lake meteorites. Ryugu's bulk porosity corresponding to the above values is 58%. Thus, the macroporosity of rubble pile asteroids may have been systematically overestimated in the past.

The bulk porosity inside rubble pile asteroids can be separated into two contributions: the first one stems from the intrinsic porosity of rocks and boulders and is termed microporosity, while the second contribution refers to voids in-between particles and is termed macroporosity (Britt et al., 2002). The latter is directly related to the geometrical arrangement of the constituent blocks, also known as the packing state, which qualitatively describes the arrangement of particles and can vary between random loose and random close packings. Macroporosity of average C-complex asteroids was estimated to be 25%-30% (Britt et al., 2002), which is generally consistent with numerical models of the reassembly of blocks after a catastrophic disruption, which result in macroporosities of 20%-40% (Wilson et al., 1999). However, simulations suffer from unrealistically large lower cutoff sizes for the considered boulder population, such that rubble pile asteroids may still exhibit lower macroporosities.
Here we investigate the macroporosity of asteroid Ryugu using semiempirical models for the porosity of multicomponent mixtures of nonspherical, cohesive particles (Zou et al., 2011). Such models predict the macroporosity of granular material given the particle size as well as the particle shape distributions applying linear mixing and using the concept of controlling mixtures (A. B. Yu & Standish, 1991) to calculate the packing state. In general, polydisperse particle mixtures can have a macroporosity which is considerably smaller than the canonical ∼36% for a random close packing or ∼42% for a random loose packing of spherical, monosized particles (Scott, 1960), and values down to 10% can be reached (Dullien, 1991).
In binary mixtures, the way particles interact depends on their size ratio, here defined as the ratio of the respective particles' volume-equivalent diameters. If this ratio is less than 0.154 (Graton & Fraser, 1935), small particles will not affect the packing state and simply fill the gaps between larger ones. In contrast to this unmixing of particles, similar sized particles will mix, creating a new packing structure (also see A. Yu & Zou (1998) for a discussion of mixing and unmixing effects). Applying these concepts to polydisperse mixtures, particle unmixing will take place for very small and very large particles, as smaller particles start filling the gaps and larger particles completely fill some regions with solid material. The component controlling the porosity of the mixture is then defined by intermediate-sized particles, which do not change their packing state by the addition of unmixing components (A. Yu & Zou, 1998). An illustration of particle mixing and unmixing is shown in Figure 1. The semiempirical models by A. Yu & Zou (1998) and Zou et al. (2011) and can be applied to particle mixtures in loose and dense packing states. They have been shown to reproduce the porosity of mixtures created using the funnel method, in which particles are gently poured into a container, as well as the porosity of mixtures tapped many times to reach maximum compaction.
It is important to note that packing is determined by the interplay of the different grain sizes present, and it can be misleading to consider individual grain sizes only. For example, while the addition of a single large block to the mixture can reduce porosity by displacing smaller particles and filling void spaces, the addition of many large blocks can increase porosity by creating large voids. Similarly, addition of some small particles may reduce porosity, while many small particles can create a large number of small voids, again increasing porosity. Therefore, the porosity finally attained by the mixture depends on the details of the size-frequency distribution (SFD) of the particles present.
In order to apply the theory of multicomponent mixtures, the size and shape distributions of boulders need to be known. Here we use the boulder size and shape distributions determined by Michikami et al. (2019), who extend the analysis in Sugita et al. (2019) using images from the Hayabusa2 optical navigation camera (Kameda et al., 2017;Suzuki et al., 2018;Tatsumi et al., 2019) which have near global coverage and were acquired at altitudes between 20 km and 6.5 km. These have spatial resolutions down to 0.65 m/pixel, and global counts were performed for boulders with diameters >2 m and a completeness limit of 5 m. In addition, smaller boulders, cobbles, and pebbles with sizes of 0.02-9.1 m were studied using close-up images of the sampling areas, where images taken at altitudes from 67 m to 620 m with resolutions down to <0.01 m/pixel are available (Michikami et al., 2019). Overall, size-frequency and shape distributions were determined in the 0.02-140 m size range. experiments on the disruption of monoliths (Michikami et al., 2016), which suggest that boulders on bodies such as Itokawa, Bennu, and Ryugu are relicts of the direct formation of those asteroids by gravitational reaccumulation following the disruption of their parent bodies (Michel et al., 2020;Michel & Richardson, 2013) rather than the result of impact events after formation has been completed. Impacts could reshape the size distribution by the production of smaller particles after reaccretion has been completed, but the importance of this process may be limited. This is due to the so-called armoring effect (Sugita et al., 2019), by which a large fraction of the impact energy is lost when the projectile contacts the first large boulder, thus producing only few fragments. Another mechanism that could be responsible for a difference between the SFDs observed on the surface and present in the interior is seismic shaking, and the Brazil Nut Effect could lead to an overrepresentation of large boulders on the surface (Maurel et al., 2017;Tancredi et al., 2015). However, the seismic efficiency of impacts in granular material appears to be low (Nishiyama et al., 2020;Yasui et al., 2019), such that surface modifications are likely localized. Nevertheless, seismic shaking could have an impact on the global boulder SFD over geological timescales. Finally, it has been argued that particle size sorting may take place during rubble pile reaccretion, with larger blocks accreting first and thus in the center (Britt & Consolmagno, 2001). These caveats need to be kept in mind when interpreting the results presented below.
A second important input parameter for the multicomponent mixing model is the material's packing state, which can vary between a random loose and random close packing. In general, little is known about the packing state of rubble pile asteroids following reaccretion, which depends on many parameters such as the distribution of angular momentum in the reaccreting system as well as the size distribution and shape GROTT ET AL.
3 of 15 10.1029/2020JE006519 Figure 1. Top left: Illustration of particle unmixing for particles with strongly disparate diameters. As small particles are added to a system of larger particles, the larger particles resist being displaced and the packing state does not change. A similar effect occurs for the addition of very large particles. Bottom left: Illustration of particle mixing for particles with similar diameters. As similar sized particles are added to a system, particles can be displaced thus changing the packing state. Right: Two dimensional illustration of the random packing structure of strongly polydisperse spheres. As compared to monodispersed configurations, porosity is reduced by the filling of void spaces. Macroporosity refers to the porosity generated by the void spaces between particles, while microporosity is caused by void spaces and cracks that formed inside individual particles. Figure adapted from Yu and Zou (1998).
of reaccreting fragments. While impact experiments indicate that shattered, elongated particles with large deviations from a spherical shape can be produced (Durda et al., 2015;Michikami et al., 2016;Nakamura & Fujiwara, 1991) the results of disruption experiments need to be interpreted with caution, as the high strain rates imposed during the experiment may not be representative for the destruction of larger blocks. Furthermore, long term seismic shaking could lead to the reduction of pore spaces. Given these unknowns, we will systematically vary the packing state in the analysis below.
In the following, we will first introduce the theory of determining asteroid macroporosity from observed size and shape distributions for the rocks and boulders. We will then derive a simple equation relating grain density to macro-and microporosity. Results of the macroporosity calculation and relevant uncertainties will then be used to estimate grain density of Ryugu's constituent material given estimates of boulder microporosity Hamm et al., 2020;Okada et al., 2020). Finally, results, assumptions, and implications will be discussed.

Particle Size and Shape Distributions
To estimate Ryugu's macroporosity, the constituent boulder's size and shape distributions need to be known. These were determined by who fitted size-frequency data using power laws. Power law exponents between 1.65 and 2.65 were obtained, with 2.65 being the best fit for the global data set. Furthermore, particles were generally found to be elongated, and axis ratios for boulders >2 m are close to 0.7 on average. The SFD of boulders on small bodies may better be described by a Weibull distribution than a power law (Schröder et al., 2020), and we have used a cumulative Weibull (Rosin-Rammler) distribution (Brown & Wohletz, 1995;Rosin, 1933;Weibull, 1951;Wingo, 1989) to represent the data provided by Michikami et al. (2019). The cumulative SFD N(D) is then given by where D is the mean horizontal diameter, and we determined the fit parameters β = 0.09495, λ = 33.78 m, and N T = 5.28 × 10 14 km −2 by a weighted least-squares approach as a practical means to obtain a good representation of the data. The resulting distribution N(D) is shown together with the uncertainty of the data in Figure 2, where uncertainty comprises the Poisson uncertainty as well as the uncertainty of particle diameters introduced by the limited image resolution. It is worth noting that representing the data using a single power law for the entire size range does not adequately represent the data.  Given the SFD N(D) as determined from surface counts of boulders, the normalized cumulative volume distribution V(D) can be calculated by numerical integration. It is given by where N tot is the total number of particles counted per unit area, D min and D max are the minimum and maximum particle sizes of the particle size distribution N(D), respectively, and c (m −1 ) is a normalization factor chosen such that V(D max ) = 1.
In addition to the Weibull distribution fit to the data represented by Equation 1, we will also consider a simple power law to systematically study the influence of the particle size distribution's power law exponent p on the obtained results. The distribution can then be expressed as where p = 2.65 represents the best fit to the global data set (Michikami et al., 2019). For the power law defined by Equation 3, Equation 2 can be integrated analytically and the volume-size distribution is then given by Michikami et al. (2019) give the shape of boulders in terms of the maximum dimensions in three mutually orthogonal planes (a ≥ b ≥ c). Here we primarily regard the horizontal axis ratio b/a, with a being the maximum and b the intermediate dimension. As reported by Michikami et al. (2019), shape of particles on Ryugu appears to be largely independent of geographical longitude, whereas some dependence on latitude may indicate boulder migration. Nevertheless, average b/a is only weakly size-dependent and close to 0.7.
In general, particle sphericity is defined as the ratio of the surface area of a sphere (with the same volume as the particle) to the surface area of the particle (Wadell, 1932). However, this is difficult to evaluate in practice, and the Krumbein (1941) or Riley (1941) simplifications are usually applied. Working with two dimensional (image) data, we define sphericity Ψ as where D i is the diameter of the largest inscribed circle and D c is the diameter of the smallest circumscribing circle for a given particle (Riley, 1941). Using Equation 5, the shape parameter b/a derived by then translates into an average sphericity of Ψ = 0.83. In addition, Michikami et al. (2019) also estimated the third axis, c/a, of 121 arbitrarily selected boulders. The mean axes ratio c/a was found to be 0.44, and the sphericity of a parallelepiped with axis ratios a:b:c of 1:0.71:0.44 is 0.796. On the other hand, sphericity of a triaxial ellipsoid with the same axis ratios is 0.913. Therefore, sphericity depends not only on axis ratios, but also on particle shape, and we will use Ψ = 0.85 ± 0.06 as an average sphericity rather than the average sphericity derived from the shape data in Michikami et al. (2019) when calculating interparticle forces and initial porosities below.

Macroporosity
The macroporosity of Ryugu can be calculated from the volume SFD (Equation 2) assuming linear mixing models (A. Yu & Zou, 1998;Zou et al., 2011). In the mixing theory, the macroporosity achieved for a given size distribution will be a function of the volume fractions X i , the initial porosity ϕ i , as well as the nominal equivalent volume diameter d i of particles in each bin. The latter represents the diameter of a volume-equivalent sphere. Furthermore, i = 1,…,n is the number of size bins used and d 1 > d 2 > … d n for convenience.
Note that the equivalent volume diameter d i of particles is not strictly identical to the mean horizontal diameter as defined by Michikami et al. (2019), but as the observed boulder axis ratios on Ryugu change only little as a function of horizontal diameter, the shape factor relating horizontal diameter to the equivalent volume diameter d i is close to constant. It can thus be factored out for the mixing model below and has a negligible effect on the Bond number.
The above formulation holds if particle sphericity is independent of particle size, which is the assumption made in the following. However, we note for completeness that the method to estimate macroporosity used here can be generalized to arbitrary sphericity-size relations Ψ(d) by introducing the equivalent packing diameter d p , which then accounts for particle shape effects, that is, mixing of particles that have different sphericities at different sizes. Then, the equivalent volume diameter d in Equation 6 needs to be replaced by the equivalent packing diameter d p , which is related to the observed equivalent volume diameter d through sphericity Ψ(d) by (A. Yu & Zou, 1998) by the empirical relation The dimensionless specific volume describing the packing state for each bin is defined as (Zou et al., 2011) and the macroporosity finally attained by the mixture will be governed by the interaction of all differently sized particles. However, there will be one intermediate-sized bin i that controls the packing structure (see A. Yu & Zou (1998), also compare Figure 1). While the size bin number i of the controlling component is not known a priori, the specific volume  V i of a particular packing can in general be expressed as where small particles have indices j = 1 … i − 1 and large particles have indices j = i + 1 … n. The functions f(d i ,d j ) and g(d i ,d j ) are referred to as interaction functions between components i and j and were derived experimentally (A. Yu et al., 1997;Zou et al., 2011). They are given by and depend on the equivalent packing diameter size ratios r ij between small and large particles of the two components. Parameters r ij can be expressed as (Zou et al., 2011) where R ij = d j /d i is the small-to-large size ratio and i < j. The empirical parameter k is 0.451 (Zou & Yu, 1996), and x ij depends on the type of particle-particle interaction (Zou et al., 2011). It is given by cri cri 0.697 / cri cri 1 0 1 1.543 In the above equation, the critical particle diameter d cri divides fine and coarse particles, that is, it is the particle diameter below which cohesion between particles starts to influence particle interactions. Under Earth gravity conditions, d cri is close to 150 μm (Zou et al., 2011), but under microgravity conditions, cohesion GROTT ET AL.
6 of 15 10.1029/2020JE006519 can be relevant even for decimeter-sizes boulders (Kiuchi & Nakamura, 2015;Scheeres et al., 2010;Zou et al., 2011). Here, we define the critical diameter based on the Bond number B, that is, the ratio between interparticle forces and the weight of a particle (Scheeres et al., 2010).
We define the Bond number assuming a cleanliness factor equal to unity and a particle separation of 1.5 × 10 −10 m (Scheeres et al., 2010). Furthermore, we calculate the cohesive force for equally sized particles and include effects of particle sphericity Ψ and roundness Ω (Powers, 1953) by adding these as multiplicative factors (Wood, 2020). The Bond number is then given by (14) where d is particle diameter, g = 0.9825 × 10 −4 m s −2 is volume averaged gravity of Ryugu (Yamamoto et al., 2020), and A = 4.1 × 10 −20 J is the Hamaker constant for olivine in high vacuum (Perko et al., 2001). While olivine is certainly not the most common mineral in carbonaceous material, we consider its Hamaker constant to be a more appropriate choice than, for example, the widely used Hamaker constant for amorphous SiO 2 . In any case, the Hamaker constant needs to be regarded as highly uncertain. This also implies that the exact choice of parameters like boulder density, sphericity, and roundness has little influence on the results presented below. We choose boulder bulk density ρ = 1420 kg m −3 to match a macroporosity of 16% and a bulk density of 1,190 kg m −3 (Watanabe et al., 2019) for consistency, where ρ was determined using an iterative approach. Furthermore, we choose a particle roundness Ω of 0.24, as appropriate for angular to subangular particles (Powers, 1953).
The resulting Bond number for parameters appropriate for Ryugu is shown in Figure 3 as a function of particle diameter. The critical diameter d cri corresponding to a critical Bond number B cri = 0.1 is indicated in blue and has been calculated using Equation 14. We use B cri = 0.1 as a baseline, that is, we assume that cohesion starts to have a noticeable effect on porosity once the interparticle forces exceed 10% of the particle weight. For Ryugu, B cri = 0.1 corresponds to d cri = 0.52 m, but the influence of varying B cri over a large range will also be discussed.
To evaluate Equation 9, we first discretize the size range between D min = 0.02 m and D max = 140 m into log(D max /D min )/log(q) logarithmically spaced bins. We use a size factor of q = 1.05 from one bin to the next, resulting in a total of 182 size bins, which turned out to be sufficient. Volume fractions X i in each size bin were calculated according to the Weibull or power law representation of the SFD as needed. Furthermore, initial specific volumes V i and therefore initial porosities ϕ i need to be prescribed. While initial porosities of coarse monosized spherical particles generally vary between 0.42 for loose random packing and 0.36 for dense random packing (Scott, 1960), cohesive forces between small particles can considerably increase porosities (Kiuchi & Nakamura, 2015;Scheeres et al., 2010). We use the empirical relation (Kiuchi & Nakamura, 2015) to determine initial porosity, where ϕ 0 is the porosity of the noncohesive particles and describes the packing state. Note that we here implicitly assume initial porosities as appropriate for spherical particles, and as for the relevant range of observed sphericities, the influence of deviations from an ideal spherical shape on initial porosity is negligible (Zou & Yu, 1996). Particle shape enters Equation 15 in the Bond number B(d i ) only, and it is a secondary effect in the analysis presented for Ryugu below. The constants α = 2.414 and γ = 0.1985 have been derived from a new fit to the data of Kiuchi and Nakamura (2015). Finally, the specific volume occupied by the mixture is obtained by calculating the maximum of all specific volumes for the different controlling mixture sizes and  Figure 3. Bond number, that is, the ratio between interparticle forces and particle weight, as a function of particle diameter, assuming parameters as appropriate for Ryugu. The diameter corresponding to a critical Bond number of B cri = 0.1 is indicated. Mixture macroporosity is then given by ϕ Macro = 1 − 1/V.
In summary, the following steps need to be performed to determine the macroporosity of a granular mixture using the model above: First, volume fractions in the individual size bins need to be calculated from the given SFD (Equations 8 and 9). Then, initial porosity in each size bin needs to be determined. This will primarily depend on the packing state. Furthermore, it also depends on particle roundness and shape, which influence cohesion (Equations 14 and 15) as well as the geometrical packing properties (not considered here). Finally, the macroporosity is determined by examining all possible particle interactions (Equation 16).

Average Grain Density
While the main goal of the present paper is a determination of the macroporosity of rubble pile asteroid Ryugu, additional information on the asteroid's average grain density can be derived. As macroporosity ϕ Macro , microporosity ϕ Micro , and bulk density ρ Bulk are related by information on grain density ρ Grain can be extracted from Equation 18 requires the macroporosity, microporosity, and bulk density to be known. While the bulk density of Ryugu was estimated to be 1190 ± 20 kg m −3 (Watanabe et al., 2019), the boulders' microporosity cannot currently be unambiguously constrained due to the difficulties associated with extrapolating meteorite thermal conductivities to porosities in excess of 20% Macke et al., 2011). However, end-member models (Flynn et al., 2018;Henke et al., 2016) suggest microporosities ϕ Micro of either 32% ± 2% or 50% ± 2% for Ryugu's dark and rugged boulders (Hamm et al., 2020) which comprise the vast majority of all boulders observed on the surface (Okada et al., 2020;Sugita et al., 2019). We will use Monte-Carlo simulations to propagate these uncertainties to the determination of Ryugu's grain density, while simultaneously taking the uncertainty associated with Ryugu's macroporosity as derived from the linear mixing theory (Section 2.2) into account.

Results
Given the parameterization of the SFD (Equation 1) for the boulders observed on the surface of Ryugu, and assuming the distribution also applies to the interior, we have first calculated the corresponding volume frequency distribution using Equation 2. Given roundness Ω, Hamaker constant A, particle bulk density ρ, and volume average gravity g (see Equation 14), we then varied the initial porosity ϕ i in each size bin (Equation 15) using a Gaussian distribution for ϕ 0 centered around 39.5% with standard deviation of 3%. In addition, particle sphericity was varied using a Gaussian distribution centered around 0.85 with standard deviation of 0.06, and 10 6 draws from these distributions were used in a Monte-Carlo simulation to calculate the resulting macroporosity according to Equation 16.
Results of the calculation are shown in the left-hand panel of Figure 4, where a histogram of the obtained macroporosities ϕ M is shown. The range of macroporosities obtained in the calculations is ϕM = 16.2% ± 2.6% (1-sigma), and thus considerably smaller than porosities of monodisperse packings. This is not surprising given the broad particle size distribution observed on the surface of Ryugu.
Given the range of macroporosities derived above as well as estimates for the boulder microporosities derived from in situ thermal inertia measurements (Grott et al., 2017Hamm et al., 2020), we calculated the range of grain densities compatible with the observed bulk porosity of Ryugu (Watanabe et al., 2019) using Equations 17 and 18. We applied two end-member models for the microporosity ϕ Micro : for the first model (Flynn et al., 2018) we use ϕ Micro = 50% ± 2%, while for the second model (Henke et al., 2016)  using Gaussian distributions centered around 50% and 32% with standard deviations of 2%, respectively. Furthermore, we varied bulk density using a Gaussian distribution centered around 1,190 kg m −3 with a standard deviation of 20 kg m −3 (Watanabe et al., 2019) and macroporosity using a Gaussian distribution centered around 16.2% with a standard deviation of 2.6%.
Results of the calculation are shown in the right hand panel of Figure 4, where the resulting histograms for the grain densities ρ Grain are shown for the two end-member models. Owing to the two different models used to estimate boulder microporosity, two separate peaks are obtained for the distribution of grain densities.  Hamm et al., 2020), and two end-member models for the microporosity have been assumed. asteroids. For smaller power law exponents, the SFD is shallower as compared to distributions with larger p, and as a result, such distributions represent surfaces with a higher ratio of large particles.
In general, the macroporosities ϕ Macro obtained using the above mixing theory show a distinct minimum at intermediate power law exponents p, whereas distributions which have too many small or too many large particles result in unfavorable mixing and larger ϕ Macro are obtained. This minimum around p = 2.5 is known as the Fuller parabola in the engineering literature and has long been known as the optimum packing size distribution for spherical particles (Fuller & Thompson, 1907). Results obtained varying the critical Bond number are shown in the left panel of Figure 5, where the critical Bond number parametrizes the particle size below which interparticle forces result in significant cohesion. As expected, low critical Bond numbers, corresponding to larger contributions from cohesive particles, result in larger macroporosities. However, the overall effect is small and in the few percent range. The low critical Bond number of 0.1 adopted above therefore results in a conservative upper limit on macroporosity. It is also worth noting that results obtained using a power law distribution with p = 2.65, which overestimates the fraction of small particles, are lower than those obtained using the Weibull representation of the data by 4%-5%, such that results obtained using global power law fits must be interpreted with caution. For comparison, results obtained neglecting cohesion are also shown in the left panel of Figure 5, and macroporosity approaches a limit of 39.5% (compare Equation 15) for large p (not shown).
The influence of varying the lower cutoff diameter D min of the SFD on the obtained macroporosity ϕ Macro is shown in the right panel of Figure 5, where ϕ Macro is shown as a function of D min for three power law exponents p. In the calculations, a critical Bond number of B cri = 0.1 has been assumed. While the minimum macroporosity that can be achieved by the packing is close to constant for small D min , predicted macroporosity drastically increases for cutoff diameters larger than a few decimeters. In this case, unfavorable mixing is a result of the sparsity of smaller rocks to fill the gaps between larger blocks. These results indicate that image data with centimeter resolution are necessary to properly characterize the packing state of rubble pile asteroids, and that results presented above are largely independent of the cutoff size of D min = 0.02 m imposed by the image data available for Ryugu.

Discussion and Conclusions
In the present study, we have used semiempirical models for the porosity of multicomponent mixtures to estimate the macroporosity of Cb-type asteroid (162173) Ryugu based on the observed SFD of boulders on the asteroid's surface and the assumption that the surface distribution of boulders is representative for the bulk asteroid. Using the concept of controlling mixtures (A. B. Yu & Standish, 1991;A. Yu & Zou, 1998;Zou et al., 2011), we estimated the macroporosity of Ryugu to be ϕ M = 16.2% ± 2.6%. Based on estimates of boulder microporosity, we furthermore constrained the average grain density of Ryugu's boulders to ρ Grain = 2848 ± 152 kg m −3 or ρ Grain = 2093 ± 96 kg m −3 , depending on the microporosity model used.
Boulder shape can affect the above mixing model by changing interparticle cohesion, by changing the geometrical arrangement between different particle sizes, and by changing the initial porosity in each individual size bin. In the modeling, we have taken the influence of shape on particle cohesion explicitly into account, while we neglected its influence on geometrical interactions and initial porosity. This is justified because for the case of Ryugu the majority of particles have axis ratios b/a in excess of 0.5 (Michikami et al., 2019), corresponding to sphericities larger than 0.7. For such particles, initial porosity is nearly independent of shape and equal to the value appropriate for spherical particles (Zou & Yu, 1996). It is also worth noting that for Ryugu all of the above are secondary effects when compared to the unknown packing state, which we address by considering the entire range stretching from a random loose to a random close packing.
For the case of Ryugu, the primary factor determining macroporosity is the boulder SFD, and while the applied model takes cohesion between particles into account, disregarding cohesion results in only a slight modification of the obtained macroporosity ϕ M . Switching off cohesion in the model by assuming a critical bond number of 10 5 results in a macroporosity of 16.1%, only 0.1% smaller than the value presented above. This is a direct consequence of the low volume fraction of small cohesive particles on Ryugu, which directly follows from the given boulder SFD. This also implies that the results presented here are robust with respect to the exact choice of parameters like boulder density, roundness, and sphericity, which enter the calculation of the Bond number. It is also worth noting that the power law representation of the data significantly overestimates the influence of cohesion on macroporosity when compared to the Weibull fit by overestimating the volume fraction of small particles.
The cratering experiment performed by Hayabusa2's small carry-on impactor resulted in the formation of a crater in the gravity-dominated regime , indicating that particle cohesion played a minor role in the crater formation process. On the other hand, particles with diameters of 0.2 m were observed in the SCI crater wall, which, according to Equation 14, have Bond numbers close to unity and should therefore interact cohesively. This apparent discrepancy is resolved by the fact that small particles appear not to be volumetrically dominant inside Ryugu. This is indicated by the shallow particle size distribution for particles smaller than 1 m on the surface (Michikami et al., 2019;Sugita et al., 2019) and inside the artificial crater   Figure S5), where the particle size distribution shows a power law exponent p ∼ 2 (also compare the volume-size distribution on the right hand side of Figure 2). Therefore, results of the cratering experiment confirm that cohesion has a small influence on Ryugu's packing state. However, cohesion may become significant for rubble pile asteroids with a steep particle size distribution, for example, power laws with p > 3, where-in contrast to Ryugu-the mixture is dominated by a high volume fraction of very small particles.
Although a full analysis using empirical fits of the cumulative boulder SFD of other small bodies has not been performed here, macroporosity results can be qualitatively compared by considering the power law exponents of their respective size distributions and assuming similar size cutoffs D min and D max . The former have been widely used to describe size distributions in the literature, and values of p = 2.9 ± 0.3 and p = 3.52 ± 0.20 have been obtained for Bennu (Lauretta et al., 2019) and Itokawa (Mazrouei et al., 2014;Michikami et al., 2008Michikami et al., , 2019, respectively. Assuming B cri = 0.1 as above, these correspond to macroporosities between 10% and 38% for Bennu and 43%-52% for Itokawa. Assuming average grain densities of 2,600 kg m −3 , Bennu's low bulk density of 1190±13 kg m −3 (Lauretta et al., 2019) implies a bulk porosity of 54%, indicating significant microporosity. For Itokawa, average grain density has been estimated based on the modal abundance of minerals in the returned samples, and densities of 3,400 kg m −3 have been obtained (Tsuchiyama et al., 2011(Tsuchiyama et al., , 2014. This implies a bulk porosity of 39% ± 6% (Abe et al., 2006;Fujiwara et al., 2006;Tsuchiyama et al., 2011), consistent with the results obtained from the mixing theory. For Eros, the power law exponent of p = 3.31 ± 0.06 (Thomas et al., 2001) implies macroporosities of 40%-45%, which is larger than the inferred bulk porosity of Eros. The latter is estimated to be 21%-33% (Wilkison et al., 2002;Yeomans et al., 2000), indicating that macroporosities derived using the presented mixing model are incompatible with the observations. However, although Eros is a heavily fractured body, there is little evidence that it was ever catastrophically disrupted and later reaccumulated into a rubble pile (Wilkison et al., 2002), such that the theory presented here can probably not be applied.
Results for Ryugu have been obtained assuming minimum and maximum particle sizes of 0.02 m and 140 m, respectively, and these results are robust with respect to the cutoff at small particle sizes D min . Only shifting the cutoff D min to values larger than 0.30 m has a noticeable effect on the macroporosity. The upper cutoff size D max was chosen to correspond to the Otohime boulder, which is the largest boulder observed on Ryugu's surface. However, boulders larger than Otohime could potentially reside in Ryugu's interior, which would decrease the obtained macroporosity through a filling of void spaces. Reasonable upper limits on monolith sizes are 200 m, as derived from observations of fast rotators in the asteroid population (Pravec & Harris, 2000) and the catastrophic disruption threshold (Benz & Asphaug, 1999;Jutzi et al., 2010). Assuming D max = 200 m reduces ϕ Macro to 15%.
One way to increase macroporosity in the above models would be an increased initial porosity in each size bin, which may for example be caused by mechanical interlocking of particles due to particle angularity. For a random loose packing, noncohesive initial porosity can increase from ∼42% for smooth frictionless particles to ∼44% for very rough particles (Jerkins et al., 2008;Onoda & Liniger, 1990). In the frame of the applied mixing model, this effect is taken into account in the chosen initial porosity (Equation 15) and shifting the applied Gaussian distribution in the performed Monte-Carlo simulations by 2% results in slightly increased macroporosities of 18.0% ± 3%. Therefore, while roughness and particle interlocking can increase macroporosity, this is likely not a significant effect.
While the obtained macroporosity may appear to be relatively low, a significant reduction with respect to the porosity of random close packings of monodisperse spheres can be expected. Even binary mixtures of particles can be arranged in packing states with porosities of 15%-20% (A. Yu et al., 1992;A. B. Yu & Standish, 1991), such that it should not be surprising to achieve similar packing densities with the broad size distributions used here. Ternary mixtures can achieve ϕ Macro <10% (A. B. Yu & Standish, 1991), and while most common loose or compact granular materials have macroporosities between 30% and 50%, almost any degree of macroporosity between 10% and 90% can be obtained for polydisperse angular particles (Dullien, 1991). Experimentally, macroporosities down to 10% have been produced in the lab (Latham et al., 2002). Therefore, the macroporosity of Ryugu obtained here falls within a reasonable range, and Ryugu's high bulk porosity is a direct consequence of the very large microporosity of Ryugu's boulders.
The average grain densities obtained here are much lower than typical grain densities of ordinary chondrites, which range from 3,520 to 3,710 kg m −3 (Flynn et al., 2018), and also lower than those of most carbonaceous chondrites, which typically have grain densities in excess of 3,360 kg m −3 (Flynn et al., 2018).
Only the CM and CI subclasses show lower grain densities, and ρ CM,Grain = 2960 ± 40 kg m −3 while ρ CI,Grain = 2420 kg m −3 (Consolmagno et al., 2008;Flynn et al., 2018;Macke et al., 2011). The Tagish Lake meteorite, an ungrouped carbonaceous chondrite, exhibits similar grain densities in the range between 2,430 and 2,840 kg m −3 (Ralchenko et al., 2014). While the larger grain densities of 2848 ± 152 kg m −3 are consistent with the CM and Tagish Lake results, the lower densities of 2093 ± 95 kg m −3 are inconsistent with those of known meteorite samples.
Estimates of grain densities discussed above indicate that extrapolating boulder porosities as a function of thermal conductivity using the model by Flynn et al. (2018) is preferred to extrapolations using the model by Henke et al. (2016). In addition, laboratory measurements of thermal conductivity  using the UTPS Tagish Lake meteorite simulant (Miyamoto et al., 2018) provide further evidence of high boulder microporosity. The UTPS simulant has a grain density of 2,813 kg m −3 and a porosity of 47.5%, while thermal conductivity was determined to be similar to that of Ryugu's rugged boulders . It therefore seems likely that boulder porosity on Ryugu falls within the high range determined by Grott et al. (2019), but more laboratory measurements of thermal conductivity at high porosity are needed to confirm these results and reduce uncertainties. If grain densities are indeed of the order of 2,850 kg m −3 , Ryugu's bulk porosity is estimated to be 58% (cf. Equation 18).
It is noted that close-up images have revealed that many boulders on Ryugu and Bennu exhibit morphologic properties consistent with a brecciated structure (Sugita et al., 2019;Walsh et al., 2019). Breccia would have much larger microporosities than pristine rocks, consistent with the large microporosities preferred here. Furthermore, the presence of breccia on Ryugu and Bennu is consistent with the fact that many carbonaceous chondrites and, in particular, all CM and CI meteorites found on Earth are known to be brecciated (Bischoff et al., 2006). However, it remains to be investigated if brecciation is the main mechanism providing microporosity or whether the boulder's highly porous structure is a result of the formation mechanisms acting in Ryugu's parent body (Neumann et al., 2014(Neumann et al., , 2015. If microporosity in typical carbonaceous asteroids is as high as predicted here for Ryugu, macroporosities of rubble pile asteroids may have been systematically overestimated in the past (e.g., Consolmagno et al., 2008). Macroporosities have been estimated based on measurements of asteroid bulk density and porosities of meteorite samples, the latter of which could have been underestimated compared to values for actual carbonaceous material on asteroids derived from in situ measurements Consolmagno et al., 2008). This bias could be the result of filtering by the Earth's atmosphere, as only the strongest, densest carbonaceous meteoroids would survive atmospheric entry, while weaker samples would break up (Grott et al., 2017. This could explain the absence of high porosity samples in our meteorite collections, where the most porous sample reported to date is the Tagish Lake meteorite, which shows porosities in the range from 26% to 36% (Popova et al., 2011). The samples to be returned from Ryugu by the Hayabusa2 mission will provide crucial information on this issue.
Results presented here assume that the SFD observed on the surface of Ryugu is representative for the entire asteroid, but as discussed in Section 1, the reaccretion process itself as well as post accretion surface modifications could influence the observed SFD. For example, meteorite impacts could increase the number of small boulders on the surface and the observed SFD would be steeper than the distribution in the interior. Therefore, macroporosity would have been overestimated in the presented model, as the interior distribution would move closer to the Fuller minimum (Fuller & Thompson, 1907). Conversely, the Brazil Nut Effect could bias the slope of the surface SFD toward smaller values, implying that macroporosity would have been underestimated. This topic can be addressed once average grain density and possibly also microporosity have been determined from the returned samples, as has been done for Itokawa (). Then, Ryugu's macroporosity can be derived given the measured bulk density (). Any significant deviation from the macroporosity value calculated here will indicate a nonhomogeneous boulder size distribution in the bulk volume of the asteroid.

Data Availability Statement
The numerical code and data have been made available