11institutetext:  Korea Astronomy and Space Science Institute (KASI), 776
Daedeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea 22institutetext:
Department of Physics and Astronomy, Seoul National University, Gwanak-ro 1,
Gwanak-gu, Seoul 08826, Republic of Korea 33institutetext: SNU Astronomy
Research Center, Department of Physics and Astronomy, Seoul National
University, Gwanak-ro 1, Gwanak-gu, Seoul 08826, Republic of Korea
44institutetext: Center for Astronomy, University of Hyogo, 407-2 Nishigaichi,
Sayo, Hyogo 679-5313, Japan 55institutetext: Japan Spaceguard Association,
Bisei Spaceguard Center 1716-3 Okura, Bisei, Ibara, Okayama 714-1411, Japan
66institutetext: Nayoro Observatory, 157-1 Nisshin, Nayoro, Hokkaido 096-0066,
Japan 77institutetext: Department of Astronomy, Graduate School of Science,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
77email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Quantitative grain size estimation on airless bodies
from the negative polarization branch

II. Dawn mission targets (4) Vesta and (1) Ceres

Yoonsoo P. Bach 11 2 2 3 3 Masateru Ishiguro 22 3 3 Jun Takahashi 44 Jooyeon
Geem 22 3 3
Daisuke Kuroda 55 Hiroyuki Naito 66 Jungmi Kwon 77

(Received ; accepted )

###### Abstract

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. Most polarimetric studies on airless bodies
thus far have focused on optical wavelengths ($\lambda\lesssim 1\,\mathrm{\mu
m}$).

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\text{--}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.

###### Key Words.:

Minor planets, asteroids: general – Minor planets, asteroids: individual: (1)
Ceres – Minor planets, asteroids: individual: (4) Vesta

##  1 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(\lambda)$, and the size parameter,
$X\propto D/\lambda$, where $D$ represents grain size and $\lambda$ 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\lesssim 10\%$ and a shallowing trend for $A\gtrsim 10\%$. Second,
regarding $D/\lambda$, we found that the NPB stays relatively unchanged for
$D/\lambda\gg 10$. However, as $\lambda$ increases or $D/\lambda$ decreases,
the WD trend starts to appear at $D/\lambda\lesssim 5\text{--}10$. By focusing
on the effect of $D/\lambda$, 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/\lambda$. 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/\lambda$ 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\sim 950\,\mathrm{km}$), the presence of ultra-porous fine
particles ($D\ll 1\,\mathrm{\mu 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\sim 530\,\mathrm{km}$) is
predominantly covered with particles of sizes $\lesssim 20-25\,\mathrm{\mu 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\,\mathrm{\mu m}$ and $<10\,\mathrm{\mu 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/\lambda$ 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 ($\lambda\sim
1.6\,\mathrm{\mu m}$). We stress that our study is the first NIR polarimetry
study of a V-type asteroid (Vesta) and the first polarimetry study of airless
bodies at $\lambda>2\,\mathrm{\mu m}$ (the $\mathrm{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.

##  2 Observations and methods

###  2.1 Instruments and observations

Our observations were conducted using the Nishiharima Infrared Camera (NIC)
attached to the 2.0-m Nayuta telescope at the Nishi-Harima Astronomical
Observatory (NHAO; 134.3356∘E, 35.0253∘ N, 449 m elevation). The telescope has
an altazimuth mount with a f/12.0 focal ratio. The system consists of three
detector arrays, each $1024\times 1024$ pixel, with two dichroic mirrors such
that the three wavelengths are imaged simultaneously. The observed wavelengths
are J (central wavelength $\lambda=1.253\,\mathrm{\mu m}$ with a full width at
half maximum $\Delta\lambda=0.157\,\mathrm{\mu m}$), H
($\lambda=1.632\,\mathrm{\mu m}$, $\Delta\lambda=0.298\,\mathrm{\mu m}$), and
$\mathrm{K_{s}}$($\lambda=2.146\,\mathrm{\mu m}$,
$\Delta\lambda=0.313\,\mathrm{\mu m}$).

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\arcsec/\mathrm{pix}$. Each image has a size of $\sim 150\times
430\,\mathrm{pixels}$ ($24\arcsec\times 69\arcsec$) 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^{\circ}$, $45^{\circ}$, $22.5^{\circ}$, and
$67.5^{\circ}$, 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 $\alpha=4.1\degr$, 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).

Table 1: Observing circumstances.

Date | Time | Object | $N$ | EXPTIME | $r_{\mathrm{h}}$ | $r_{\mathrm{o}}$ | $\alpha$ | Skymonitor record
---|---|---|---|---|---|---|---|---
| UT |  |  | sec | au | au | ° |
2019-10-22 | 15:27-16:04 | Vesta | 21 | 2 | 2.5 | 1.6 | 10.4 | https://youtu.be/PaU3ign_sdY
2019-11-08(a)(a)footnotemark: $(a)$ | 12:24-17:27 | Vesta | 145 | 2, 3, 4, 5, 10, 30 | 2.5 | 1.6 | 4.1 | https://youtu.be/dg5ejH7fPd8
2019-11-21 | 12:22-16:37 | Vesta | 165 | 2 | 2.6 | 1.6 | 5.7 | https://youtu.be/_gHz8hH7eXQ
2019-12-18 | 13:20-14:57 | Vesta | 67 | 2, 3 | 2.6 | 1.8 | 15.6 | https://youtu.be/ceSwHR12S-s
2020-01-10 | 13:18-13:59 | Vesta | 31 | 2 | 2.6 | 2.0 | 20.7 | https://youtu.be/3Wkn83F4Qzg
2020-02-13 | 09:04-10:00 | Vesta | 30 | 2 | 2.6 | 2.5 | 22.5 | https://youtu.be/VzwG5y-18ME
2020-06-21 | 18:30-19:20 | Ceres | 18 | 20 | 3.0 | 2.5 | 18.9 | https://youtu.be/Y4H6QNyWcp0
2020-10-02(b)(b)footnotemark: $(b)$ | 09:19-10:30 | Ceres | 20 | 20 | 3.0 | 2.2 | 13.1 | https://youtu.be/muyaqzYcnb4
2020-10-18 | 08:50-09:59 | Ceres | 22 | 20 | 3.0 | 2.3 | 16.5 | https://youtu.be/6-SD0WrY1No
2021-09-09(c)(c)footnotemark: $(c)$ | 16:15-17:45 | Ceres | 25 | 20 | 2.8 | 2.5 | 20.8 | https://youtu.be/fUuu9HiQT8M
2021-11-07 | 12:34-22:40 | Ceres | 31 | 5 | 2.8 | 1.8 | 8.6 | https://youtu.be/W7xduRk5tSM

111 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_{\mathrm{h}}$: The
heliocentric distance to the target. $r_{\mathrm{o}}$: The observer-centric
distance to the target. $\alpha$: The phase angle. Skymonitor record: The NHAO
skymonitor camera record archived on YouTube (also available in Appendix D).
$a$$a$footnotetext: Unstable weather (refer to the skymonitor record).
$b$$b$footnotetext: Most H-band frames are saturated. $c$$c$footnotetext: All
H-band frames are saturated.

###  2.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 software222 https://github.com/ysBach/NICpolpy,
https://pypi.org/project/NICpolpy/  (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.

###  2.3 Polarization parameter calculation

The polarization degree of solar system objects is often described by the
properly rotated polarization degree (Lyot 1929):

$P_{\mathrm{r}}\coloneqq
P\cos(2\theta_{\mathrm{r}})\approx\frac{I_{\perp}-I_{\parallel}}{I_{\perp}+I_{\parallel}}\leavevmode\nobreak\
,$ (1)

where $P$ is the total linear polarization degree, $\theta_{\mathrm{r}}$ is
the position angle of the strongest electric vector relative to the normal
vector of the scattering plane (e.g., $0$ or $180^{\circ}$ indicates it is
normal to the scattering plane), and $I_{\perp}$ and $I_{\parallel}$ 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; $\alpha\coloneqq
180\degr-\mathrm{scattering\leavevmode\nobreak\ angle}$),
$P_{\mathrm{r}}(\alpha)$, is denoted the polarization phase curve (PPC). A PPC
at small $\alpha$ can be expressed with a three-parameter linear-exponential
function (Muinonen et al. 2002; Kaasalainen et al. 2003):

$P_{\mathrm{r}}(\alpha)=\theta_{1}\left(e^{-\alpha/\theta_{2}}-1\right)+\theta_{3}\alpha\leavevmode\nobreak\
,$ (2)

for free parameters $\theta_{1,\,2,\,3}$. The inversion angle $\alpha_{0}$ can
be solved numerically from $P_{\mathrm{r}}(\alpha_{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\coloneqq\left.\frac{dP_{\mathrm{r}}}{d\alpha}\right|_{\alpha=\alpha_{0}}=P_{\mathrm{r}}^{\prime}(\alpha=\alpha_{0})$,
inversion angle $\alpha_{0}$, $\alpha_{\mathrm{min}}$ and
$P_{\mathrm{min}}=P_{\mathrm{r}}(\alpha_{\mathrm{min}})$. By rearranging the
equation, we can generate the fitting parameters with more intuitive meanings
and facilitate subsequent statistical analyses of the parameters:

$P_{\mathrm{r}}(\alpha)=P_{\mathrm{r}}(\alpha;h,\,\alpha_{0},\,k)=h\frac{(1-e^{-\alpha_{0}/k})\alpha-(1-e^{-\alpha/k})\alpha_{0}}{1-(1+\alpha_{0}/k)e^{-\alpha_{0}/k}}\leavevmode\nobreak\
,$ (3)

where $k\,[\degr]$ is a scale parameter. Having rearranged the equations such
that two out of the three fitting parameters ($h$, $\alpha_{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 $\alpha_{\mathrm{min}}$ through the condition
$P_{\mathrm{r}}^{\prime}(\alpha=\alpha_{\mathrm{min}})=0$ as follows:

$\alpha_{\mathrm{min}}=-k\ln\left\\{\frac{k}{\alpha_{0}}\left(e^{-\alpha_{0}/k}-1\right)\right\\}\leavevmode\nobreak\
.$ (4)

Therefore, four important parameters, $h$, $\alpha_{0}$,
$\alpha_{\mathrm{min}}$, and
$P_{\mathrm{min}}=P_{\mathrm{r}}(\alpha_{\mathrm{min}})$, are obtained.

##  3 Results

###  3.1 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=\sqrt{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, $(\sum\sigma_{i}^{-2})^{-1/2}$, where
$\sigma_{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\theta_{P}=d\theta_{\mathrm{r}}$.
For $dP_{\mathrm{r}}$, we applied error propagation to Equation (1), while
neglecting the covariance between $dP$ and $d\theta_{\mathrm{r}}$. The maximum
value between this propagated uncertainty and $dP$ is obtained as
$dP_{\mathrm{r}}$.

Table 2: Polarimetric data reduction results.

Object | Filter | $n$ | $\alpha$ | $P$ | $\theta_{P}$ | $P_{\mathrm{r}}$ | $dP_{\mathrm{r}}$ | $\theta_{\mathrm{r}}$ | $d\theta_{\mathrm{r}}$
---|---|---|---|---|---|---|---|---|---
|  |  | $\degr$ | % | $\degr$ | % | % | $\degr$ | $\degr$
Ceres(a)(a)footnotemark: $(a)$ | J(a)(a)footnotemark: $(a)$ | $-$ | $8.3$ | $-$ | $-$ |  $-1.740$(a)(a)footnotemark: $(a)$ |  $0.020$(a)(a)footnotemark: $(a)$ |  $83.40$(a)(a)footnotemark: $(a)$ |  $0.10$(a)(a)footnotemark: $(a)$
Ceres(a)(a)footnotemark: $(a)$ | J(a)(a)footnotemark: $(a)$ | $-$ | $21.0$ | $-$ | $-$ |  $+0.840$(a)(a)footnotemark: $(a)$ |  $0.020$(a)(a)footnotemark: $(a)$ |  $22.30$(a)(a)footnotemark: $(a)$ |  $0.50$(a)(a)footnotemark: $(a)$
Ceres | J | 31 | $8.6$ | $1.620$ | $-2.37$ | $-1.619$ | $0.038$ | $88.91$ | $0.67$
Ceres | J | 17 | $13.1$ | $1.054$ | $+39.70$ | $-1.053$ | $0.023$ | $88.97$ | $0.61$
Ceres | J | 14 | $16.5$ | $0.388$ | $+35.13$ | $-0.383$ | $0.029$ | $85.44$ | $2.11$
Ceres | J | 17 | $18.9$ | $0.289$ | $-68.86$ | $+0.288$ | $0.032$ | $1.58$ | $3.15$
Ceres | J | 25 | $20.8$ | $0.898$ | $-83.42$ | $+0.898$ | $0.031$ | $1.01$ | $1.00$
Ceres(a)(a)footnotemark: $(a)$ | H(a)(a)footnotemark: $(a)$ | $-$ | $8.3$ | $-$ | $-$ |  $-1.610$(a)(a)footnotemark: $(a)$ |  $0.030$(a)(a)footnotemark: $(a)$ |  $86.60$(a)(a)footnotemark: $(a)$ |  $0.10$(a)(a)footnotemark: $(a)$
Ceres(a)(a)footnotemark: $(a)$ | H(a)(a)footnotemark: $(a)$ | $-$ | $21.0$ | $-$ | $-$ |  $+0.880$(a)(a)footnotemark: $(a)$ |  $0.020$(a)(a)footnotemark: $(a)$ |  $16.20$(a)(a)footnotemark: $(a)$ |  $0.20$(a)(a)footnotemark: $(a)$
Ceres | H | 31 | $8.6$ | $1.558$ | $-2.82$ | $-1.558$ | $0.046$ | $89.36$ | $0.85$
Ceres | H | 4 | $13.1$ | $1.155$ | $+40.00$ | $-1.153$ | $0.037$ | $88.68$ | $0.92$
Ceres | H | 15 | $16.5$ | $0.273$ | $+32.77$ | $-0.273$ | $0.028$ | $87.80$ | $2.97$
Ceres | H | 17 | $18.9$ | $0.356$ | $-69.82$ | $+0.356$ | $0.029$ | $0.61$ | $2.35$
Ceres | $\mathrm{K_{s}}$ | 31 | $8.6$ | $1.626$ | $-1.92$ | $-1.624$ | $0.039$ | $88.46$ | $0.68$
Ceres | $\mathrm{K_{s}}$ | 20 | $13.1$ | $1.063$ | $+40.65$ | $-1.061$ | $0.049$ | $88.03$ | $1.33$
Ceres | $\mathrm{K_{s}}$ | 14 | $16.5$ | $0.358$ | $+32.52$ | $-0.357$ | $0.046$ | $88.05$ | $3.72$
Ceres | $\mathrm{K_{s}}$ | 13 | $18.9$ | $0.535$ | $-70.83$ | $+0.535$ | $0.050$ | $0.40$ | $2.69$
Ceres | $\mathrm{K_{s}}$ | 24 | $20.8$ | $1.105$ | $-80.27$ | $+1.102$ | $0.041$ | $2.14$ | $1.06$
Vesta | J | 55 | $4.1$ | $0.657$ | $-49.51$ | $-0.657$ | $0.038$ | $89.53$ | $1.63$
Vesta | J | 158 | $5.7$ | $0.659$ | $+52.90$ | $-0.659$ | $0.028$ | $89.42$ | $1.22$
Vesta | J | 20 | $10.4$ | $0.674$ | $-7.78$ | $-0.672$ | $0.089$ | $87.92$ | $3.77$
Vesta | J | 67 | $15.6$ | $0.413$ | $+22.18$ | $-0.411$ | $0.033$ | $87.13$ | $2.28$
Vesta | J | 31 | $20.7$ | $0.128$ | $+11.06$ | $-0.122$ | $0.033$ | $81.19$ | $7.43$
Vesta | J | 30 | $22.5$ | $0.069$ | $-28.91$ | $-0.002$ | $0.047$ | $45.69$ | $19.78$
Vesta | H | 16 | $4.1$ | $0.636$ | $-52.04$ | $-0.632$ | $0.255$ | $87.05$ | $11.48$
Vesta | H | 99 | $5.7$ | $0.773$ | $+52.92$ | $-0.773$ | $0.037$ | $89.48$ | $1.36$
Vesta | H | 18 | $10.4$ | $0.792$ | $+1.05$ | $-0.770$ | $0.078$ | $83.25$ | $2.82$
Vesta | H | 67 | $15.6$ | $0.582$ | $+21.48$ | $-0.578$ | $0.033$ | $86.43$ | $1.60$
Vesta | H | 31 | $20.7$ | $0.197$ | $+15.94$ | $-0.195$ | $0.033$ | $86.08$ | $4.77$
Vesta | H | 30 | $22.5$ | $0.077$ | $+4.72$ | $-0.072$ | $0.041$ | $79.32$ | $15.30$
Vesta | $\mathrm{K_{s}}$ | 45 | $4.1$ | $0.906$ | $-49.01$ | $-0.906$ | $0.048$ | $89.97$ | $1.51$
Vesta | $\mathrm{K_{s}}$ | 158 | $5.7$ | $0.984$ | $+51.64$ | $-0.982$ | $0.047$ | $88.16$ | $1.37$
Vesta | $\mathrm{K_{s}}$ | 20 | $10.4$ | $1.017$ | $-1.04$ | $-1.003$ | $0.093$ | $85.34$ | $2.63$
Vesta | $\mathrm{K_{s}}$ | 66 | $15.6$ | $0.710$ | $+20.74$ | $-0.702$ | $0.037$ | $85.69$ | $1.51$
Vesta | $\mathrm{K_{s}}$ | 29 | $20.7$ | $0.366$ | $+8.43$ | $-0.337$ | $0.044$ | $78.56$ | $3.46$
Vesta | $\mathrm{K_{s}}$ | 29 | $22.5$ | $0.160$ | $-15.27$ | $-0.077$ | $0.055$ | $59.33$ | $9.81$
333 $n$: Number of used polarimetric datasets after data filtering and outlier
rejections (one set corresponds to four exposures). $\alpha$: Solar phase
angle. $P$ and $dP$: The total linear polarization degree and its uncertainty.
$\theta_{P}$: The position angle of the strongest electric field vector.
$P_{\mathrm{r}}$ and $dP_{\mathrm{r}}$: The proper linear polarization degree
referring to the scattering plan and its uncertainty. $\theta_{\mathrm{r}}$
and $d\theta_{\mathrm{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$$a$footnotetext: Data from
Masiero et al. (2022).

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 $\approx\sqrt{n}$. In their work, Masiero et al. (2023) accounted
for this discrepancy by quadratically adding $0.1\,\mathrm{\%p}$ scatter to
their reported $dP_{\mathrm{r}}$ values. Additionally, these uncertainties do
not include systematic errors (see Appendix B).

As shown in Table 2, we have confirmed that our $\theta_{\mathrm{r}}$ values
are generally close to integer multiples of $90\degr$, as expected from
symmetry considerations (Lyot 1929; Cellino et al. 2018). Consequently,
$\cos(2\theta_{\mathrm{r}})\approx 1$. The extreme exception is Vesta at
$\alpha=22.5\degr$, where its polarization degree is near zero; thus,
$\theta_{\mathrm{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 $\theta_{\mathrm{r}}$ values of
asteroids to measure the instrumental offset (assuming it must be an integer
multiple of $90\degr$) when $|P_{\mathrm{r}}|$ is large. For example, using
Ceres data at $\alpha=8.6\degr$, we find that the instrumental position angle
offset is $<2\degr$ 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 the figures, and the results of our
observations closely match the corresponding fitting curves. The shaded
regions in the figures represent the approximate 1$\sigma$ 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)444 We used
$P_{\mathrm{r}}=-P$ if neither $P_{\mathrm{r}}$ nor $\theta_{\mathrm{r}}$ was
available since such data are from $\alpha\ll\alpha_{0}$. We did not use R2
filter ($1/\lambda=1.20\,\mathrm{\mu 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/\lambda=1.05\,\mathrm{\mu
m}^{-1}$). Additionally, two data points in the $10\,\mathrm{\mu m}$ N band
(Johnson et al. 1983) were mistakenly labeled as the $0.333\,\mathrm{\mu m}$ N
band; therefore, we did not use these data.  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 $\alpha$s), after taking 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\,\mathrm{\%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 $\alpha_{0}$ affect the $h$ values (i.e., they have high leverage
values). Finally, we excluded data points with $>0.25\,\mathrm{\%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\text{--}0.5\,\mathrm{\%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.

Figure 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$\sigma$ bounds.
Figure 2: Same as Fig. 1 but with optical data (Lupishko 2022). The solid
lines are the least squares solutions for each wavelength range of the data,
and the cross markers (“x”) indicate the masked data points.

###  3.2 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 $\ll 1\sigma$. For Ceres, we
also included data from Masiero et al. (2022). More detailed information on
the MC simulation is provided in Appendix C.

Table 3: Fitting results. Object | Filter |  $h$ [$\%/\degr$] |  $\alpha_{0}$ [$\degr$] |  $k$ [$\degr$] |  $\alpha_{\mathrm{min}}$ [$\degr$] |  $P_{\mathrm{min}}$ [%]
---|---|---|---|---|---|---
Ceres | J | $0.2570^{+0.0043}_{-0.0042}$ | $17.80^{+0.06}_{-0.05}$ | $6.544^{+0.522}_{-0.435}$ | $6.99^{+0.14}_{-0.13}$ | $-1.754^{+0.025}_{-0.023}$
Ceres | H | $0.2567^{+0.0061}_{-0.0047}$ | $17.72^{+0.07}_{-0.06}$ | $7.728^{+0.967}_{-0.628}$ | $7.23^{+0.19}_{-0.14}$ | $-1.641^{+0.031}_{-0.027}$
Ceres | $\mathrm{K_{s}}$ | $0.3008^{+0.0126}_{-0.0064}$ | $17.46^{+0.13}_{-0.07}$ | $13.087^{+6.304}_{-1.567}$ | $7.77^{+0.35}_{-0.14}$ | $-1.636^{+0.057}_{-0.029}$
Vesta | J | $0.0588^{+0.0098}_{-0.0075}$ | $22.72^{+0.65}_{-0.54}$ | $4.124^{+1.124}_{-0.746}$ | $7.05^{+0.64}_{-0.55}$ | $-0.700^{+0.023}_{-0.023}$
Vesta | H | $0.0821^{+0.0140}_{-0.0071}$ | $23.26^{+0.41}_{-0.49}$ | $6.283^{+2.391}_{-0.855}$ | $8.38^{+0.70}_{-0.36}$ | $-0.831^{+0.030}_{-0.027}$
Vesta | $\mathrm{K_{s}}$ | $0.0905^{+0.0135}_{-0.0092}$ | $23.96^{+0.58}_{-0.63}$ | $5.067^{+1.175}_{-0.656}$ | $7.92^{+0.51}_{-0.37}$ | $-1.058^{+0.032}_{-0.030}$

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 $\mathrm{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 spectra555 http://smass.mit.edu/catalog.php , which covers
$\lambda\sim 0.45\text{--}2.5\,\mathrm{\mu 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 $\alpha=0\degr$, implying that the backscattering
amplitude is assumed to be identical over the wavelength range. The geometric
albedo at a wavelength of $\lambda=0.55\,\mathrm{\mu m}$ (FC2 filter of the
Dawn spacecraft) was adopted as $p_{\mathrm{0.55\,\mathrm{\mu m}}}=0.094\pm
0.007$ for Ceres (Ciarniello et al. 2017) and $0.38\pm 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\,\mathrm{\mu m}$ bezels (e.g., for
$0.45\text{--}0.60$, wavelengths from $0.45+0.01=0.46\,\mathrm{\mu m}$ to
$0.60-0.01=0.59\,\mathrm{\mu m}$). We used albedos only for
$\lambda>0.4\,\mathrm{\mu m}$.

Table 4: Geometric albedos and uncertainties.  Object | Wavelength($\mathrm{\mu m}$)/Filter | Albedo | Uncertainty
---|---|---|---
Ceres | 0.45–0.60 | 0.0931 | 0.0069
Ceres | 0.60–0.75 | 0.0957 | 0.0071
Ceres | 0.75–1.00 | 0.0943 | 0.0070
Ceres | J | 0.0897 | 0.0067
Ceres | H | 0.0921 | 0.0069
Ceres | $\mathrm{K_{s}}$ | 0.0942 | 0.0070
Vesta | 0.45–0.60 | 0.372 | 0.039
Vesta | 0.60–0.75 | 0.403 | 0.042
Vesta | 0.75–1.00 | 0.322 | 0.034
Vesta | J | 0.450 | 0.047
Vesta | H | 0.434 | 0.046
Vesta | $\mathrm{K_{s}}$ | 0.404 | 0.043

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_{\mathrm{min}}$–$h$–$\alpha_{0}$ space. Among these parameters,
the Umov laws ($h$–albedo and $P_{\mathrm{min}}$–albedo relations) are shown
in Figures 3–5 of Paper I. For optical wavelengths, rotational variations in
amplitude $<0.1\,\mathrm{\%p}$ and the Umov law were studied for Vesta
(Cellino et al. 2016). However, we did not consider the rotational variation.

In the $P_{\mathrm{min}}$–$\alpha_{0}$ space (Figure 7 in Paper I), the change
in the location of Vesta, likely starting in the $\mathrm{K_{s}}$ band, is
noteworthy. Vesta is near S-complex asteroids (Belskaya et al. 2017) but moves
toward the lunar fine region in the $\mathrm{K_{s}}$ band. For the Moon, most
regions exhibit a WD trend (because of the albedo effect and/or $D/\lambda$
effect). A similar WD trend is also observed for Vesta. In contrast, the NPB
shows almost no change in Ceres in the $P_{\mathrm{min}}$–$\alpha_{0}$ space.

##  4 Discussion

###  4.1 The Umov laws ($h$, $P_{\mathrm{min}}$, and albedo)

The $h$–albedo and $P_{\mathrm{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/\lambda$) 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_{\mathrm{min}}$) only in the $\mathrm{K_{s}}$ band, while the
albedo and $P_{\mathrm{min}}$ remain nearly constant.

###  4.2 $P_{\mathrm{min}}$–$\alpha_{0}$ space

The NPB of Vesta clearly shows a trend of WD increasing $\lambda$, as shown
in, for example, Figure 6 of Paper I. However, Ceres stays at a fixed location
over the investigated wavelengths ($<0.5\,\mathrm{\mu m}$ to $>2\,\mathrm{\mu
m}$). We therefore discuss the particle sizes on these objects based on these
observations.

####  4.2.1 Determining Vesta’s particle size

Vesta has albedo $\gg 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/\lambda$. We
argue that $D/\lambda\sim 5\text{--}10$ is reached in the H or
$\mathrm{K_{s}}$ band ($\lambda\sim 2\,\mathrm{\mu m}$), so Vesta is likely
covered with particles of size $D\gtrsim 10\text{--}20\,\mathrm{\mu m}$. Finer
particles ($\sim 1\,\mathrm{\mu 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, $\lesssim 20\,\mathrm{\mu m}$
particles were found to be necessary to reproduce the observed $\alpha_{0}$ of
Vesta. However, finer particles resulted in too high an albedo666 They wrote
“$0.27$ to $0.4$ in the green filter.”  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\,\mathrm{\mu m}$) covered with fine
($<10\,\mathrm{\mu m}$) particles (denoted Bereba B) can match both the
$\alpha_{0}$ and albedo at the same time.

However, based on the Dawn mission, Vesta’s geometric albedo is updated to
$p_{\mathrm{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
$\alpha=5\degr$ is $\gtrsim 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\,\mathrm{\mu m}$) of
HED powders best reproduced the Vesta $0.3\text{--}2.6\,\mathrm{\mu m}$
spectra. Hence, they also concluded that very fine ($<25\,\mathrm{\mu 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\,\mathrm{\mu 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\,\mathrm{\mu m}$.

Finally, utilizing the thermal conduction modeling algorithm (MacLennan &
Emery 2022b), we found two solutions for Vesta, $D\lesssim 10\,\mathrm{\mu m}$
and $D\gtrsim 50\,\mathrm{\mu m}$, both of which require a porosity $\gtrsim
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\sim 10\text{--}20\,\mathrm{\mu
m}$ particles.

####  4.2.2 The case of Ceres

Because Ceres shows no change in the $P_{\mathrm{min}}$–$\alpha_{0}$ space,
one could argue that $D/\lambda\gg 10$ until the $\mathrm{K_{s}}$ band ($D\gg
20\,\mathrm{\mu 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\lesssim 10\%$). Therefore, another possible
hypothesis is that the fact that $D/\lambda\ll 1\text{--}2$, even at the
shortest wavelength (i.e., $D\ll 0.5\text{--}1\,\mathrm{\mu 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 $\lesssim 1\,\mathrm{\mu 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 $\ll 1\,\mathrm{\mu m}$ and filament separation $\gg
1\,\mathrm{\mu m}$. These can support the small particle scenario.

Utilizing a thermal modeling algorithm similar to Vesta (MacLennan & Emery
2022b), we found two solutions, $D\lesssim 2\,\mathrm{\mu m}$ and $D\gtrsim
30\,\mathrm{\mu 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 modeling, 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 $\mathrm{K_{s}}$ band. This may indicate the existence of a fine
structure. However, as shown in Paper I, the deviation for low-albedo
($A\lesssim 10\%$) materials has been investigated less than that for higher-
albedo materials; thus, further investigations are needed.

###  4.3 Revisiting several previous works

Belskaya et al. (1987) demonstrated similar trends in $P_{\mathrm{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_{\mathrm{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 possible777 Belskaya & Shevchenko (2000) shows the geometric albedo in
the V band $A_{V}(\alpha=0\degr)\sim 0.21$ and the intensity ratio
$I(\alpha=0.3\degr)/I(\alpha=5\degr)=1.44$ for the 10 S-type samples; thus,
$A_{V}(\alpha=5\degr)<A_{V}(0\degr)\times I(5\degr)/I(0.3\degr)=15\%$. Hence,
they are at the boundary of low and high albedo regimes investigated in Paper
I. . 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_{\mathrm{r}}(\lambda;\alpha)|$ decreased with
increasing reflectance, which matches the Umov law, while (236) Honoria showed
the opposite trend over $\lambda$. 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.

We emphasize that if the ND trend is present on these Barbarian asteroids,
this is an expected behavior. Honoria and Luisa were observed at
$\alpha=7.1^{\circ}$ (where $P_{\mathrm{r}}$ gets deeper with increasing
$\lambda$) and $\alpha=26.9^{\circ}$ (slightly below $\alpha_{0}$, where
$P_{\mathrm{r}}$ gets shallower with increasing $\lambda$), 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/\lambda<1\text{--}2$). This example hints that the Umov law
should be applied with care; for example, one may use $h$ or
$P_{\mathrm{max}}$ ($P_{\mathrm{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 $\alpha_{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\,\mathrm{km}$ already possess fine
particles that can induce the ND trend, so $\alpha_{0}$ decreases as $D$
decreases (i.e., a decreasing $\alpha_{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.

##  5 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/\lambda$ 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. 1.

The Umov laws (slope–albedo and $P_{\mathrm{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/\lambda$.

  2. 2.

Vesta showed a characteristic change in the NPB when $\lambda$ was increased
to the H and $\mathrm{K_{s}}$ bands ($\lambda\gtrsim 2\,\mathrm{\mu m}$).
Combining our data with those of previous studies, we find that the optically
significant grain size on Vesta is $\sim 10\text{--}20\,\mathrm{\mu 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. 3.

For Ceres, there was no significant change in the NPB ($P_{\mathrm{min}}$,
$\alpha_{0}$). This indicates the potential existence of either very fine
particles ($D\ll 1\,\mathrm{\mu m}$) or large particles ($D\gg 20\,\mathrm{\mu
m}$).

  4. 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
$\alpha_{\mathrm{min}}$ can reveal valuable hints about the parameters that
change the NPB. A change in the polarimetric slope ($h$) of Ceres in the
$\mathrm{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.

###### Acknowledgements.

This work was supported by a National Research Foundation of Korea (NRF) grant
funded by the Korean government (MEST) (No. 2023R1A2C1006180). YPB would like
to thank Eric MacLennan for helpful cross-checking of the grain size
calculations and engaging discussions about Ceres. The authors want to thank
Joseph Masiero for providing constructive comments as a referee.

## References

  * Akitaya et al. (2014) Akitaya, H., Moritani, Y., Ui, T., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91474O
  * Bach et al. (2022) Bach, Y. P., Ishiguro, M., Takahashi, J., & Geem, J. 2022, Stars and Galaxies (arXiv:2212.14167), 5, 4
  * Bagnulo et al. (2015) Bagnulo, S., Cellino, A., & Sterzik, M. F. 2015, Monthly Notices of the Royal Astronomical Society, 446, L11
  * Bagnulo et al. (2009) Bagnulo, S., Landolfi, M., Landstreet, J. D., et al. 2009, Publications of the Astronomical Society of the Pacific, 121, 993
  * Belskaya et al. (2017) Belskaya, I. N., Fornasier, S., Tozzi, G. P., et al. 2017, Icarus, 284, 30
  * Belskaya et al. (2009) Belskaya, I. N., Levasseur-Regourd, A.-C., Cellino, A., et al. 2009, Icarus, 199, 97
  * Belskaya et al. (1987) Belskaya, I. N., Lupishko, D. F., & Shakhovskoi, N. M. 1987, Soviet Astronomy Letters, 13, 219
  * Belskaya & Shevchenko (2000) Belskaya, I. N. & Shevchenko, V. G. 2000, Icarus, 147, 94
  * Bendjoya et al. (2022) Bendjoya, P., Cellino, A., Rivet, J. P., et al. 2022, Astronomy and Astrophysics, 665, A66
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, Astronomy and Astrophysics, Supplement, 117, 393
  * Capria et al. (2014) Capria, M. T., Tosi, F., De Sanctis, M. C., et al. 2014, Geophysics Research Letters, 41, 1438
  * Cellino et al. (2016) Cellino, A., Ammannito, E., Magni, G., et al. 2016, Monthly Notices of the Royal Astronomical Society, 456, 248
  * Cellino et al. (2018) Cellino, A., Bagnulo, S., Belskaya, I. N., & Christou, A. A. 2018, Monthly Notices of the Royal Astronomical Society, 481, L49
  * Ciarniello et al. (2017) Ciarniello, M., De Sanctis, M. C., Ammannito, E., et al. 2017, Astronomy and Astrophysics, 598, A130
  * Delbo’ et al. (2007) Delbo’, M., dell’Oro, A., Harris, A. W., Mottola, S., & Mueller, M. 2007, Icarus, 190, 236
  * Delbo et al. (2015) Delbo, M., Mueller, M., Emery, J. P., Rozitis, B., & Capria, M. T. 2015, in Asteroids IV, 107–128
  * Doodson (1917) Doodson, A. T. 1917, Biometrika, 11, 425
  * Foreman-Mackey et al. (2019) Foreman-Mackey, D., Farr, W., Sinha, M., et al. 2019, The Journal of Open Source Software, 4, 1864
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Gil-Hutton (2017) Gil-Hutton, R. 2017, in ACM (Asteroid, Comets, Meteors) 2017, Montevideo, Poster1.d.17
  * Gundlach & Blum (2013) Gundlach, B. & Blum, J. 2013, Icarus, 223, 479
  * Hapke (2012) Hapke, B. 2012, Theory of reflectance and emittance spectroscopy (Cambridge university press)
  * Hasegawa et al. (2019) Hasegawa, S., Hiroi, T., Ohtsuka, K., et al. 2019, Publications of the Astronomical Society of Japan, 71, 103
  * Hiroi et al. (1994) Hiroi, T., Pieters, C. M., & Takeda, H. 1994, Meteoritics, 29, 394
  * Ishiguro et al. (2011) Ishiguro, M., Takahashi, J., Zenno, T., Tokimasa, N., & Kuroda, T. 2011, Annual Report of the Nishi-Harima Astronomical Observatory, 521, 13
  * Johnson et al. (1983) Johnson, P. E., Kemp, J. C., Lebofsky, M. J., & Rieke, G. H. 1983, Icarus, 56, 381
  * Kaasalainen et al. (2003) Kaasalainen, S., Piironen, J., Kaasalainen, M., et al. 2003, Icarus, 161, 34
  * Kwon et al. (2023) Kwon, Y. G., Bagnulo, S., & Cellino, A. 2023, Astronomy and Astrophysics, 677, A146
  * Lauretta et al. (2017) Lauretta, D. S., Balram-Knutson, S. S., Beshore, E., et al. 2017, Space Science Reviews, 212, 925
  * Le Bertre & Zellner (1980) Le Bertre, T. & Zellner, B. 1980, Icarus, 43, 172
  * Li et al. (2011) Li, J.-Y., Bodewits, D., Feaga, L. M., et al. 2011, Icarus, 216, 640
  * Li et al. (2013) Li, J.-Y., Le Corre, L., Schröder, S. E., et al. 2013, Icarus, 226, 1252
  * Li et al. (2019) Li, J.-Y., Schröder, S. E., Mottola, S., et al. 2019, Icarus, 322, 144
  * Lupishko (2022) Lupishko, D. 2022, NASA Planetary Data System, 1
  * Lyot (1929) Lyot, B. 1929, PhD thesis, Universite Pierre et Marie Curie (Paris VI), France
  * MacLennan & Emery (2022a) MacLennan, E. M. & Emery, J. P. 2022a, Planetary Science Journal, 3, 263
  * MacLennan & Emery (2022b) MacLennan, E. M. & Emery, J. P. 2022b, Planetary Science Journal, 3, 47
  * Martikainen et al. (2019) Martikainen, J., Penttilä, A., Gritsevich, M., Videen, G., & Muinonen, K. 2019, Monthly Notices of the Royal Astronomical Society, 483, 1952
  * Masiero et al. (2023) Masiero, J. R., Devogèle, M., Macias, I., Castaneda Jaimes, J., & Cellino, A. 2023, Planetary Science Journal, 4, 93
  * Masiero et al. (2022) Masiero, J. R., Tinyanont, S., & Millar-Blanchaer, M. A. 2022, Planetary Science Journal, 3, 90
  * Morrison (1977) Morrison, D. 1977, Icarus, 31, 185
  * Muinonen et al. (2002) Muinonen, K., Piironen, J., Kaasalainen, S., & Cellino, A. 2002, Mem. Societa Astronomica Italiana, 73, 716
  * Nakamura et al. (2023) Nakamura, T., Matsumoto, M., Amano, K., et al. 2023, Science, 379, abn8671
  * Park et al. (2008) Park, J., Liu, Y., Kihm, K. D., & Taylor, L. A. 2008, Journal of Aerospace Engineering, 21, 266
  * Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical recipes 3rd edition: The art of scientific computing (Cambridge university press)
  * Reddy et al. (2012) Reddy, V., Nathues, A., Le Corre, L., et al. 2012, Science, 336, 700
  * Rosner (1983) Rosner, B. 1983, Technometrics, 25, 165
  * Schröder et al. (2017) Schröder, S. E., Mottola, S., Carsenty, U., et al. 2017, Icarus, 288, 201
  * Schröder et al. (2021) Schröder, S. E., Poch, O., Ferrari, M., et al. 2021, Nature Communications, 12, 274
  * Serkowski et al. (1975) Serkowski, K., Mathewson, D. S., & Ford, V. L. 1975, Astrophysical Journal, 196, 261
  * Stetson (1987) Stetson, P. B. 1987, Publications of the Astronomical Society of the Pacific, 99, 191
  * Takahashi (2019) Takahashi, J. 2019, Stars and Galaxies, 2, 3
  * Takahashi et al. (2021) Takahashi, J., Itoh, Y., Matsuo, T., et al. 2021, Astronomy and Astrophysics, 653, A99
  * Takahashi et al. (2018) Takahashi, J., Zenno, T., Saito, T., & Itoh, Y. 2018, Stars and Galaxies, 1, 17
  * Zellner et al. (1974) Zellner, B., Gehrels, T., & Gradie, J. 1974, Astronomical Journal, 79, 1100
  * Zellner et al. (1977) Zellner, B., Leake, M., Lebertre, T., Duseaux, M., & Dollfus, A. 1977, Lunar and Planetary Science Conference Proceedings, 1, 1091

##  Appendix A Details of the polarimetry

A set of observations contains four exposures at HWP angles of
$\theta=0^{\circ}$, $45^{\circ}$, $22.5^{\circ}$, and $67.5^{\circ}$. In this
work, we call the first two and the last two the $q$-subset and $u$-subset, as
they are used to obtain Stokes $q\coloneqq Q/I$ and $u\coloneqq U/I$,
respectively. For each subset, four intensities were measured (the o- and
e-ray for each exposure). The measured intensities of the target are $a$ for
the o-ray of the first frame, $b$ for the e-ray of the first frame, $c$ for
the o-ray of the second frame, and $d$ for the e-ray of the second frame. The
Stokes parameter $s$ (either $q$ or $u$) is estimated as

$s=\frac{1-R}{1+R}\leavevmode\nobreak\ ,$ (5)

where

$R=\sqrt{ad/bc}\leavevmode\nobreak\ .$ (6)

To obtain $a$ to $d$, we used simple circular aperture photometry since the
point-spread function of the NIC is nearly circular. We always fixed the sky-
estimating region as a circular annulus with inner and outer radii of 70 and
90 pixels, respectively. This is determined so that it is large enough to
avoid any flux from our object even when the target is the brightest
throughout our whole observational period. For some selected cases, by
changing this annulus we obtain a result that changes less than the final
1$\sigma$ error bar for the $P_{\mathrm{r}}$ value. The sky value was
estimated based on the SExtrator-like algorithm (Bertin & Arnouts 1996).
First, 3-sigma clipping was performed for a maximum of 5 iterations to reject
extreme pixel values. Then, for the sigma-clipped mean, median, and standard
deviation ($s_{\mathrm{sky}}$), the modal value is
$m_{\mathrm{sky}}=(2.5\times\mathrm{median})-(1.5\times\mathrm{mean})$ if
$(\mathrm{mean}-\mathrm{median})/s_{\mathrm{sky}}\leq 0.3$; otherwise,
$m_{\mathrm{sky}}=\mathrm{median}$. This is a slightly modified version of the
widely used mode estimation from a unimodal distribution (Doodson 1917). Since
the size of the final image for each of the o- and e-rays is $130\times 410$
pixels (Bach et al. 2022), the regions outside these regions are regarded as
blank pixels in sky estimation.

In addition to $a$ to $d$, we also calculate their uncertainties by

$\Delta a=\sqrt{\frac{(1+\delta_{F})(a+m_{\mathrm{sky}})}{g}+\pi
r_{\mathrm{ap}}^{2}\left[\left(\frac{R_{e}}{g}\right)^{2}+s_{\mathrm{sky}}^{2}\right]}\leavevmode\nobreak\
,$ (7)

and similar for $b$ to $d$. Here, $\delta_{F}$ is a factor for coping with the
uncertainty in flat fielding, which we fixed to 0.02 (Stetson 1987, and
DAOPHOT4 uses the default 0.0075, Stetson, P. B. 2023, priv. comm.), $g$ is
the gain factor of the detector in electron/ADU (Bach et al. 2022),
$r_{\mathrm{ap}}$ is the aperture radius (Sect. A.2), and $R_{e}$ is the
readout noise in the electron unit. The S/N is calculated as
$\mathrm{S/N}_{a}=a/\Delta a$, and so on. The uncertainty in Stokes $s$ is

$\Delta s=\frac{R}{(R+1)^{2}}\sqrt{\left(\frac{\Delta
a}{a}\right)^{2}+\left(\frac{\Delta b}{b}\right)^{2}+\left(\frac{\Delta
c}{c}\right)^{2}+\left(\frac{\Delta d}{d}\right)^{2}}\leavevmode\nobreak\ .$
(8)

###  A.1 Effect of outliers

For each subset, we calculate four S/N values:
$\\{\mathrm{S/N}_{a},\,\mathrm{S/N}_{b},\,\mathrm{S/N}_{c},\,\mathrm{S/N}_{d}\\}$.
Then, the minimum of these is saved as the representative $\mathrm{S/N}$, and
the peak-to-peak difference (maximum minus minimum) is called $\mathrm{dS/N}$.
Next, we selected sets only if all the following criteria were met:

  1. 1.

$\mathrm{S/N}_{s}>\texttt{SNR\\_min}$ for both $s=q$ and $u$

  2. 2.

$\mathrm{S/N}_{s}/d\mathrm{S/N}_{s}>\texttt{SNR\\_factor}$ for both $s=q$ and
$u$

  3. 3.

$\Delta s<\texttt{ds\\_max}$ for both $s=q$ and $u$.

We selected $\texttt{SNR\\_min}=20$, $\texttt{SNR\\_factor}=5$, and
$\texttt{ds\\_max}=0.5\%$. When we change these, some data points are
affected, but the effects on $P_{\mathrm{min}}$ and $\alpha_{0}$ are marginal.
We introduce this screening method to remove the observations that were made
when thin clouds passed through the field of view (i.e., when the S/N was too
small or changed rapidly within two consecutive exposures). The choice of
these parameters generally does not change the $P_{\mathrm{r}}$ value at all
(much smaller than the 1$\sigma$ error bar). For the case of the H band of
Vesta on UT 2019-11-08, we found that different combinations of SNR_min and
SNR_factor (e.g., 30 and 4) can change it by the order of the
$0.1\text{--}0.2\,\mathrm{\%p}$ level, but this change is still within the
order of the error bar. Even if this data point at $\alpha=4.1\degr$ is very
uncertain, it is far from $\alpha_{\mathrm{min}}$ and $\alpha_{0}$, so the
estimated locations of Vesta in the $P_{\mathrm{min}}$–$\alpha_{0}$ space are
not affected.

After this first screening, we applied the generalized extreme Studentized
deviate test (Rosner 1983, implemented in scikit-
posthocs888https://github.com/maximtrp/scikit-posthocs, specifically
implemented in this pull request, https://github.com/maximtrp/scikit-
posthocs/pull/56) to remove an unknown number of outliers in the $q$-$u$
plane. This rejection was based on the distance of the data points with
respect to the average point in the $q$-$u$ plane, with a significance level
of 0.95. We did not find a strong deviation from the final results due to this
process.

###  A.2 Effect of aperture modes

Changing the centroiding algorithm can change the total flux by a small
amount, but this change is much smaller than the change in the sky modal
values. Thus, the remaining important factor that can affect the final result
is how we select the aperture radius. We tested multiple algorithms to
determine the best aperture size for estimating the incident flux.

Considering that the typical full width at half maximum of NIC is
approximately 15 pixels, we first calculated the aperture sum using a fixed
aperture radius of $r_{\mathrm{ap}}=25,\,30,\,35$ and $40$ pixels (denoted as
codes 25, 30, 35, and 40). Since the full widths at half maximum clearly vary
each night, this results in either a loss of flux (small radius) or a larger
scatter (large radius). What makes this more challenging is that we sometimes
defocused the image intentionally to avoid saturation. Thus, a single-aperture
size is sometimes not appropriate, even over a single night.

Another possible method is to change the aperture size for each subset. Since
each exposure results in two images of the target (the o-ray and the e-ray),
there are four images. For each of these images, we obtained the modal pixel
values as a function of the radial distance from the centroid of our target
(for distances of 5 to 60 pixels). The modal values were calculated using an
algorithm identical to that used for sky modal value estimation by using a
circular annulus of 1 pixel thickness. Then, we selected the maximum radius
that has a modal value $>m_{\mathrm{sky}}+1\times s_{\mathrm{sky}}$ (denote it
$r_{\mathrm{ap0}}$). If the object is too faint and no such radius exists, we
use $r_{\mathrm{ap0}}=5$. For the four $r_{\mathrm{ap0}}$ values obtained in
this way, we can choose a fixed $r_{\mathrm{ap}}$ for all four images using
either the minimum or maximum of the $r_{\mathrm{ap0}}$’s (denoted as code min
and max, respectively). Finally, the point spread function of the o- and
e-rays is expected to be different for NIC since it uses a beam displacer.
Thus, we can also try using separate $r_{\mathrm{ap}}$ for the o-ray and the
e-ray, $r_{\mathrm{ap,o}}$ and $r_{\mathrm{ap,e}}$, respectively.
$r_{\mathrm{ap,o|e}}$ can be the smallest among the two $r_{\mathrm{ap0,o|e}}$
values (denoted as code minoe) or the largest (maxoe).

For example, if we assume the $r_{\mathrm{ap0}}$ values are 15 and 16 for the
first frame and 17 and 18 for the second frame (o- and e-ray, respectively),
under our notation $r_{\mathrm{ap}}$ would be $25$ for the 25 code, regardless
of $r_{\mathrm{ap0}}$ values; $r_{\mathrm{ap}}=15$ for min;
$r_{\mathrm{ap}}=18$ for max; $r_{\mathrm{ap,o|e}}=15|16$ for minoe; and
$r_{\mathrm{ap,o|e}}=17|18$ for maxoe.

We find our final results to be extremely robust to all of these different
aperture models (Figs. 3 and 4). For Vesta, (i) the inversion angle
$\alpha_{0}$ might be slightly uncertain, especially because there are no data
at $\alpha>\alpha_{0}$, and (ii) the error bar of the smallest $\alpha$
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.

Figure 3: Effect of the choice of aperture size in the Ceres case. Different
colors correspond to different bands (J, H, and $\mathrm{K_{s}}$), while
different markers correspond to different aperture modes. The solid lines are
linear-exponential fits to each band and aperture mode. Figure 4: Same as Fig.
3 but for Vesta.

##  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 $\mathrm{K_{s}}$ band,
with values of $(q_{\mathrm{off}},u_{\mathrm{off}})=(-0.02\pm 0.30,-0.07\pm
0.31)\,\mathrm{\%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 $\leq 0.03\,\mathrm{\%p}$. Therefore, the instrumental
polarization degree offset is expected to be at most of the same order as the
random scatter ($dP_{\mathrm{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_{\mathrm{r}}$ is also expected to be very small. According to Eq. (1), an
extreme systematic offset of $\Delta\theta_{\mathrm{r}}=10\degr$ from
$\theta_{\mathrm{r}}=0\degr$ would result in only $\Delta
P_{\mathrm{r}}/P_{\mathrm{r}}=0.06$. Since the actual offset is
$\Delta\theta_{\mathrm{r}}\sim 1\degr$ and our data are
$|P_{\mathrm{r}}|<2\,\%$, we expect $\Delta P_{\mathrm{r}}\ll
0.1\,\mathrm{\%p}$ for all our data.

Another potential source of systematic offset comes from the polarization
efficiency, $p_{\mathrm{eff}}$. As mentioned by Takahashi (2019), the NIC has
an efficiency of $p_{\mathrm{eff}}>0.9$, with values of approximately $100\%$
within the error bars: $(98,\,95,\,92)\pm(6,7,12)\%$ for the J, H, and
$\mathrm{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_{\mathrm{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_{\mathrm{r}}$ in Eq. (1). This means that
$P_{\mathrm{r}}$ might be systematically offset (multiplied) by $\lesssim
0.1\,\mathrm{\%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 $\Delta P=+0.042$ and $+0.088\,\mathrm{\%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\,\mathrm{\mu m}$ for calibrating the instrumental offsets.
Likely stemming from the different calibration processes, our data are
slightly offset from theirs (e.g., $0.1\,\mathrm{\%p}$ level near
$\alpha_{\mathrm{min}}$ for Ceres). We again emphasize that this difference
does not affect our discussion.

##  Appendix C Monte Carlo simulation and uncertainties

We adopted the Markov chain MC package emcee Foreman-Mackey et al. (2013,
2019) with an affine-invariant ensemble sampler for our fitting process. We
first calculated the least-square solution for each object and each filter.
The priors for $h$, $\alpha_{0}$, and $k$ in Eq. (3) are chosen to be uniform
in the range of $h\in[0.01,10]$ (in %/°), $\alpha_{0}\in[10,30]$, and
$k\in[10^{-5},100]$ (both in degrees). The initial positions for 32 chains are
randomly selected near this position in the three-dimensional parameter space.
Then, with the stretch-move algorithm, 3000 samples per chain are drawn. The
first 1000 samples from all the samples were discarded; thus, 64000 MC samples
were left per object per filter. The $\alpha_{\mathrm{min}}$ and corresponding
$P_{\mathrm{min}}$ values are calculated for each sample after the sampling
process (Eq. 4). We assumed that the uncertainty of each observation was
purely Gaussian, so the log-likelihood method used the usual figure of merit:

$\mathcal{L}=-\frac{1}{2}\sum_{i=1}^{N}\left(\frac{P_{\mathrm{r,obs,i}}-P_{\mathrm{r}}(\alpha;h,\alpha_{0},k)}{dP_{\mathrm{r,obs,i}}}\right)^{2}\leavevmode\nobreak\
.$ (9)

Here, $P_{\mathrm{r,obs,i}}$ and $dP_{\mathrm{r,obs,i}}$ are the observed
polarization degree of the $i$th data point and its uncertainty, respectively.
The log posterior probability function is calculated by adding the log of the
prior function to $\mathcal{L}$.

Then, we discarded samples if they were outside any of the $\pm 5\sigma$
ranges of the three parameters (the asymmetry $\sigma$ for each parameter was
calculated from the 16th, 50th, and 84th percentiles of the MC samples). This
process discarded up to a few hundred samples. We checked that the randomness
of the MC sampling did not change the results by more than $0.1\sigma$.

To obtain the 1$\sigma$ uncertainty range of the curves in Figs. 1 and 2, we
first calculated the log-likelihood of each sample and converted it to the
usual $\chi^{2}$ statistic ($\chi^{2}=-2\mathcal{L}$). Then, the samples with

$\chi^{2}<\mathrm{min}\\{\chi^{2}\\}+\chi^{2}_{\mathrm{ppf}}(0.6827,\nu=3)$
(10)

are selected. Here, $\mathrm{min}\\{\chi^{2}\\}$ is the lowest $\chi^{2}$
value of all the samples (which should be nearly identical to that of the
least-square solution), and $\chi^{2}_{\mathrm{ppf}}$ is the percent point
function, that is to say, the inverse of the cumulative distribution function
of the $\chi^{2}$ distribution of the degrees of freedom $\nu$. The constant
0.6827 is selected specifically for the “1$\sigma$” range. In our case, Eq.
(3) contains three parameters of interest, so $\nu=3$. Model parameters within
the 1$\sigma$ uncertainty range were selected in the one-tailed test, which is
the method suggested in Press et al. (2007, Sect.
15.6)999http://numerical.recipes/book/book.html.

##  Appendix D Data and code availability

We obtained approximately 8 GB and 13 GB of raw data for Ceres and Vesta,
respectively. There are approximately 16 GB of dome flats, and some standard
star observations are approximately 16 GB. Then, the preprocessed data (levels
1-4, Bach et al. 2022) are approximately 3 times the total size (100+ GB of
data). Due to the large size of the total dataset, we only uploaded (i) the
raw (level 0) data of Ceres and Vesta, (ii) the final (level 4) data of Ceres
and Vesta, (iii) the corresponding log files produced from NICpolpy (Bach et
al. 2022), (iv) the used calibration frames, (v) the automatically generated
observational logs in PDF format and (vi) the skymonitor record videos. They
are available via Zenodo101010https://doi.org/10.5281/zenodo.7788352. Part of
the codes and data used for this publication are available via
GitHub111111https://github.com/ysBach/BachYP_etal_CeresVesta_NHAO