© ESO 2021

1 Introduction

Knowledge of the surface properties of asteroids provides insights into the nature of the surface materials and structures, which in turn are fundamental for different reasons. These properties are, for instance, important in modeling the formation of regolith-like materials (Delbó et al. 2014), the weakening and degradation of boulders on asteroids through thermal shocks (Molaro et al. 2017), and space-weather processes (Brunetto et al. 2015). Moreover, surface properties are also crucial for planning spacecraft interactions with the surface of an asteroid, including deflection missions (Bruck Syal et al. 2016).

The surface thermal inertia informs us about porosity and cohesion. Low values are indicative that a thermal insulting layer coats the surface (see Alí-Lagoa et al. 2020, and references therein), consistent with a dusty, porous regolith, while high values are specific for rock-like materials. Thermal inertia depends on density ρ, heat capacity C, and thermal conductivity K. However, while density and heat capacity typically vary within a factor of a few, the uncertainty on thermal conductivity spans a range of more than four orders of magnitude (Delbó et al. 2015). Thermal conductivity therefore is the most important parameter to constrain.

Before in-situ explorations of asteroids, it was thought that only very large objects could be covered by regolith (see Murdoch et al. 2015, and references therein). However, the visit of the Galileo and NEAR-Shoemaker NASA missions to asteroids (243) Ida, (433) Eros, and (951) Gaspra revealed that even smaller asteroids can retain regolith-like materials on their surfaces. Furthermore, with the visit of the JAXA Hayabusa mission to asteroid (25143) Itokawa in 2005, it became evident from its in-situ observations and measurements that despite the low gravitational environment, even subkilometer bodies can retain gravel particles, boulders, or regolith on their surfaces (Miyamoto et al. 2007). More recently, data collected in situ on asteroids (101955) Bennu and (162173) Ryugu by the OSIRIS-REx and the Hayabusa 2 missions, respectively, did not reveal fine dust on their surfaces, but still showed grain particles and boulders of different sizes (Lauretta et al. 2019a; Michikami et al. 2019; Morota et al. 2020; Sugita et al. 2020; Susorney et al. 2020).

These findings changed our perspective on the ability of kilometer- and subkilometer-sized asteroids to preserve regolith-like materials on their surfaces. This new perspective raises the question whether even smaller objects, from a few dozen to a few hundred meters, can also be covered by a dust layer. The extremely low-gravity environment on such bodies would suggest that a dust layer isnot likely. When we further restrict our considerations to small superfast rotators, this possibility seems even less probable.

Pravec & Harris (2000) found that objects larger than about 150 m very rarely have rotation periods shorter than about 2.2 h. This barrier can be explained by the fact that for diameters between about 150 m and 10 km, every Solar System body has a rubble-pile internal structure (Walsh 2018) that is not gravitationally strong enough to counteract the centripetal force due to the fast rotation, thus causing a breakup of the object. On the other hand, there is no such limit for the spin period of smaller objects, so that they are allowed (and often seen) to rotate very fast, even with periods of a few minutes. Asteroids smaller than 150 m in diameter are therefore considered to be mostly monolithic objects (see, e.g., Pravec & Harris 2000; Whiteley et al. 2002; Polishook 2013), with strong internal forces that hold the body intact even at very fast spin rates. Figure 1 shows the measurements of the spin period of near-Earth objects (NEOs) as a function of the diameter of the body, taken from the Asteroid Lightcurve Database¹ (LCDB) developed by Warner et al. (2009). We can clearly see the spin barrier for diameters larger than 150 m, while at a few dozen meters, the spin periods can be as short as few to a few dozen minutes (possibly even shorter, but typical observation scenarios would not detect them). Such fast spin rates, coupled with the very low gravitational environment due to the small size, are expected to cause the ejection of material on the surface, thus preventing regolith grains from being retained.

Determining values of thermal parameters for asteroids is a challenge because direct measurements from Earth cannot be performed. Consequently, reliable thermal inertia (conductivity) estimates are currently available for only a relatively small number of asteroids (see, e.g., Delbó et al. 2015; Harris & Drube 2016; Marciniak et al. 2019). Most often, thermal inertia is derived from an analysis of observations in the infrared band. Unless the observing facility is space borne, such data for very small objects are rare.

The dynamics of bodies up to about 30 km in diameter is known to be subject to thermal effects that are caused by solar radiation. Objects are heated by the Sun, so that they reradiate away the energy in the thermal waveband, creating a small thrust and causing a drift in semimajor axis that is perceived over decades: this physical phenomenon is known as the Yarkovsky effect (Rubincam 1995, 1998; Farinella et al. 1998; Bottke et al. 2006; Vokrouhlický et al. 2015). The Yarkovsky effect perturbs the trajectory, but it is generated by a thermal force that depends on material properties and the internal structure of asteroids. Therefore the Yarkovsky effect is a phenomenon in which orbital dynamics is linked to asteroid composition and physical properties.

The detection of the Yarkovsky effect has become possible as part of orbit determination (i.e., the osculating orbital elements of epoch) ona regular basis when high-precision astrometry observations are available. So far, this applies uniquely to NEOs. An alternative approach to estimating the drift induced by the Yarkovsky effect in the orbital semimajor axis has been formulated for the first time by Nesvorný & Bottke (2004), who used backward orbit propagation of the Karin cluster members and monitored convergence of the secular angles. The same method has been later applied to several other young asteroid families. However, this approach provides less accurate results than the direct orbit determination mentioned above. Additionally, data for small and fast-rotating asteroids, relevant in context of our paper, are not available for the main-belt families.

Depending on other available information, measurements of the Yarkovsky effect allow constraining the physical properties of small objects. Chesley et al. (2003) were able to deduce an interval of values for the thermal conductivity of (6489) Golevka, compatible with the existence of thin regolith. Additional examples include Farnocchia et al. (2013), who were able to constrain the values of the thermal inertia of six NEOs for which the obliquity was known. Similarly, using thermal infrared observations, we can estimate the thermal inertia of an asteroid using a suitable thermophysical model (see Delbó et al. 2015, and references therein). Combining these measurements with the Yarkovsky drift estimation, we can break the degeneracy between thermal inertia and density, and can directly derive the bulk density of the object (Chesley et al. 2014; Farnocchia et al. 2014; Mommert et al. 2014a,b; Rozitis et al. 2014; Reddy et al. 2016).

In this paper we develop a statistical method for constraining the thermal conductivity of NEOs. Our approach is based on the comparison of the measured and theoretically predicted Yarkovsky drift. Using suitable distributions for the parameters that determine the magnitude of the drift, we perform a Monte Carlo simulation, and produce a probability density function for the thermal conductivity K. This probability informs us about the most likely values of K, and it can be used to exclude some other values. We applied our model to asteroid (499998) 2011 PT, a superfast rotating NEO with a diameter of about 35 m. Interestingly, we found that the measured semimajor axis drift is large enough such that it can be achieved only for low values of thermal conductivity. Two most likely values are about 0.0001 and 0.005 W m⁻¹ K⁻¹, suggesting that a thermal insulating layer coats the surface of 2011 PT.

Fig. 1 Measured spin periods for NEOs as of September 2020 and taken from LCDB (see Warner et al. 2009).

2 Methods

The semimajor axis drift da∕dt caused by the Yarkovsky effect for a spherical body placed on a circular orbit around the Sun can be expressed through analytical formulas, and it is given by the sum of the seasonal and diurnal effects (see Appendix A). It depends on the orbital parameters and on the physical and thermal characteristics of the body, that is, the semimajor axis a, the diameter D, the density ρ, the thermal conductivity K, the heat capacity C, the obliquity γ, the rotation period P, the absorption coefficient α, and the emissivity ε. Of these parameters, the thermal conductivity K is the most uncertain parameter because it strongly depends on the composition of the surface and on the porosity, and it can vary by several orders of magnitude. If a measurement ${(\text{d}a/\text{d}t)}_{\text{m}}$ of the Yarkovsky effect is available from astrometry (see, e.g., Farnocchia et al. 2013; Del Vigna et al. 2018; Greenberg et al. 2020), solving the relation of model versus observed drift-rate $$\left(\frac{\text{d}a}{\text{d}t}\right)(a,D,\rho ,K,C,\gamma ,P,\alpha ,\epsilon )={(\frac{\text{d}a}{\text{d}t})}_{\text{m}},$$(1)

for K allows us to estimate the thermal conductivity, which can be examined statistically using a Monte Carlo method. To this purpose, in Eq. (1) the measured parameters are either fixed to a specific value when their uncertainty is negligible, or they can be modeled assuming a Gaussian error distribution. The parameters for which measurements are not available or that are subject to large uncertainties can be modeled using a population-based probability density function (PDF).

2.1 Probabilistic model of the parameters of 2011 PT

Asteroid (499998) 2011 PT is an NEO discovered with the Pan-STARRS 1 survey telescope on August 1, 2011. With a semimajor axis of a = 1.3123 au and an eccentricity of about e = 0.2147, it belongs to the Amor group, meaning that its perihelion distance q = a(1 − e) satisfies 1.017au < q < 1.3au. Hereafter we describe the probabilistic model of the orbital and physical parameters of 2011 PT that we used to compute the drift induced on the semimajor axis from the theoretical model of the Yarkovsky effect, that is, the left-hand side of Eq. (1).

2.1.1 Determining the absolute magnitude

The absolute magnitude H does not play a role in the analytical model of the Yarkovsky effect of Eq. (1) (see also Appendix A), but it is related to the size of the object and it is needed to model the albedo distribution (see Sect. 2.1.2). The JPL Small-Body Database² (SBDB) reports a value of H = 23.9 mag, but it is known that orbital catalogs all suffer uncertainties in H, including systematic effects (Pravec et al. 2012). To overcome this issue and improve the accuracy of the results, we remeasured all available Pan-STARRS images and recomputed an absolute magnitude of H = 24.07 ± 0.42 mag using the phase angle computed from the orbital geometry at each observation and the corrections to V band as given by Denneau et al. (2013). We note that the available data do not allow a more accurate determination of H.

2.1.2 Albedo distribution

Although the albedo does not enter in the formulation of the Yarkovsky effect as given by Eq. (1) directly, it is related to both the density and the diameter of the object. The value p_V of the albedo gives a first indication of the taxonomic type, which determines the density, and it provides an estimate of the diameter by the conversion formula (see, e.g., Bowell et al. 1989; Pravec & Harris 2007), $$D=\frac{1329\text{km}}{\sqrt{p_V}}{10}^{-H/5},$$(2)

where H is the absolutemagnitude. A PDF of the albedo of a specific NEO can be constructed using the NEO population model by Granvik et al. (2018) and the NEO albedo distribution by Morbidelli et al. (2020).

Given the orbital elements (a, e, i) and the absolute magnitude H, the population model by Granvik et al. (2018) provides the probability for an NEO to come from each main-belt source region, which are the ν₆ secular resonance, the 3:1, 5:2, and 2:1 Jupiter mean-motion resonances, the Hungaria region, the Phocaea region, and the Jupiter-family comets (JFC). For the sake of notation, we denote this probability with P_s (a, e, i, H), where the subscript s is used to denote a specific source route. Numerical values corresponding to (499998) 2011 PT are given in Table 1.

In the NEO albedo distribution by Morbidelli et al. (2020), objects are binned into three albedo categories: category (1) p_V ≤0.1; category (2) 0.1 <p_V ≤ 0.3; category (3) 0.3 < p_V. For each source region, the model provides the fraction of delivered bodies that belong to each albedo category c₁, c₂, and c₃. For the sake of notation, we denote the fractions with p_s(c_i), i = 1, 2, 3, where s again indicates the source region. The exact numerical values are reported in Morbidelli et al. (2020). The Hungaria region mostly provides high-albedo objects, in agreement with the fact that this part of the main belt contains a large number of E-type asteroids (DeMeo & Carry 2013, 2014). The inner belt regions, that is, the ν₆ secular resonance and the 3:1 Jupiter mean motion resonance, mostly provide asteroids in the S-complex. The only exception is the Phocaearegion, which supplies an approximately equal number of dark (category 1) and bright asteroids (categories 2 and 3). The larger fraction of dark asteroids coming from the Phocaea region compared to the rest of the inner main belt is probably due to a dark asteroid family that is located therein (Novaković et al. 2017). Nevertheless, the probability for 2011 PT to originate in this region is very small. The middle and outer belt regions, that is, the 5:2 and 2:1 Jupiter mean motion resonances, account for a higher percentage of C-complex asteroids, in agreement with the fact that primordial carbonaceous bodies are found to be more abundant at higher semimajor axis values (DeMeo & Carry 2013, 2014).

For each escape route, we defined an albedo probability function p_s (p_V) assigning a uniform density in the first and second categories, and an exponentially decaying density in the third category, as done in Morbidelli et al. (2020), $$p_s(p_V)=\{\begin{array}{ll}\frac{p_s(c_1)}{0.1},\hfill & p_V\le 0.1,\hfill \\ [2ex]\frac{p_s(c_2)}{0.2},\hfill & 0.1<p_V\le 0.3,\hfill \\ [2ex]p_s(c_3)\frac{{2.6}^{-\frac{p_V-0.3}{0.1}}}{{\displaystyle \underset{0.3}{\overset{1}{\int }}2}{.6}^{-\frac{x-0.3}{0.1}}\text{d}x},\hfill & p_V>\mathrm{0.3.}\hfill \end{array}$$(3)

Given the probabilities P_s(a, e, i, H) for the object to come from each source region, the albedo PDF p(p_V) for a specific NEO is defined by $$p(p_V)={\sum }_{s=1}^7{P}_s(a,e,i,H)p_s(p_V).$$(4)

The resulting albedo PDF for 2011 PT is shown in the top left panel of Fig. 2.

Fig. 2 Probability density function for the albedo pV (top left panel), the density ρ (top right panel), the diameter D (bottom left panel), and the obliquity γ (bottom right panel) computed for asteroid (499998) 2011 PT.

2.1.3 Density distribution

In order to obtain the PDF of the density of 2011 PT, we first divided asteroids into three main complexes: (1) the C-complex, composed ofdark carbonaceous objects, which includes C-, P-, D-, B-, and T-types, (2) the S-complex, composed of stony objects, which gathers S-, Q-, V-, A-, and R-types, and (3) the X-complex, composed of objects with different compositions, namely the intermediate-albedo M-type, and the high-albedo E-type.

For the S-complex we take as a reference the density of the S-type because they are the most common and numerous type (DeMeo & Carry 2014; DeMeo et al. 2015). The density is reported to be ρ = 2720 ± 540kg m⁻³ (Carry 2012).

On the other hand, C-type asteroids are the most common in the C-complex, therefore we assumedtheir density to model this population. Carry (2012) reported a density of ρ = 1330 ± 580kg m⁻³ for the C-type, which seemsto be slightly high when compared to the precise measurements performed on asteroids (101955) Bennu and (162173) Ryugu, however. (101955) Bennu is a B-type asteroid with a composition similar to the aqueously altered CM carbonaceous chondrite, and its density is 1190 ± 13 kg m⁻³ (Scheeres et al. 2019), while (162173) Ryugu is a Cb-type asteroid with a composition similar to carbonaceous chondrites, and a density of 1190 ± 20 kg m⁻³ (Watanabe et al. 2019). These two precise measurements suggest that the value of the density of C-type asteroids proposed in Carry (2012) may be too high, and this deviation might be due to large errors introduced by the method used to produce the estimates. The difference may also arise because asteroids tend to have higher density with increasing diameter (see, e.g., Carry 2012; Milani et al. 2014). For this reason, we assumed a lower value for the C-complex, setting the density to ρ = 1200 ± 300kg m⁻³.

As previously said, the X-complex shows a different variety of densities and compositions, but the E-type, which is connected to Aubrite meteorites, is the most common type delivered by the Hungaria region (Ćuk et al. 2014; Binzel et al. 2019), which is one of the main NEO source regions (Granvik et al. 2018). We therefore assumed the density of E-type asteroids as a reference for the X-complex. Carry (2012) reported only one measurement with an accuracy better than 20% for this asteroid type, and it is equal to 2600 ± 200 kg m⁻³. However, drawing a statistical picture from one single measurement may not be correct: for this reason, we established other values. To this purpose, we took into account the grain density of the connected Aubrite meteorites, which is reported to be ρ = 3200 ± 80kg m⁻³ (Ostrowski & Bryson 2019), and we corrected it taking the porosity into account. The porosity distribution was deduced from the values reported in Carry (2012): we took all the asteroids with porosities higher than zero and relative accuracy better than 50% in the density measurements. Then we fit these data with a gamma distribution, see Fig. 3.

Assuming a Gaussian distribution for the density of Aubrite, we corrected it for the porosity using the formula $$\rho ={\rho }_{\text{m}}(1-\frac{P(\%)}{100}),$$(5)

from Carry (2012), where P(%) is the porosity and ρ_m is the density of the meteorite. In this manner, we found a density of ρ = 2350 ± 520 kg m⁻³, which we used for the purpose of this work.

The three taxonomic complexes could be also linked to the three albedo categories defined by Eq. (3). In this respect, we assumed that category 1, which refers to dark asteroids, corresponds to the C-complex, category 2 corresponds to the S-complex because it refers to moderate albedo values, and category 3 is associated with the X-complex that contains higher albedo objects. This can be used to construct the PDF for the density of 2011 PT.

To generate a distribution of the density, we first generated a random value of the albedo according to the distribution defined by Eq. (4). Based on the generated albedo, we then determined the albedo category and associated asteroid complex. Finally, we randomly generated a value of the density according to the density distribution of the corresponding asteroid complex, which is assumed to be log-normal (see, e.g., Spoto et al. 2014; Tardioli et al. 2017). The PDF for ρ obtained with this procedure, smoothed using the kernel density estimation, is reported in the top right panel of Fig. 2. As expected, for 2011 PT, higher density solutions (> 2000 kg m⁻³) are much more likely than low-density solutions (<1500 kg m⁻³).

Fig. 3 Fit of the asteroid porosity with a gamma distribution using data from Carry (2012).

2.1.4 Size distribution

The size of an object is estimated using its absolute magnitude H and its albedo p_V by the conversion formula of Eq. (2). Because we determined an absolute magnitude of H = 24.07 ± 0.42 mag, the diameter is mostly affected by the uncertainty on the albedo, and it can vary by even an order of magnitude when the range of p_V is wide.

The distribution of the diameter D is obtained again by conversion of the albedo distribution of Eq. (4): random values of p_V are generated first, and then converted into diameter using Eq. (2). The uncertainty in the determination of the absolute magnitude H is taken into account assuming a Gaussian distribution of the errors. The PDF for the diameter of 2011 PT obtained with this procedure, smoothed using the kernel density estimation, is reported in the bottom left panel of Fig. 2. The most likely size suggested by this model is about 35 m.

2.1.5 Obliquity distribution

The obliquity of an asteroid can be determined from photometric measurements using the light-curve inversion method (Kaasalainen & Torppa 2001; Kaasalainen et al. 2001, 2002; Ďurech et al. 2015), or from radar observations (Hudson 1994; Hudson et al. 2000; Magri et al. 2007). When measurements are available, it is possible to model this parameter using a Gaussian distribution, with a mean corresponding to the nominal value and a standard deviation corresponding to the error. On the other hand, when it is unknown, which is the most common situation, we have to rely on general properties of the NEO population.

It is known that the NEO population includes an excess of retrograde rotators (La Spina et al. 2004), caused by the transport mechanism of objects from the main belt, driven by the Yarkovsky effect (Bottke et al. 2002; Granvik et al. 2018). A first attempt to derive an obliquity distribution for NEO was made by Farnocchia et al. (2013), who used a four-bin probability density function to match the observed obliquities. A more sophisticated approach was used in Tardioli et al. (2017), who tested different distributions. The best-fit solution is obtained for a quadratic distribution of cos γ, $$p(\cos \gamma )=a{\cos }^2\gamma +b\cos \gamma +c,$$(6)

where the parameters are a = 1.12, b = −0.32, c = 0.13. The PDF ofEq. (6) is shown in the bottom right panel of Fig. 2.

2.1.6 Fixed parameters

Semimajor axis

The uncertainty on the semimajor axis a of 2011 PT determined from astrometry, as reported in the JPL SBDB, is about 9.3 × 10⁻¹⁰ au: small changes of a within this uncertainty interval produce negligible changes in the predicted Yarkovsky drift, therefore this parameter was kept fixed to the nominal value.

Heat capacity

The heat capacity C depends on the physical characteristics of the object, but its value may vary within a factor of several (Delbó et al. 2015). Typical assumed values for kilometer-sized rocky or regolith-covered main-belt asteroids are in the range 600–700 J kg⁻¹ K⁻¹, while for iron-rich objects the value is lower, about 500 J kg⁻¹ K⁻¹ (see, e.g., Farinella et al. 1998). However, the heat capacity increases for increasing temperature, and its value may be higher for NEOs, but is assumed not to exceed ~1200 J kg⁻¹ K⁻¹. Ostrowski & Bryson (2019) reported the measured heat capacities for different types of meteorites. From their data, we can draw the conclusion that C is larger than 500 J kg⁻¹ K⁻¹ for stony meteorites, while it is smaller than this threshold value if there is a high component of metallic material. However, it is not clear how a distribution can be produced from the typically assumed value for kilometer-sized asteroid. Moreover, data on meteorites are not enough to produce a clear statistical distribution, they might be biased, and it is not clear how to correct their value when we consider small asteroids. For all these reasons, we kept the value of the heat capacity C fixed and produced the results for a few different reasonable values.

Emissivity and absorption coefficient

Ostrowski & Bryson (2019) reported the measurements of the emissivity for 61 meteorites, and all the objects but one have an emissivity between 0.9 and 1, with an average value of 0.984. Because we do not expect the Yarkovsky effect to vary much with respect to this parameter, we fixed ε to the mean value computed with meteorite measurements. The absorption coefficient was set to α = 1 in all our calculations.

2.1.7 Parameters modeled with errors in measurements

Rotation period

Kikwaya Eluo & Hergenrother (2015) determined the rotation period of 2011 PT using photometric observations, and it is P = 0.17 ± 0.05h, corresponding to about 11 min. Although the small uncertainty might produce very small variations on the predicted Yarkovsky drift, we took it into account by assuming a Gaussian distribution of the errors.

Yarkovsky drift measurement

The Yarkovsky drift has been determined by Greenberg et al. (2020) with a nominal value and an uncertainty of $${(\frac{\text{d}a}{\text{d}t})}_{\text{m}}=(-88.44\pm 14.6)\times {10}^{-4}{\text{auMy}}^{-1},$$

which we used to perform the Monte Carlo simulation. This parameter was also modeled assuming a Gaussian distribution of the errors. The JPL SBDB provides an automatic solution of the nongravitational transverse acceleration A₂ = (−2.237 ± 0.302) × 10⁻¹³ au d⁻², corresponding to (d a∕dt) = (−86.9 ± 11.7) × 10⁻⁴ au My⁻¹, while Del Vigna et al. (2018) reported a value of (d a∕dt) = (−91.3 ± 14.44) × 10⁻⁴ au My⁻¹, both nicely consistent with the value given by Greenberg et al. (2020).

2.2 Monte Carlo simulation

We generated samples $$\begin{array}{l}{\{(D_h,{\rho }_h)\}}_{h=1,\dots ,n_{p_V}},\\ {\left\{{\gamma }_j\right\}}_{j=1,\dots n_{\gamma }},\\ {\left\{P_k\right\}}_{k=1,\dots n_P},\\ {\left\{{(\text{d}a/\text{d}t)}_i\right\}}_{i=1,\dots n_{\text{Yarko}}},\end{array}$$

for the parameters modeled with a distribution, following the procedure described in the previous section, where $n_{p_V},n_{\gamma },n_P,n_{\text{Yarko}}$ denote the number of points in each corresponding sample. Diameter and density were generated in pairs using the same albedo values, which takes their correlation into account. In this manner, we avoided producing nonphysical combinations such as a small diameter (obtained by converting a high albedo) and low density (obtained by converting a low albedo). For a given Yarkovsky drift ${(\text{d}a/\text{d}t)}_i$ and a possible combinationof parameters (D_h, ρ_h, γ_j, P_k), we numerically inverted the relation $$\left(\frac{\text{d}a}{\text{d}t}\right)(a,D_h,{\rho }_h,K,C,{\gamma }_j,P_k,\alpha ,\epsilon )={(\frac{\text{d}a}{\text{d}t})}_i,$$(7)

computing the values of the thermal conductivity K for which the measured drift was achieved. In this manner, for each combination of parameters, we obtain $$K_1^{(h,j,k,i)},\dots ,K_n^{(h,j,k,i)},$$(8)

where n is the number of solutions of Eq. (7). Collecting all the solutions, we obtain a distribution for the possible values of the thermal conductivity K of the specific object, from which we can compute the corresponding probability density function using the kernel density estimator.

3 Results

3.1 Preliminary constraints on the thermal conductivity

A first indication that the thermal conductivity of 2011 PT might be low is given by simply computing the Yarkovsky drift for the nominal values of the parameters. For instance, fixing the diameter to D = 35 m and the heat capacity to C = 680 J kg⁻¹ K⁻¹, we computed the maximum negative Yarkovsky drift (thus selecting γ = 180°) as a function of the density ρ and the thermal conductivity K. The resulting drift values together with the three level curves corresponding to the measured nominal value and the 1σ levels are reported in Fig. 4.

Because we computed the maximum negative Yarkovsky drift, we can conclude from the plot of Fig. 4 that the values of (ρ, K), residing inside the area determined by the external level curve at 1σ level are all compatible with the drift measurements. For all the values of the density we obtained a low value of the thermal conductivity K, which is always lower than 0.3 W m⁻¹ K⁻¹, and decreases with increasing density. For instance, for ρ = 1500 kg m⁻³, we have $$0.00001{\text{Wm}}^{-1}{\text{K}}^{-1}\lesssim K\lesssim 0.074{\text{Wm}}^{-1}{\text{K}}^{-1},$$

while, doubling the density to ρ = 3000 kg m⁻³, we obtain $$0.00005{\text{Wm}}^{-1}{\text{K}}^{-1}\lesssim K\lesssim 0.0055{\text{Wm}}^{-1}{\text{K}}^{-1}.$$

These considerations already suggest a low thermal conductivity of the surface, which has to be analyzed in more depth taking the uncertainties in the parameters involved into account and by performing the Monte Carlo simulations explained in the previous section.

We note that the fast rotation of (499998) 2011 PT is a key element of the solution. This is because the maximum drift-rate d a∕dt is obtained when the diurnal F value is maximum in Eq. (A.2). The maximum occurs when Θ_d ≃ 1, that is, when the corresponding diurnal thermal parameter is on the order of unity. Because Θ_d depends on the ratio K∕P (see Eq. (A.6)), the optimal conductivity K scales as K ∝ P. This would mean that the whole pattern in Fig. 4 would slide toward higher conductivity values in the ordinate for more slowly rotating bodies. In the case of the very fast-rotating asteroid 2011 PT, K is constrained to values lower than ~0.1 W m⁻¹ K⁻¹.

Fig. 4 Estimated maximum negative Yarkovsky drift for 2011 PT obtained for C = 680 J kg−1 K−1, D = 35 m, and P = 0.17 h. The red level curves correspond to the nominal Yarkovsky drift (central level curve), and the drift obtained taking the uncertainties into account (external level curves).

Fig. 5 Derived thermal parameters of 2011 PT obtained from the Monte Carlo simulations for different values of the heat capacity C. The plots in the first column show the distribution of the thermal conductivity K, and the plots in the second column show the distribution of the thermal inertia Γ. The thermal inertia is obtained through $\Gamma =\sqrt{\rho KC}$, where ρ is the corresponding value of the density used to invert Eq. (7).

3.2 Results of the Monte Carlo simulations

To better quantify the low thermal conductivity of the surface of (499998) 2011 PT, we performed a Monte Carlo simulation with the method described in Sect. 2. We generated samples $\{(D_h,{\rho }_h)\},\left\{{\gamma }_j\right\},\left\{P_k\right\},\text{and}{\{\text{d}a/\text{d}t\}}_i$ of 50 000 points for diameterand density, and 10 000 points for obliquity, rotation period, and Yarkovsky drift. To perform the Monte Carlo experiment and generate theresults, we randomly chose one million combinations of the parameters. The output distribution for the thermal conductivity K was then fit using the kernel density estimation to obtain a smooth probability density function. We ran different simulations for some values of the heat capacity C, $$C=500,680,800,1000,\text{and}1200{\text{Jkg}}^{-1}{\text{K}}^{-1}.$$

Using these values, we explored the expected range of C for NEOs, from which we extrapolated an eventual trend. The results of the Monte Carlo simulations are shown in Fig. 5, together with the corresponding distribution obtained for the thermal inertia, $$\Gamma =\sqrt{\rho KC}.$$(9)

From the plots we note that we always obtained a bimodal distribution, with a first peak in the interval $$0.00001{\text{Wm}}^{-1}{\text{K}}^{-1}<K<0.001{\text{Wm}}^{-1}{\text{K}}^{-1},$$

and the second peak in the interval $$0.001{\text{Wm}}^{-1}{\text{K}}^{-1}<K<0.1{\text{Wm}}^{-1}{\text{K}}^{-1}.$$

The distributions of the thermal inertia Γ are also bimodal, with a first peak in the interval $$5{\text{Jm}}^{-2}{\text{K}}^{-1}{\text{s}}^{-1/2}<\Gamma <20{\text{Jm}}^{-2}{\text{K}}^{-1}{\text{s}}^{-1/2},$$

and the second peak in the interval $$75{\text{Jm}}^{-2}{\text{K}}^{-1}{\text{s}}^{-1/2}<\Gamma <150{\text{Jm}}^{-2}{\text{K}}^{-1}{\text{s}}^{-1/2}.$$

Moreover, as the heat capacity C increases, a third peak started to grow at high values of thermal conductivity, maintaining a low probability, however. The bimodal distribution is explained by the fact that for each set of parameters $\{(D_h,{\rho }_h)\},\left\{{\gamma }_j\right\},\left\{P_k\right\},\left\{{(\text{d}a/\text{d}t)}_i\right\}$, there are always two low thermal conductivity solutions of Eq. (1), as can be easily noted from Fig. 4. Furthermore, because the input distributions all have only one most likely value, most combinations of the input parameters selected during the Monte Carlo simulations are gathered nearby a most likely value, and this produces two clearly distinct peaks of low K. On the other hand, the third peak with a low probability is caused by other high thermal conductivity solutions of Eq. (1) appearing at high heat capacity, for a large diameter (around ~60 m) and a low density (around ~1100 kg m⁻³). These combinations of density and diameter are generated with a low probability because they correspond to the low-albedo case.

To give quantitative boundaries for the most likely values of thermal conductivity, we computed the probabilities $$\begin{array}{ll}{P}_1\hfill & =P(0.00001<K<0.001),\hfill \\ P_2\hfill & =P(0.001<K<0.1),\hfill \end{array}$$(10)

for the different simulations, as well as their sum and the two local maxima K₁, K₂. The results are reported in Table 2. From these values and from the plots of Fig. 5, we can conclude that the probability that the thermal conductivity is higher than 0.1 W m⁻¹ K⁻¹ was always almost zero, with slightly higher values for higher heat capacity values. The highest peak K₁ always lies in the first interval and is around 0.0001 W m⁻¹ K⁻¹, while the lowest peak K₂ belongs to the second, being ~0.005 W m⁻¹ K⁻¹. Moreover, we note that they both move to lower values of the thermal conductivity as the heat capacity C increases. The probability P₁ is always greater than the probability P₂, but only with a difference of a few percent. The values of P₁ and P₂ as a function of the heat capacity C can be linearly fit, and the results are show in Fig. 6, together with the 99% confidence interval region of the fit. As the heat capacity increases, the probability P₁ increases, while P₂ decreases with a more pronounced slope than P₁, which results in a decreasing trend in their sum (see Table 2). This fact was expected, however, because the third peak grows at high thermal conductivity.

Fig. 6 Linearfits of the probabilities P1 and P2 with respect to the heat capacity C. The dashed lines represent the 99% confidence interval region.

4 Discussion

4.1 Robustness and limitations of the model

The mathematical model of the Yarkovsky effect described in Appendix A is based on some assumptions that may affect the resulting thermal conductivity distribution. We discuss here the robustness and the limitations of the model, quantifying the effects of the orbital and physical features of 2011 PT that the model does not take into account.

4.1.1 Uncertainty due to the orbital eccentricity

The semimajor axis drift caused by the Yarkovsky effect that we describe analytically in Appendix A refers to an object revolving around the Sun on a circular orbit. However, the Yarkovsky effect may be increased on elliptic orbits (Spitale & Greenberg 2001), and it may cause an object to reach the measured migration rate even for higher values of the thermal conductivity. While significant changes are expected for orbits with high eccentricity, it is interesting to test the role of a moderately low orbital eccentricity in the case of (499998) 2011 PT.

The osculating semimajor axis drift caused by the Yarkovsky effect for an asteroid placed on an eccentric orbit is given by $$\frac{\text{d}a}{\text{d}t}=\frac{2}{n^{2}a}{f}_Y\cdot v,$$(11)

where n is the mean motion, v is the heliocentric orbital velocity, and f_Y is the instantaneous value of the Yarkovsky acceleration. The acceleration f_Y can be computed analytically, for instance, using the formulation given in Vokrouhlický et al. (2017), including both the diurnal and the seasonal effects (again with the simplifying assumption of linear boundary conditions). To compute the total semimajor axis drift, we averaged the values of the instantaneous drift rate of Eq. (11) over an orbital period numerically (rather than analytically, as before for the circular orbit).

To test the effects of the orbital eccentricity on our results, we compared the mean drift rates obtained for the circular orbit and the eccentric orbit by again performing a Monte Carlo simulation. We chose the eccentricity ~0.215 appropriate for (499998) 2011 PT. We used the same distributions as in Sect. 2 for the diameter, the density, the obliquity, and the rotation period. Moreover, we produced a sample of the thermal conductivity using the distributions obtained in Sect. 3. Using one million combinations of the parameters, we evaluated the two different drift rates and computed the relative difference.

We applied this procedure to (499998) 2011 PT using a fixed heat capacity of C = 680 J kg⁻¹ K⁻¹. The distribution of the percent absolute differences had two peaks, one around 1% and the other at 10%, and almost all the combinations gave a relative difference smaller than about 10%. Moreover, the peak at 1% was found to be corresponding to the lower peak in thermal conductivity, while the peak at about 10% was associated with the second peak in K, meaning that a higher conductivity produces larger changes in the Yarkovsky drift when the eccentricity of the orbit is taken into account. The contribution of the eccentricity to the total semimajor axis drift is therefore limited, and the circular orbit approximation should not produce misleading results for the thermal conductivity.

4.1.2 Further issues in modeling the Yarkovsky effect

We also briefly discuss further issues related to the Yarkovsky effect model that we outline in Appendix A and later used inSect. 3. The Yarkovsky effect in real-world objects is more complicated in several aspects that we did not explore above. These aspects include (i) nonlinearity of the boundary conditions for heat conduction at the surface, (ii) nonisotropy of the thermal emission by the surface facets (“thermal beaming”), and (iii) nonsphericity of the shape. While we did not test anyof them in detail, we note that some of them fortunately may affect the result in opposite ways, thus potentially compensating for each other (at least partially), while others may play a minimal role for (499998) 2011 PT.

As an example we note that the nonlinearity effects decrease the theoretically predicted semimajor axis drift rate d a∕dt in the relevant range of thermal conductivity values (see already Rubincam 1995; Čapek 2007). The reduction factor ranges between 0.7 and 0.9. On the other hand, thermal beaming effects always increase the semimajor axis drift rate by a factor ranging between 1.1 to 1.5 (see, e.g., Rozitis & Green 2012). The two effects thus have a tendency to compensate for each other, although in particular cases, some residual factor clearly may remain. At high obliquity, which we consider the most probable for 2011 PT, the nonsphericity decreases the predicted semimajor axis drift rate (see, e.g., Vokrouhlický 1998a). In this configuration the body presents a smaller cross section to the sunlight and thus processes less thermal energy than a spherical body. The reduction factor may be similar to that characterising the nonlinearity effects. While we do not have a shape model of 2011 PT, the low amplitude of the available photometric light curve in Kikwaya Eluo & Hergenrother (2015) suggests a near-spherical shape; or, better, taken the opposite way, there is no suggestion that the shape of 2011 PT is largely nonspherical. Thus it is possible that nonsphericity effects in the evaluation of the Yarkovsky drift rate are only small (modifying the predicted semimajor axis drift by not more than a factor ~ (0.8−0.9); see Fig. 3 in Vokrouhlický 1998a).

Another set of real-world details that were not included in our simple analytical model of the Yarkovsky effect has to do with physical parameters. In particular, they were assumed to be uniform and constant. In reality, things may be more complicated. For instance, the thermal conductivity, and to a lesser extent, also the heat capacity, may be temperature dependent. This is especially true for powdered mixtures of minerals (e.g., the proposed regolith-like surface) because energy transfer by conduction at the contact of the grains is efficiently accompanied by radiative transfer through cavities. The orbital eccentricity of (499998) 2011 PT fortunately is small, such that the thermal conditions at perihelion and aphelion do not change dramatically. For a low enough conductivity they may not exceed the temperature variations of the daily cycle. For instance, the equilibrium subsolar temperatures at pericenter and apocenter are ~320 and ~390 K, respectively. Data in Gundlach & Blum (2013) indicate that the conductivity would typically change in this temperature range by a factor ~1.5-1.7 only. Our assumed constant values therefore represent a median within this variation range. While the true effects of varying conductivity are not modeled here, we expect their effect to remain small.

Yet another issue is the assumption of homogeneity of the physical parameters. Consider, for instance, the density ρ discussed in Sect. 2.1.3. The way in which we introduced it corresponds to the volume-average (bulk) density. However, when we assume the existence of the regolith layer, we should ideally adopt a more complex model of the body: a sphere with a core with the bulk density and a lower density layer at the surface. For the fast rotation, the dominant diurnal effect remains constrained to the skin processes. Therefore the density in the right-hand side of Eqs. (A.4)–(A.6) should be the surface density, rather than the bulk density. Because the surface density is assumed to be lower (or equal at most) to the bulk density, the same thermal lag (and thus also dynamical effect) would need higher conductivity. The scaling is not clear a priori. When the regolith depth is much larger than the penetration depths of diurnal and seasonal thermal wave l_d and l_s, K recalibrates proportionally to the ratio of the bulk density over the surface density (this is because the thermal parameter in Eqs. (A.5) and (A.6) contains a Kρ factor). Because the plausible density ratiomay peak at ~2, the requiredconductivity values would be skewed to somewhat higher values. For instance, the nicely symmetric solution pattern seen in Fig. 4 would shift to higher conductivity values for densities above ~2000 kg m⁻³. However, even this approach may be too simple when l_d and, especially l_s, become comparable to or larger than the regolith depth. In this case, a full-fledged two-layer model would be needed. Vokrouhlický & Brož (1999) developed this approach for the seasonal component, but a complete formulation for the diurnal component is currently not available. At this moment we can therefore only conclude by noting that a fully consistent thermal model with the low-density regolith layer included may result in a moderate increase in the surface conductivity value.

Despite the modeling issues discussed above, it is worth mentioning that Chesley et al. (2014) predicted the bulk density of asteroid (101955) Bennu by combining a model for the Yarkovsky effect, based on known physical properties, and the measurements of the semimajor axis drift obtained from orbit determination. The more accurate bulk density value obtained during the OSIRIS-REx mission (Scheeres et al. 2019) revealed that the estimate by Chesley et al. (2014) was indeed correct. This successful prediction could suggest that the model is accurate enough for the analysis performed on 2011 PT.

4.2 Low thermal conductivity of (499998) 2011 PT

The results obtained with the Monte Carlo simulations suggest that the Yarkovsky drift determined for (499998) 2011 PT can be achieved preferentially for low values of thermal conductivity K, with two most likely values corresponding to the local maxima: the first around 0.0001 W m⁻¹ K⁻¹, and the second around 0.005 W m⁻¹ K⁻¹ (see Table 2). This conclusion can be explained in plain words by the combination of the fast rotation and the relatively strong Yarkovsky drift. If the thermal conductivity were high, the fast rotation would largely spread the thermal gradient on the surface of the asteroid, and a strong semimajor axis drift would be impossible to achieve. On the other hand, a low thermal conductivity would be able to keep a fair difference in temperature between the day-side and the night-side of the asteroid, that would be able to produce stronger semimajor axis drifts. The results obtained here indicate that values of thermal conductivity higher than 0.1 W m⁻¹ K⁻¹ are very unlikely, and the admitted values are instead fully compatible with the presence of a thermal insulating layer on the surface. This suggests that the surface is covered by a regolith-like material or by very thin dust. As a comparison, the values of the thermal conductivity of lunar dusty regolith, obtained during the Apollo missions, are between 0.0012 and 0.0035 W m⁻¹ K⁻¹ (Cremers & Birkebak 1971; Cremers & Hsia 1974), which are comparable with the values of the second peaks reported in Table 2. Assuming that thermal properties of the surface give indication about the grain size, Gundlach & Blum (2013) found that the radius of lunar regolith is estimated to be about 50 μm, which is compatible with the size distribution of sample lunar soil (Heiken et al. 1991). According to this, the value of thermal conductivity of the second peak suggests the presence of dusty regolith grains with a diameter of several dozen μm on the surface of (499998) 2011 PT.

At the same time, the thermal conductivity generally decreases with decreasing particle size, and thus the lower values obtained for the first peak might suggest an even smaller grain size. However, the dust layer in practice consists of particles of different sizes, but the thermal properties of particle size mixtures have not yet been well studied. Recently, Ryan et al. (2020) developed a 3D model to study the thermal conductivity of regolith by simulating heat flow in randomly packed spheres. These authors found that results predicted by simpler theoretical models are not reliable when the particles themselves are made of a material that itself has relatively low thermal conductivity, which could result in significant underestimation of particle sizes on asteroid surfaces.

It is worth noting that typical values for the thermal conductivity of bare rock are around 2–5 W m⁻¹ K⁻¹, but low values of about 0.1 W m⁻¹ K⁻¹ might be obtained also for a very high porosity of the material because measurements performed on meteorites show that K decreases as the porosity increases (Opeil et al. 2010, 2012; Ostrowski & Bryson 2019). Kikwaya Eluo & Hergenrother (2015), using photometry data, suggested a classification as a X-type asteroid for (499998) 2011 PT, which would imply a stony composition. Moreover, the high spin rate rules out the possibility that this object is a pure rubble-pile (Pravec & Harris 2000) because it is lower than the cohesionless spin-barrier of 2.2 h, hence the high-porosity hypothesis does not seem very likely. On the other hand, it is not fully clear what minimum level of cohesion forces are able to withstand moderate stresses existing on 2011 PT, what the associated level of maximum porosity is, and what this implies for thermal inertia of the surface. These issues deserve further analysis in the future.

We also note that the third peak at high thermal conductivity (see Fig. 5, bottom panels), arising at high heat capacity C, would either suggest a bare-rock composition or the presence of irony material on the surface. On the other hand, the presence of iron on the surface would cause a decrease in the heat capacity, as shown by measurements performed on meteorites collected by Ostrowski & Bryson (2019), which contradicts the fact that this peak grows at high C, leaving therefore only the bare-rock composition compatible with the third-peak solution. However, the occurrence of a high thermal conductivityassociated with the third peak is very unlikely because the probability is low. According to the above facts, the hypothesis that dusty regolith-like material is present appears to be most probable.

Harris & Drube (2016) analyzed the published values of thermal inertia Γ estimated through thermophysical models (Delbó et al. 2015), and suggested a possible correlation between the thermal inertia and the rotation period of NEOs, such that Γ decreases for increasing spin rates. On the other hand, Marciniak et al. (2019) found no evidence of this trend at very slow rotation rates. Assuming that the trend of Harris & Drube (2016) is valid also for shorter rotation periods, our result would support this hypothesis, although we underline that the data used by Harris & Drube (2016) and Marciniak et al. (2019) refer to objects with spin periods longer than 2 h, corresponding to the spin barrier of rubble piles. Moreover, it is difficult to draw any reliable conclusion about the trend for superfast rotators at the moment because there is evidence that such rotators could also have high thermal inertia, as the reported estimate for asteroid (54509) YORP indicates (Delbó et al. 2015).

4.3 Regolith production and retention on small bodies

The low thermal conductivity we found for the surface of (499998) 2011 PT could be explained with the presence of regolith. We discuss here the usual mechanisms that are capable to produce such material, and give some general hypotheses of how it may be retained on the surface of a fast rotator. It is generally thought that regolith can be produced in different ways. One of the main hypotheses is by micrometeoroid impacts, which produce small grains on the surface by means of fallback ejecta and by breaking up boulders (Hoerz et al. 1975; Horz & Cintala 1996). However, the velocity of crater ejecta is much higher than the escape velocity of small asteroids (Housen & Holsapple 2011), making this process not very feasible in this context. A second possibility is by thermal fatigue, which was proposed by Delbó et al. (2014). The temperature on the surface of an asteroid follows a diurnal cycle, and the stress resulting from sudden changes in temperature can cause damage to the surface material. This mechanism of rock weathering and cracking without ejection could be a possible explanation for the presence of regolith on small bodies. Moreover, thermal fatigue is a process that is independent of the size of the asteroid, as opposed to the micrometeoroid impact process. The cracking induced by the constant day-night cycle has recently been observed on asteroid (101955) Bennu with in-situ images taken during the OSIRIS-REx mission (Molaro et al. 2020). Another hypothesis is that the asteroid itself is active. An active asteroid holds some mechanism that causes mass shedding and ejection of particles from the surface, such as volatile releasing, rotational disruption, ice sublimation, and electrostatic lofting (Jewitt et al. 2015). Particles may be ejected from the surface with low enough velocity and consequently impact again on the asteroid, contributing to the cracking of superficial material due to micro-impacts. This process of particle ejection and fallback has been observed during the OSIRIS-REx mission (Lauretta et al. 2019b; McMahon et al. 2020), but it is not clear whether this activity is limited to primitive bodies of spectral types within the C-complex (which would not be applicable for 2011 PT).

To explain the presence of regolith on small and fast rotators, Sánchez & Scheeres (2020) developed a model to study the rotational conditions under which the surface material is retained, partially lost, or completely lost. The authors found that the covering material is preferentially lost in some regions of the body, and the area grows as the spin rate increases. Nevertheless, regions at high latitudes are able to retain regolith at arbitrarily high spins. This poses a clue that even small and very fast rotators can retain regolith-like grains and dust on their surface. However, it remains to be seen what happens to these polar seed-regions when the Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect (Vokrouhlický et al. 2015) causes changes in the rotation rate and potentially also orientation of the spin axis in space. A robust retention process of the formed grains on the surface may be needed.

Such a retention phenomenon may be provided by some type of cohesive forces, such as van der Waals forces, produced by the direct interaction of solar plasma with surface material, creating a superficial electric field. Cohesion may be able to hold fine grains together and stick them to the surface, and they have been demonstrated to be very high in lunar regolith (Mitchell et al. 1974). The strength of the cohesive force is represented by the bond number B, which is given by $$B={10}^{-5}{g}_{\text{A}}^{-1}{d}^{-2},$$(12)

where g_A is the ambient gravity and d is the grain diameter (Scheeres et al. 2010; Rozitis et al. 2014). The stability of surface material needs a bond number larger than ten (Rozitis et al. 2014), hence we can set upper constraints on the diameter of grains. Assuming the nominal rotation period and an obliquity of 170°, we computed the ambient gravity g_A as a function of the diameter D of the asteroid andthe density ρ, using a spherical model for the shape. Then we computed the upper limit for the grain diameter d setting a bond number equal to ten in Eq. (12), obtaining a function of R and ρ. The result is shown in Fig. 7.

The maximum grain diameter obtained is between 1 and 4 mm. Moreover, these bounds are essentially independent of the obliquity γ between 90° and 180°. These estimates give upper bounds for the size of particles that can be held on the surface by means of cohesive forces, and they are well above the ~10–100 μm size of the lunar regolith expected from the second peak of thermal conductivity. Therefore, the hypothesis that cohesive forces exist is a reasonable explanation, which could further be explored using better shape and thermophysical models when additional appropriate observational data are acquired.

Fig. 7 Maximum grain size as a function of the density ρ and the diameter D. Assuming a bond number equal to ten in Eq. (12), the maximum grain size allowed refers to the value obtained for the maximum ambient gravity gA over the asteroid surface.

5 Summary and conclusions

We developed a statistical method for constraining the thermal conductivity K of (499998) 2011 PT, anNEO with a diameter of about 35 m that rotates very fast, with a period of 11 min. The method is based on the comparison of the measured and theoretically predicted Yarkovsky drift, which depends on the orbital and physical parameters of the asteroids. Modeling distributions for the input parameters, we used a Monte Carlo simulation to produce a value distribution of the thermal conductivity that is compatible with the measured Yarkovsky drift. When the sample was obtained, a probability density function for K was computed, and it was used to deduce the most likely values. The method developed here is very flexible with respect to the known and unknown input parameters, and it can easily be applied to asteroids with different sizes and rotation periods. It may be alsocustomized on a case-by-case basis, in order to obtain the best possible solution in each individual case.

The obtained probability density functions for the thermal conductivity of 2011 PT always have two local maxima, one around 0.0001 W m⁻¹ K⁻¹, and the another one around 0.005 W m⁻¹ K⁻¹. These two values are the most likely ones, with the first having a slightly higher probability. With the high probability of at least 95%, we constrained the thermal conductivity to be in the range 0.00001 W m⁻¹ K⁻¹ < K < 0.1 W m⁻¹ K⁻¹. This is the first time that the low thermal conductivity solution for a small and superfast rotating asteroid has been obtained with such high probability. The low-thermal conductivity can be interpreted as the presence of regolith-like material or thin dust on the surface of (499998) 2011 PT, opening many related scientific questions, and in particular, how these super-fast rotators are able to hold this material on their surface.

To better explore this possibility, however, more observations and characterizations of small asteroids are needed. A more precise answer can be obtained by planning space missions to one or more such NEOs, eventually developping a technique for sampling some surface material of fast rotators that can be brought back to Earth for an accurate analysis. In this respect, the proposed extension of the JAXA Hayabusa 2 mission to visit the asteroid³ 1998 KY26, which is about 30 m in diameter and rotates with a period of about 10 min, may provide very interesting insights on the structure and the composition of very small and fast-rotating objects.

Acknowledgements

We thank the anonymous referee for the suggestions that helped us to improve the manuscript. M.F. and B.N. have been supported by the MSCA-ITN Stardust-R, Grant Agreement n. 813644 under the European Union H2020 research and innovation program. The work of D.V. was partially funded by the Czech Science Foundation (grant 18-06083S). This research has made use of data and/or services provided by the International Astronomical Union’s Minor Planet Center. Pan-STARRS is supported by the National Aeronautics and Space Administration under Grant No. 80NSSC18K0971 issued through the SSO Near Earth Object Observations Program.

Appendix A Analytical formulation of the Yarkovsky effect

We recall here the analytical formulation of the Yarkovsky effect given in, for instance, Vokrouhlický (1998a, 1999). We assume a spherical body moving on a circular orbit around the Sun that rotates along a fixed axis. By linearizing the boundary condition of the heat diffusion problem, the drift in semimajor axis a caused by the Yarkovsky effect is given by two distinct components: the seasonal effect $${(\frac{\text{d}a}{\text{d}t})}_{\text{s}}=\frac{4\alpha }{9}\frac{\Phi }{{\omega }_{\text{rev}}}F({R^{\prime }}_{\text{s}},{\Theta }_{\text{s}}){\sin }^2\gamma ,$$(A.1)

and the diurnal effect $${(\frac{\text{d}a}{\text{d}t})}_{\text{d}}=-\frac{8\alpha }{9}\frac{\Phi }{{\omega }_{\text{rev}}}F({R^{\prime }}_{\text{d}},{\Theta }_{\text{d}})\cos \gamma .$$(A.2)

In the above equations, α is the absorption coefficient of the surface, Φ is the radiation pressure coefficient (e.g., Vokrouhlický et al. 2015), ω_rev is the orbital frequency (mean motion), and γ is the spin-axis obliquity. ${R^{\prime }}_{\text{s}}$ and ${R^{\prime }}_{\text{d}}$ are the scaled (nondimensional) values of the radius R, defined as $${R^{\prime }}_{\text{s}}=\frac{R}{l_{\text{s}}},{R^{\prime }}_{\text{d}}=\frac{R}{l_{\text{d}}},$$(A.3)

where l_s, l_d are the penetration depths of the seasonal and diurnal thermal wave given by $$l_{\text{s}}=\sqrt{\frac{K}{\rho C{\omega }_{\text{rev}}}},l_{\text{d}}=\sqrt{\frac{K}{\rho C{\omega }_{\text{rot}}}}.$$(A.4)

We note that l_s and l_d depend on the thermal conductivity K, the heat capacity C, and the density ρ of the asteroid. Additionally, the two length scales depend on these frequencies: (i) the spin frequency ω_rot in the case of the diurnal effect, and (ii) the orbital frequency ω_rev in the case of the seasonal effect. The thermal parameters Θ_sand Θ_d also depend on the physical and thermal characteristics of the object, and they are defined as $$\begin{array}{lll}{\Theta }_{\text{s}}\hfill & =\hfill & \frac{\sqrt{\rho KC{\omega }_{\text{rev}}}}{\epsilon \sigma T_{\star }^3},\hfill \\ {\Theta }_{\text{d}}\hfill & =\hfill & \frac{\sqrt{\rho KC{\omega }_{\text{rot}}}}{\epsilon \sigma T_{\star }^3},\hfill \end{array}$$

where σ is the Stefan-Boltzmann constant, ε is the emissivity and T_⋆ is the subsolar temperature, defined by $\epsilon \sigma T_{\star }^4=\alpha {\mathcal{E}}_{\star }$, with ${\mathcal{E}}_{\star }$ being the solar radiation flux at the distance of the body. The function F in Eqs. (A.1) and (A.2) depends on both the corresponding scaled radius and the thermal parameter, and it determines the total magnitude of the Yarkovsky drift. The definition of F is given by $$F(R^{\prime },\Theta )=-\frac{k_1(R^{\prime })\Theta }{1+2k_2(R^{\prime })\Theta +k_3(R^{\prime }){\Theta }^2}.$$(A.7)

The coefficients k₁, k₂, and k₃ are positive analytical functions of the scaled radius, and their precise definition can be found, for instance, in Vokrouhlický (1998b, 1999). It is worth noting that the seasonal component always produces an inward migration, while the direction of migration for the diurnal component depends on the obliquity γ.

References

All Tables

All Figures

	Fig. 1 Measured spin periods for NEOs as of September 2020 and taken from LCDB (see Warner et al. 2009).
In the text

	Fig. 2 Probability density function for the albedo pV (top left panel), the density ρ (top right panel), the diameter D (bottom left panel), and the obliquity γ (bottom right panel) computed for asteroid (499998) 2011 PT.
In the text

	Fig. 3 Fit of the asteroid porosity with a gamma distribution using data from Carry (2012).
In the text

	Fig. 4 Estimated maximum negative Yarkovsky drift for 2011 PT obtained for C = 680 J kg−1 K−1, D = 35 m, and P = 0.17 h. The red level curves correspond to the nominal Yarkovsky drift (central level curve), and the drift obtained taking the uncertainties into account (external level curves).
In the text

Fig. 5 Derived thermal parameters of 2011 PT obtained from the Monte Carlo simulations for different values of the heat capacity C. The plots in the first column show the distribution of the thermal conductivity K, and the plots in the second column show the distribution of the thermal inertia Γ. The thermal inertia is obtained through $\Gamma =\sqrt{\rho KC}$, where ρ is the corresponding value of the density used to invert Eq. (7).

In the text

	Fig. 6 Linearfits of the probabilities P1 and P2 with respect to the heat capacity C. The dashed lines represent the 99% confidence interval region.
In the text

	Fig. 7 Maximum grain size as a function of the density ρ and the diameter D. Assuming a bond number equal to ten in Eq. (12), the maximum grain size allowed refers to the value obtained for the maximum ambient gravity gA over the asteroid surface.
In the text