Quantitative grain size estimation on airless bodies from the negative polarization branch. II. Dawn mission targets (4) Vesta and (1) Ceres

Context. Sunlight scattered from the surface of an airless body is generally partially polarized, and the corresponding polarization state includes information about the scattering surface, such as albedo, surface grain sizes, composition, and taxonomic types. Aims. We conducted polarimetry of two large airless bodies, the Dawn mission targets (1) Ceres and (4) Vesta, in the near-infrared region. We further investigated the change in the polarimetric phase curves over the wavelengths expected from previous works. Methods. We used the Nishiharima Infrared Camera (NIC) installed at the Nishi-Harima Astronomical Observatory (NHAO) to observe these objects at multiple geometric configurations in the J, H, and $\mathrm{K_s}$ bands ($ \lambda \sim 1.2\mathrm{-}2.3 \mathrm{\mu m} $). Results. Polarimetric parameters were determined and compared with previously reported experimental results. In particular, Vesta exhibits a characteristic change in the negative polarization branch as the wavelength increases to the $\mathrm{K_s}$ band, which we interpret as an indication of the dominant existence of $D \sim 10\mathrm{-}20 \mathrm{\mu m}$ particles. Our approach is supported by empirical reasoning and coincides well with an independent, theory-driven approach based on thermal modeling. Conclusions. This work demonstrates how near-infrared polarimetry can be utilized to quantitatively determine the particle size of airless objects. This finding will have important implications for asteroid taxonomy and regolith evolution.


Introduction
The presence of a regolith on an airless body is the result of its initial physical and chemical conditions and specific history of evolution.Among the parameters used to characterize the regolith are grain size, which is related to thermal parameters (Gundlach & Blum 2013), and scattering processes (Hapke 2012), which can change the observable parameters.Thus, quantification of these parameters is one of the important steps in modeling observed quantities.However, the fine nature of grains, as seen with the submicron-scale fine particles retrieved from the Moon (Park et al. 2008), makes it difficult to quantify their sizes, even in this space exploration era.
In our previous work (Bach et al. 2024;Paper I hereafter), our extensive investigation established a methodology for quantifying grain sizes on airless bodies through multiwavelength polarimetry.We focused on exploring the widening-and-deepening (WD) trend within the negative polarization branch (NPB) as triggered by changes in albedo, A(λ), and the size parameter, X ∝ D/λ, where D represents grain size and λ denotes the observation wavelength.We identified two significant effects influencing the NPB across wavelengths.First, as A increases, we observe a WD-like trend for A ≲ 10% and a shallowing trend for A ≳ 10%.Second, regarding D/λ, we found that the NPB stays relatively unchanged for D/λ ≫ 10.However, as λ increases or D/λ decreases, the WD trend starts to appear at D/λ ≲ 5-10.By focusing on the effect of D/λ, we can estimate the particle size that contributes to the scattering process.
Although it is a useful approach, interpreting the WD trend remains challenging due to the potential for a convoluted influence of both albedo and D/λ.Thus, observations of (i) airless objects (ii) covered with fine grains and (iii) exhibiting relatively constant reflectance across visible to near-infrared (NIR) wavelengths are necessary for more reliable interpretations when focusing on the D/λ effect for cases free of the effect of albedo.
The particle size on asteroids is roughly related to the size of the parent asteroid; the larger the asteroid is, the smaller the particle size (Delbo' et al. 2007;however, see MacLennan & Emery 2022b,a).This is generally understood to be due to the stronger gravity of larger objects (Capria et al. 2014) and/or the shorter lifetime of smaller objects (Delbo et al. 2015).For (1) Ceres (size d ∼ 950 km), the presence of ultra-porous fine particles (D ≪ 1 µm) may be required to explain the observations (Li et al. 2019;Schröder et al. 2021;Sect. 4.2.2).Similarly, spectroscopic analyses suggest that (4) Vesta (d ∼ 530 km) is predominantly covered with particles of sizes ≲ 20 − 25 µm (Hiroi et al. 1994;Li et al. 2011;Martikainen et al. 2019), and optical polarimetry has indicated the possibility of a mixture of eucrite grains with sizes > 74 µm and < 10 µm Le Bertre & Zellner 1980; also see Sect. 4.2.1).
Therefore, in our study we chose to investigate Ceres and Vesta because (i) they are among the brightest airless objects, (ii) they are likely to possess fine particles that can display the D/λ effect (Paper I), (iii) both have reflectances that suggest they are not subject to a strong albedo effect (Paper I), and (iv) they are targets of the Dawn mission, so a large amount of information is available.
Our NIR polarimetry investigation of airless bodies appears to have commenced nearly simultaneously with that described in Masiero et al. (2022), where the first observations in the J and H bands from UT 2019-03-17 are reported.To the best of our knowledge, Masiero et al. (2022) is the first peer-reviewed publication of NIR polarimetry of airless bodies; it includes our own observations that have been conducted since UT 2019-10-22, at which time we were unaware of their pioneering work.In their work, the albedo and refractive indices, but not the grain size of asteroids, are discussed.In a later study (Masiero et al. 2023), they proposed a qualitative description of the existence of fine-grained inclusions on Barbarian asteroids.Both publications studied wavelengths up to the H band (λ ∼ 1.6 µm).We stress that our study is the first NIR polarimetry study of a Vtype asteroid (Vesta) and the first polarimetry study of airless bodies at λ > 2 µm (the K s band), a range that contains critical information (Sect. 4.2.1).
In this study we conducted a comprehensive investigation utilizing our new data and previous experimental results.The observational circumstances, data reduction methods, and polarimetric parameter finding methods used are summarized in Sect. 2. The results of our observations are presented in Sect.3. In Sect. 4 we discuss how our results can be interpreted and how the grain sizes of these targets can be obtained.Finally, the appendix provides details on our data reduction and statistical analyses as well as information regarding the data and code availability.
In the optical path, a polarimetric mask, a rotating half-wave plate (HWP), and a beam displacer are configured (Takahashi et al. 2018;Takahashi 2019).As a result, two images, ordinary and extraordinary rays (hereafter, o-and e-rays), appear on the same detector at a pixel scale of 0.16 ′′ /pix.Each image has a size of ∼ 150×430 pixels (24 ′′ ×69 ′′ ) and is displaced horizontally by approximately 200 pixels.Due to the beam displacer, the point spread functions are expected to be different for o-and e-rays Each set of polarimetric observations consists of four exposures for which with the HWP angle was set to 0 • , 45 • , 22.5 • , and 67.5 • , respectively.The four exposures together are denoted as a "set," with the first two and last two exposures denoted as the q and u subsets, respectively.These parameters are used to determine the normalized Stokes q = Q/I and u = U/I parameters, respectively.We conducted five-and six-night polarimetric observations for Ceres and Vesta, respectively.The observational conditions are summarized in Table 1.On UT2019-11-08, when Vesta was at α = 4.1 • , the weather conditions were very unfavorable.This circumstance led to widely varying signal-to-noise ratios (S/Ns) even between consecutive exposures.Consequently, the uncertainty, particularly in the H-band measurements, was substantially elevated during this period (Table 2).

Data reduction process
The data reduction process follows a similar approach to that of our previous publication (Appendix B of Takahashi et al. 2021), with some modifications implemented using the NICPolpy software1 (Bach et al. 2022).The instrumental gain and readout noise parameters were also updated from Ishiguro et al. (2011).Two polarimetric dome flats were used, one taken on UT 2018-05-07 and the other taken on UT 2020-06-03.For data taken after June 2020, the second master flat was utilized.
The preprocessing steps are as follows: First, artifacts such as vertical and Fourier-like patterns in the frames are removed.Then, standard dark subtraction and flat-fielding are applied.The values at bad pixels, identified using a pre-made bad pixel mask, are replaced with interpolated values from nearby pixels.Each frame is split into two FITS files, which represent the o-ray and e-ray images.Finally, cosmic-ray rejection is performed.Detailed explanations and the implemented algorithms are presented in Bach et al. (2022).
After preprocessing, signal extraction was performed for each exposure.The Stokes parameters were estimated using the ratio method (see, e.g., Bagnulo et al. 2009), considering each subset of exposures.To ensure the robustness of our polarization degree estimation, we extensively tested our results against multiple signal extraction techniques and outlier rejection schemes (see Appendix A).We also conducted thorough verifications to confirm the stability of our final results against different parameter values used in preprocessing, such as those related to artifact removal, the use of different flats, bad-pixel interpolation, and cosmic-ray rejection parameters.These methods only marginally changed the intermediate results, such that the final results (i.e., those describing the polarimetric parameters) do not vary beyond the measurement uncertainty range.

Polarization parameter calculation
The polarization degree of solar system objects is often described by the properly rotated polarization degree (Lyot 1929):  (a) Unstable weather (refer to the skymonitor record). (b) Most H-band frames are saturated. (c) All H-band frames are saturated.
where P is the total linear polarization degree, θ r is the position angle of the strongest electric vector relative to the normal vector of the scattering plane (e.g., 0 or 180 • indicates it is normal to the scattering plane), and I ⊥ and I ∥ denote the measured intensities of scattered light along the directions perpendicular and parallel to the scattering plane, respectively.The curve describing this phenomenon as a function of the phase angle (light source-target-observer angle; α 180 • − scattering angle), P r (α), is denoted the polarization phase curve (PPC).A PPC at small α can be expressed with a three-parameter linear-exponential function (Muinonen et al. 2002;Kaasalainen et al. 2003): for free parameters θ 1, 2, 3 .The inversion angle α 0 can be solved numerically from P r (α 0 ) = 0, as in Gil-Hutton (2017).We emphasize that the three parameters used above are eventually converted to polarimetric parameters of interest, for example the polarimetric slope h dP r dα α=α 0 = P ′ r (α = α 0 ), inversion angle α 0 , α min and P min = P r (α min ).By rearranging the equation, we can generate the fitting parameters with more intuitive meanings and facilitate subsequent statistical analyses of the parameters: where k [ • ] is a scale parameter.Having rearranged the equations such that two out of the three fitting parameters (h, α 0 ) are the ones of interest, the error analyses become much more straightforward.Importantly, the function in Eq. ( 3) describes the identical curve as in Eq. ( 2).Furthermore, we can obtain α min through the condition P ′ r (α = α min ) = 0 as follows: Therefore, four important parameters, h, α 0 , α min , and P min = P r (α min ), are obtained.

PPC of targets
The polarization degrees derived from our data are summarized in Table 2.These polarization degrees (P) are obtained using the formula P = q 2 + u 2 , where q and u are the instrumental polarization efficiency and the position-angle corrected values, respectively (Takahashi 2019).To calculate the uncertainties in q and u, we first computed the weighted average and its uncertainty, ( σ , where σ i is the uncertainty of the i-th data point.The final uncertainties, dq and du, were determined as the larger value recorded between the weighted average uncertainty and the standard error (i.e., the sample standard deviation divided by the square root of the number of samples).Because the position angle of the scattering plane at the time of observation has negligible uncertainty, dθ P = dθ r .For dP r , we applied error propagation to Equation (1), while neglecting the covariance between dP and dθ r .The maximum value between this propagated uncertainty and dP is obtained as dP r .
Notably, these uncertainties are comparable to the standard error from the central limit theorem, which is inherently much less than the random scatter (standard deviation of data points or the error bar of a single data point) by a factor of ≈ √ n.In their work, Masiero et al. (2023) accounted for this discrepancy by quadratically adding 0.1 %p scatter to their reported dP r values.Additionally, these uncertainties do not include systematic errors (see Appendix B).
As shown in Table 2, we have confirmed that our θ r values are generally close to integer multiples of 90 • , as expected from symmetry considerations (Lyot 1929;Cellino et al. 2018).Consequently, cos(2θ r ) ≈ 1.The extreme exception is Vesta at α = 22.5 • , where its polarization degree is near zero; thus, θ r is not well defined in this particular case.This finding strongly supports the credibility of our data reduction process.Conversely, we can further utilize our θ r values of asteroids to measure the instrumental offset (assuming it must be an integer multiple of 90 • ) when |P r | is large.For example, using Ceres data at α = 8.6 • , we find that the instrumental position angle offset is < 2 • for all the filters, which is consistent with Takahashi (2019).
The PPCs for Ceres and Vesta are derived and plotted in Fig. 1.For Ceres, two datasets from a previous publication were included (Masiero et al. 2022).We have overplotted these data in

Object
Filter a) 0.020 (a) 83.40 (a) 0.10 (a) Ceres (a)  J (a) − 21.0 − − +0.840 (a) 0.020 (a) 22.30 (a)  Notes.n: Number of used polarimetric datasets after data filtering and outlier rejections (one set corresponds to four exposures).α: Solar phase angle.P and dP: The total linear polarization degree and its uncertainty.θ P : The position angle of the strongest electric field vector.P r and dP r : The proper linear polarization degree referring to the scattering plan and its uncertainty.θ r and dθ r : Position angle of the strongest electric vector, with respect to the scattering plane normal vector, and its uncertainty.See Sect.3.1 for the details of error calculations. (a) Data from Masiero et al. (2022).
the figures, and the results of our observations closely match the corresponding fitting curves.The shaded regions in the figures represent the approximate 1σ range of the fitting models obtained from Monte Carlo (MC) simulations (Appendix C).
In Figure 2, we present the same data and curves as shown above, but with the inclusion of optical observations from Lupishko (2022) 2 and Bendjoya et al. (2022).For Vesta, we also used the data provided in Cellino et al. (2016), which were not included in previous datasets (at eight different αs), after taking 2 We used P r = −P if neither P r nor θ r was available since such data are from α ≪ α 0 .We did not use R2 filter (1/λ = 1.20 µm −1 ) data (2 observations for both targets were adopted from Zellner et al. 1974) because they are largely scattered from other data in the I filter (1/λ = 1.05 µm −1 ).Additionally, two data points in the 10 µm N band (Johnson et al. 1983) were mistakenly labeled as the 0.333 µm N band; therefore, we did not use these data.
the nightly median (3-sigma-clipped).When fitting the optical data, we used data only when the reported error bar for the polarization degree was < 0.15 %p; we neglected any uncertainties in the data because of the inhomogeneity in the data across multiple sources.Notably, we expect optical data to have higher uncertainties.For example, a few noisy Ceres data points near α 0 affect the h values (i.e., they have high leverage values).Finally, we excluded data points with > 0.25 %p deviation from the fitting line.We explored various processing methods, such as using only the data with available uncertainties, changing the error bar threshold, and/or adding a range of uncertainties (0.02-0.5 %p) in quadratures based on different criteria.The fitted polarimetric parameters change for these optical data; thus, the uncertainties of these parameters were difficult to quantify; however, our main results, discussed in Sect.4, were not affected by these variations.

Polarimetric parameters
The final fitting results are presented in Table 3.The best fit parameter values are the least-square solution, while the error bars are from the 16th and 84th percentiles of the MC samples.The MC median is very close to that of the least-square solutions, differing by only ≪ 1σ.For Ceres, we also included data from Masiero et al. (2022).More detailed information on the MC simulation is provided in Appendix C.
Several observations can be made from the abovementioned results: (i) the Ceres NPB does not appear to vary significantly with wavelength, and (ii) Vesta shows a clear trend of changing the NPB with wavelength, including a deepening trend in the K s band.The interpretations and discussions of these results are presented in Section 4.
In addition to the polarimetric parameters, the albedo value of the objects plays a crucial role in understanding their properties.To calculate the albedos of Ceres and Vesta, we utilized the catalog of visible to NIR reflectance spectra3 , which covers λ ∼ 0.45-2.5 µm, along with the filter profile information of the NIC.We performed trapezoidal numerical integration of the spectra using the filter profiles (but ignoring the reflectance and transmission of the dichroic mirrors) and then converted the integrated value back to the physical albedo.This process implicitly assumes that the spectrum used is the same at α = 0 • , implying that the backscattering amplitude is assumed to be identical over the wavelength range.The geometric albedo at a wavelength of λ = 0.55 µm (FC2 filter of the Dawn spacecraft) was adopted as p 0.55 µm = 0.094 ± 0.007 for Ceres (Ciarniello et al. 2017) and 0.38 ± 0.04 for Vesta (Li et al. 2013).The resulting albedo values and their uncertainties are summarized in Table 4.For the optical observations, we used a uniform filter profile spanning the wavelength range with 0.01 µm bezels (e.g., for 0.45-0.60,wavelengths from 0.45 + 0.01 = 0.46 µm to 0.60 − 0.01 = 0.59 µm).
For comparison, the figures in Paper I present Ceres and Vesta overlaid on various parameter spaces along with laboratory samples and other asteroids in the albedo-P min -h-α 0 space.Among these parameters, the Umov laws (h-albedo and P minalbedo relations) are shown in Figures 3-5 of Paper I. For optical wavelengths, rotational variations in amplitude < 0.1 %p and the Umov law were studied for Vesta (Cellino et al. 2016).However, we did not consider the rotational variation.
In the P min -α 0 space (Figure 7 in Paper I), the change in the location of Vesta, likely starting in the K s band, is noteworthy.Vesta is near S-complex asteroids (Belskaya et al. 2017) but moves toward the lunar fine region in the K s band.For the Moon, most regions exhibit a WD trend (because of the albedo effect and/or D/λ effect).A similar WD trend is also observed for Vesta.In contrast, the NPB shows almost no change in Ceres in the P min -α 0 space.

Discussion
4.1.The Umov laws (h, P min , and albedo) The h-albedo and P min -albedo relations are shown in Figures 3-5 of Paper I. First, the two objects in the NIR region roughly follow Umov laws.To the best of our knowledge, this is the first report on the location of small airless bodies in these parameter spaces at NIR wavelengths, even though further investigation is needed to verify its generality.Second, Vesta changes its location toward the fine particles (small D/λ) in both parameter spaces, which supports the existence of fine particles on Vesta.As mentioned, the Ceres h optical data are susceptible to a few noisy data points near the inversion.Ceres also shows a noticeable trend of increasing h (but not P min ) only in the K s band, while the albedo and P min remain nearly constant.

P min -α 0 space
The NPB of Vesta clearly shows a trend of WD increasing λ, as shown in, for example, Figure 6 of Paper I. However, Ceres stays at a fixed location over the investigated wavelengths (< 0.5 µm to > 2 µm).We therefore discuss the particle sizes on these objects based on these observations.

Determining Vesta's particle size
Vesta has albedo ≫ 10% all over the wavelengths of interest (Table 4), which allows us to conclude that WD of the NPB by the albedo effect (Paper I) is unlikely.Therefore, the likely cause is the effect of D/λ.We argue that D/λ ∼ 5-10 is reached in the H or K s band (λ ∼ 2 µm), so Vesta is likely covered with particles of size D ≳ 10-20 µm.Finer particles (∼ 1 µm) should not dominate the surface; if they do, that indicates that we must have detected a WD trend (or a narrowing-and-deepening trend) in the UV-to-J band observations.This particle size is consistent with previous speculations.Most importantly, Le Bertre & Zellner (1980) was one of the first studies to investigate optical NPB to determine the grain size of Vesta.After testing multiple Howardite, eucrite, and diogenite (HED) meteoritic samples, ≲ 20 µm particles were found to be necessary to reproduce the observed α 0 of Vesta.However, finer particles resulted in too high an albedo4 compared to Vesta's albedo known at the time (0.22;Morrison 1977;Zellner et al. 1977).Because of this mismatch, coarse particles are required to reduce the albedo.Finally, taking all these factors into account, they concluded that the eucrite Bereba sample of coarse particles (> 76 µm) covered with fine (< 10 µm) particles (denoted Bereba B) can match both the α 0 and albedo at the same time.
However, based on the Dawn mission, Vesta's geometric albedo is updated to p V = 0.42 (Li et al. 2013), and the phase function in the F2 filter (Figure 19(a) in Li et al. 2013) shows that the albedo at α = 5 • is ≳ 0.3, compared to earlier values 0.22.Thus, the logic that led Le Bertre & Zellner (1980) to a large fraction of coarse particles should be reconsidered.Instead, the updated data may support a larger fraction of finer particles, which is consistent with our results.
Later, Hiroi et al. (1994) used an independent method (reflectance spectrum analysis) and found that the finest particle size (< 25 µm) of HED powders best reproduced the Vesta 0.3-2.6 µm spectra.Hence, they also concluded that very fine (< 25 µm) particles are necessary to explain Vesta's reflectance spectrum.Recently, Li et al. (2011) reached a similar conclusion that a single, fine-grained Howardite with a size of < 25 µm suitably fits the Hubble Space Telescope, Swift, and International Ultraviolet Explorer (IUE) observations.In addition, Martikainen et al. (2019) showed that, using ray-optic simulation calculations, Vesta's surface is dominated (> 75%) by fine Howardite particles of size < 25 µm.
Finally, utilizing the thermal conduction modeling algorithm (MacLennan & Emery 2022b), we found two solutions for Vesta, D ≲ 10 µm and D ≳ 50 µm, both of which require a porosity ≳ 80%.Remarkably, this theoretical modeling, independent of polarimetry, yielded results that align quite closely with our findings.The slight disparity between our findings and those of the theoretical thermal model can be attributed to multiple factors, including that polarimetry detects only the topmost layer of the scattering process, whereas thermal modeling probes are significantly deeper (on the order of the thermal skin depth, which varies from the millimeter scale to as high as the meter scale).
Combining the lower bound from our study and the upper bounds from previous works, we propose that Vesta is covered with D ∼ 10-20 µm particles.

The case of Ceres
Because Ceres shows no change in the P min -α 0 space, one could argue that D/λ ≫ 10 until the K s band (D ≫ 20 µm).On the other hand, as noted in Paper I, the narrowing-and-deepening (ND) trend is not clearly detected across the literature, because of the inherent practical difficulty.Moreover, further experiments are needed to investigate when the ND trend stops and how this trend can differ at low albedos (A ≲ 10%).Therefore, another possible hypothesis is that the fact that D/λ ≪ 1-2, even at the shortest wavelength (i.e., D ≪ 0.5-1 µm), is the cause of Ceres's behavior.
From the Dawn observations, the wavelength dependence of the single-particle phase function is found (Li et al. 2019).They found that a considerable fraction of ≲ 1 µm particles can explain this dependence.Some regions near young impact craters on Ceres show blue spectra (Schröder et al. 2017).Later, Schröder et al. (2021) showed that these locations are likely to contain hyperporous, foam-like phyllosilicate structures, which effectively act as Rayleigh scatterers with a characteristic filament size ≪ 1 µm and filament separation ≫ 1 µm.These can support the small particle scenario.
Utilizing a thermal modeling algorithm similar to Vesta (MacLennan & Emery 2022b), we found two solutions, D ≲ 2 µm and D ≳ 30 µm, both with porosities > 80% (cross-checked with MacLennan, E. 2023, priv. comm.).We can reject neither the small nor large particle scenarios.However, we again stress that both possible scenarios from this work and thermal mod-eling, which are completely independent approaches, expect for surprisingly consistent particle sizes.
Because the NPB behavior of Ceres seems to be typical of C-complex objects (Masiero et al. 2022), further exploration of the change in the NPB for low-albedo materials is key for understanding C-complex objects.Moreover, these findings will naturally deepen our understanding of the thermal behavior of these materials by cross-validating both approaches.
Finally, we note that Ceres's location in the h-albedo space is changed in the K s band.This may indicate the existence of a fine structure.However, as shown in Paper I, the deviation for low-albedo (A ≲ 10%) materials has been investigated less than that for higher-albedo materials; thus, further investigations are needed.Belskaya et al. (1987) demonstrated similar trends in P min for both Vesta and ( 16) Psyche across UBVRI bands, despite differences in their reflectances.A later study showed that S-and M-types generally exhibit the opposite trend to that of the Umov law, which posits an inverse correlation between the polarization degree (|P min |) and albedo (Belskaya et al. 2009).This finding suggested that the Umov law is not the sole mechanism influencing the shape of the NPB.This trend aligns with the WD trend, while the trend for low-albedo objects (WD-like as albedo increases; Paper I) is also possible5 .This finding highlights the importance of Vesta because it has a sufficiently high albedo (Table 4) and allows us to rule out the presence of the albedo effect with a higher confidence.Bagnulo et al. (2015) mentioned that "the Umov law may be violated" based on the spectropolarimetry of the Barbarian asteroids.They found that (7) Iris and (599) Luisa showed that |P r (λ; α)| decreased with increasing reflectance, which matches the Umov law, while (236) Honoria showed the opposite trend over λ.They further discussed that the difference in polarimetric colors between Honoria and Luisa is puzzling because both are Barbarians; over these wavelengths, the NPB becomes deeper for Honoria and shallower for Luisa.

Revisiting several previous works
We emphasize that if the ND trend is present on these Barbarian asteroids, this is an expected behavior.Honoria and Luisa were observed at α = 7.1 • (where P r gets deeper with increasing λ) and α = 26.9• (slightly below α 0 , where P r gets shallower with increasing λ), respectively.Indeed, Masiero et al. (2023) found that L-types show an ND trend over the wavelength in the NIR region.
We agree with their interpretation that the NPB of L-type plants is affected by fine grains (D/λ < 1-2).This example hints that the Umov law should be applied with care; for example, one may use h or P max (P min may be less accurate).It is natural that the inverse Umov law or even the independence of the polarization degree from the wavelength may be observed at certain phase angles.Similarly, the change in NPB for S-types may indicate the presence of smaller particles on them, which may be related to the explanation in Hasegawa et al. (2019).A possible test observation would include NIR polarimetry of Q-types because they would not show a clear WD trend as S-types if the WD trend was indeed driven by fine particles.
In Kwon et al. (2023), an anticorrelation was found between α 0 and the asteroid diameter of the C-complex.These authors mentioned that this contradicts the fact that larger asteroids may contain finer particles.We hypothesize that this can also be interpreted as the ND trend being observed: the C-complex objects with diameters d > 100 km already possess fine particles that can induce the ND trend, so α 0 decreases as D decreases (i.e., a decreasing α 0 increases d).We admit that this should be considered a primitive hypothesis because there was no clear ND trend observed for a single object.

Conclusions and future work
In this work we have shown that NIR polarimetry can be used as a tool for quantifying grain size.We observed two of the largest airless bodies, (1) Ceres and (4) Vesta, which are Dawn mission targets.These airless objects were chosen because they are very likely to possess very fine particles on their surfaces; hence, D/λ is expected to be small enough to indicate the change in the NPB over NIR wavelengths.The data were processed with extreme care to achieve the best possible accuracy.Conclusions from this work are summarized as follows: 1.The Umov laws (slope-albedo and P min -albedo correlations) for airless bodies apparently hold at NIR wavelengths.At the longest wavelengths, however, they begin to deviate, which is attributed to the decreased D/λ. 2. Vesta showed a characteristic change in the NPB when λ was increased to the H and K s bands (λ ≳ 2 µm).Combining our data with those of previous studies, we find that the optically significant grain size on Vesta is ∼ 10-20 µm.This finding agrees with that of multiple previous independent reports.To the best of our knowledge, this report is the first quantification of grain size on airless bodies based on NIR polarimetry.3.For Ceres, there was no significant change in the NPB (P min , α 0 ).This indicates the potential existence of either very fine particles (D ≪ 1 µm) or large particles ≫ 20 µm). 4. Our empirical estimation of particle sizes using multiwavelength polarimetry strongly aligns with an independent theoretical modeling approach (thermal modeling).This successful congruence between the two methodologies represents a significant cross-validation of the two scientific approaches.
Because Vesta is known to have high heterogeneity in terms of albedo and composition (Reddy et al. 2012), disk-resolved polarimetry of Vesta near α min can reveal valuable hints about the parameters that change the NPB.A change in the polarimetric slope (h) of Ceres in the K s band may provide important information about its surface structure.A further experimental exploration of darker materials, for example using the samples from (101955) Bennu obtained by OSIRIS-REx (Lauretta et al. 2017) or (162173) Ryugu by Hayabusa 2 (Nakamura et al. 2023), may further advance our understanding of the case of dark C-complex and related objects.For example, if we assume the r ap0 values are 15 and 16 for the first frame and 17 and 18 for the second frame (o-and eray, respectively), under our notation r ap would be 25 for the 25 code, regardless of r ap0 values; r ap = 15 for min; r ap = 18 for max; r ap,o|e = 15|16 for minoe; and r ap,o|e = 17|18 for maxoe.
We find our final results to be extremely robust to all of these different aperture models (Figs.A.1 and A.2).For Vesta, (i) the inversion angle α 0 might be slightly uncertain, especially because there are no data at α > α 0 , and (ii) the error bar of the smallest α observation (UT2019-11-08) varies depending on the aperture mode.However, as seen from the figures, this does not change the logic of this paper.All the main results presented in this work are derived using the min aperture model.

Appendix B: Possible systematic uncertainties
As usual in any observation or experiment, our observational results may be subject to systematic offsets.Here, we describe possible underlying factors.According to our observations, there are two potential sources of systematic offsets: additive instrumental offsets in terms of the polarization degree and position angle.
An investigation conducted by Takahashi et al. (2018) demonstrated that the instrumental polarization degree offset (the so-called instrumental polarization degree) is very small for NIC and is effectively zero within the uncertainty level.The maximum offset occurs in the K s band, with values of (q off , u off ) = (−0.02± 0.30, −0.07 ± 0.31) %p, where the error bars represent the standard deviation of 180 sets of three unpolarized standard stars.For all the other cases, the offsets are ≤ 0.03 %p.Therefore, the instrumental polarization degree offset is expected to be at most of the same order as the random scatter (dP r ).Thus, we decided to skip this offset correction since it is negligible compared to random scatter.
The systematic offset in the instrumental position angle is estimated to be on the order of a few degrees (Takahashi 2019).The effect of this offset on P r is also expected to be very small.According to Eq. ( 1), an extreme systematic offset of ∆θ r = 10 • from θ r = 0 • would result in only ∆P r /P r = 0.06.Since the actual offset is ∆θ r ∼ 1 • and our data are |P r | < 2 %, we expect ∆P r ≪ 0.1 %p for all our data.
Another potential source of systematic offset comes from the polarization efficiency, p eff .As mentioned by Takahashi (2019), the NIC has an efficiency of p eff > 0.9, with values of approximately 100% within the error bars: (98, 95, 92) ± (6, 7, 12)% for the J, H, and K s bands, respectively.We applied these efficiency corrections to our data.The uncertainties provided are the standard deviations, not the standard errors.Unlike the previous two offsets discussed, p eff is generally considered to have a multiplicative effect (Akitaya et al. 2014), as it is caused mainly by the imperfect extinction ratio of the polarizer.Therefore, it linearly affects P r in Eq. ( 1).This means that P r might be systematically offset (multiplied) by ≲ 0.1 %p in our case (Table 2).We emphasize, however, that this does not affect our overall logic.
In the work by Masiero et al. (2022), they also found a hint of instrumental systematic offset with ∆P = +0.042and +0.088 %p in the J and H bands, respectively.In their work, they used an empirical law (Serkowski et al. 1975) with the assumption that the maximum polarization occurs at 0.55 µm for calibrating the instrumental offsets.Likely stemming from the different calibration processes, our data are slightly offset from theirs (e.g., 0.1 %p level near α min for Ceres).We again emphasize that this difference does not affect our discussion.

Fig. 1 .Fig. 2 .
Fig.1.PPCs of two objects in the NIR data from this study (star) and a previous publication (circle;Masiero et al. 2022).The solid lines are the least-square solutions, and the shaded areas show the MC 1σ bounds.
Fig. A.1.Effect of the choice of aperture size in the Ceres case.Different colors correspond to different bands (J, H, and K s ), while different markers correspond to different aperture modes.The solid lines are linear-exponential fits to each band and aperture mode.
Fig. A.2. Same as Fig. A.1 but for Vesta.

Table 1 .
Observing circumstances.Date: The UT date in YYYY-MM-DD.Time: The starting and ending UT time in HH:mm.N: The total number of obtained polarimetric sets (one set consists of four exposures).EXPTIME: The exposure times used.r h : The heliocentric distance to the target.r o : The observer-centric distance to the target.α: The phase angle.Skymonitor record: The NHAO skymonitor camera record archived on YouTube (also available in Appendix D).