TP1

Summary:

We investigated the Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect [Rubincam 2000] acting on asteroid (162173) Ryugu. To accurately simulate the thermal condition on Ryugu, we employed a thermophysical model recently developed by Nakano and Hirabayashi [2021], which considers self-heating effect and scattering of radiation and solves the 3-dimensional heat equation via Finite Element Modeling (FEM) approach. Our result indicates that on 1 August 2018 TDB, the thermal torque (YORP torque) was acting in a direction that accelerates Ryugu’s spin rate.

Introduction:

Dynamics of small bodies among the near-Earth and main belt asteroids are controlled by secular effects due to thermal radiation forces and torques, known as the Yarkovsky and YORP effects [Rubincam 2000; Vokrouhlický and Čapek, 2002, Bottke et al., 2006]. Recent increase in the number of observed top-like shape objects suggests that the YORP effect is a key mechanism for their spin rates acceleration (or deceleration), which can induce mass-shedding or catastrophic disruptions [e.g., Holsapple, 2010; Hirabayashi, 2015]. Thus, modeling the Yarkovsky and YORP effects based on a proper thermophysical technique is one of the crucial steps to understand their formation and evolution mechanisms. Here, we investigate the YORP effect on asteroid (162173) Ryugu by using a thermophysical model recently developed by Nakano and Hirabayashi [2021].

Thermophysical modeling of Ryugu:

Unlike earlier thermophysical models which typically solve the 1-dimensional heat equation, the FEM approach thermophysical model solves the 3-dimensional heat equation. It considers the self-heating and scattering of radiation, similar to other models [e.g., Rozits and Green, 2011; Davidsson and Rickman, 2014]. The surface boundary condition includes: direct solar flux U, diffuse solar radiation flux W, direct self-heating u, and diffuse self-heating w [Davidsson and Rickman, 2014]. To efficiently describe the surface/subsurface thermal evolution, “layer” FEM shape model, where only the layers of tetrahedral meshes exist for the surface and subsurface, and the inside is left as cavity, is used. (See Nakano and Hirabayashi [2021] for more detail.)

We retrieved Ryugu’s ephemeris and the sun’s position vector as seen from Ryugu on 1 August 2018 TDB, using a SPICE kernel for Ryugu. The orbital elements and physical parameters used in this investigation are listed in Table 1. The layer FEM shape model was constructed based on SFM_49k_v20180804 [Watanabe et al., 2019], but the number of facets was decreased from 49k to 4.9k. Figure 1 shows the surface temperature distribution on Ryugu obtained after 25 spins (≈190 hr) at the fixed distance from the Sun. The surface temperature distribution is comparable to that observed by Hayabusa2 spacecraft [Okada et al., 2020].

The YORP effect on Ryugu:

Once the surface temperature is obtained, a recoil force df due to reflected solar and thermally radiated photons for each facet can be computed [Rozits and Green, 2012]. Assuming Lambert (isotropic) emission from the surface, it is given by [Vokrouhlický and Čapek, 2002; Rozits and Green, 2012]:

where ε is the emissivity of the surface, σ is the Stefan-Boltzmann constant, T is the surface temperature, G is reflected scattered flux, c is the speed of light, and dS is the facet area. The net torque τ acting on Ryugu can be computed by taking a cross product of the facet position r and the recoil force df for each facet and integrating it over the whole surface:

The net torque computed by the above equation can be decomposed into three meaningful YORP torque components of our interest: a spin period change component τ_ω, a spin pole obliquity change component τ_ε, and a precession of rotation τ_ψ (see e.g., Bottke et al., 2006 for the equations). Figure 2 shows the time history of YORP torque components over one Ryugu’s spin. The one-spin averaged YORP torque components are found to be:

Discussion:

Due to its top-like asymmetric shape, Ryugu should be susceptible to the YORP effect. Our investigation with the 3-dimensional thermophysical modeling technique clearly showed that Ryugu does experience non-negligible YORP effect. While the YORP torque components change their signs during one spin (Figure 2), the one-spin averaged values were all found to be positive. In particular, the positive  implies an increase of Ryugu’s spin rate. This contradicts with the current hypothesis that Ryugu spun faster in the past and later slowed down to the current spin state [Watanabe et al., 2019]. In this investigation, however, we only explored Ryugu’s thermal condition at a particular epoch. Therefore, we cannot infer Ryugu’s secular dynamical evolution. The YORP effect may be highly sensitive to small topographic features on the surface [Statler, 2009]. The surface roughness may also affect thermal evolutions [Shimaki et al., 2020]. We plan to expand our modeling technique to investigate such effects and to further constrain the Yarkovsky and YORP effects on small bodies.

Table 1. List of parameters used in this investigation

Figure 1. Surface temperature distribution on 1 August 2018 TDB

Figure 2. YORP torque components versus time

Acknowledgments: RN and MH acknowledge support from NASA/Solar System Workings (NNH17ZDA001N/80NSSC19K0548) and Auburn University/Intramural Grant Program.

References:

Bottke et al. (2006) Annu. Rev. Earth Planet. Sci., 34; Davidsson and Rickman (2014), Icarus, 243; Hirabayashi (2015) MNRAS, 454; Holsapple (2001) Icarus, 154; Nakano and Hirabayashi (2021) 52^nd LPSC, 1292; Okada et al. (2020) Nature, 579; Rozits and Green (2011) MNRAS, 415; Rozits and Green (2012) MNRAS, 423; Rubincam (2000) Icarus, 148; Shimaki et al. (2020) Icarus, 348; Statler (2009) Icarus, 202; Vokrouhlický and Čapek (2002) Icarus, 159; Watanabe et al. (2019) Science, 364

How to cite: Nakano, R. and Hirabayashi, M.: Investigation of the YORP effect on asteroid (162173) Ryugu – An application of FEM approach thermophysical model, Europlanet Science Congress 2021, online, 13–24 Sep 2021, EPSC2021-441, https://doi.org/10.5194/epsc2021-441, 2021.

Abstract

We present here a method for the retrieval of the internal density distribution of a small body, via a least squares inversion of its gravity field and rotational state. Multiple solutions are averaged in order to limit the effect of the a priori density distribution. Simulations show successful density map estimations even with little to no initial information on the interior structure of the body.

1. Introduction

The gravity field of a body is a property among the most readily obtainable by a spacecraft mission. Be it by employing radio-tracking, optical, or gravimetric data, several missions have provided and will provide gravity expansions of degree 2 and higher for small bodies. In turn, the extended gravity field is an effect of the distribution of the mass inside the body. However, the inverse problem of estimating the density distribution from the gravity field is not a trivial one, and generally with non-unique solutions [1].

Here, we propose a method to obtain a unique solution for the internal mass distribution of a small body by restricting our estimation to zones of constant density, which are limited in number by the accuracy of the gravity model [2]. Where available, additional observables such as the libration amplitudes could help in the estimation of the interior structure. Averaging of the solutions obtained for different subdivisions of zones provides then the estimated density map [3].

2. Methods

The code builds from available software which allows to compute the gravity harmonics of a body approximated by a collection of cubes, each with a given density [4]. Here, the cubes are grouped in subregions of constant density, as many as there are measured coefficients. Each subregion is assigned an initial density, from which the software computes gravity coefficients matching in degree and order those observed. The difference between the coefficients thus computed and the observed ones is then minimized in a least-squares sense [5]. The solution of the regularized least-squares inversion is the density for each of the initial zones, along with its uncertainty. However, this solution is heavily dependent on the initial subdivision of the body into zones. Hence, a wide range of possible initial zones is explored, and the corresponding density solutions averaged together. 

From such an average solution, a new set of subregions is generated, with the help of a blob detection algorithm [6]. Finally, the gravity coefficients are inverted using this new zone subdivision. This step is necessary because, while providing a good distribution of the density zones, the density values of the averaged map are biased by solutions coming from inaccurate initial zones subdivisions. 

The data used in the inversion are the spherical harmonics coefficients of the gravity field expansion. The advantage of these observables is that they can be measured without the need for specific payload. Although the method can be easily adapted to process surface gravity measurements, the larger errors on the surface gravity coming from the cubes discretization make it less appealing for these type of data.

3. Discussion and outlook

Figure 1 shows the results obtained when applying the method to a fictitious density distribution for asteroid Bennu. The nominal density distribution, shown in Figure 1a, was used to produce a synthetic gravity field expansion of degree 11. An optimistic noise profile was assumed for these simulated observables, with a 1% error on each gravity coefficient. The spherical harmonics coefficients were then inverted via the iterative process described in Section 2, resulting in the density distribution shown in Figure 1b. 

With a gravity field of this resolution, the inversion method is clearly able to distinguish the 4 density anomalies inside the body. Still, the comparison between the nominal and the estimated density models reveals the current limitations of the method in accurately reconstructing the size of the anomalies or their density — which is equivalent, given the constraint on the total mass. 

Once a better agreement between the retrieved and the nominal models will be achieved with the current (high) resolution of the gravity harmonics, the synthetic gravity fields will be truncated at lower degrees, making them more realistic as outputs of spacecraft missions to small bodies. Future simulation studies will also make use of more realistic noise profiles, with the errors on the observed coefficient increasing as a function of the degree.

(a)

(b)

Figure 1: (a) density model used to simulate a degree-11 gravity field; (b) density distribution retrieved from said gravity field, without any initial assumptions about the true distribution

Acknowledgements

This work was financially supported by the French community of Belgium within the framework of a FRIA grant, and by the Belgian PRODEX program, managed by the European Space Agency in collaboration with the Belgian Federal Science Policy Office.

References

[1]  P. Tricarico. Global gravity inversion of bodies with arbitrary shape. Geophysical Journal International, 195(1):260–275, 2013.
[2]  Y. Takahashi and D. J. Scheeres. Morphology driven density distribution estimation for small bodies. Icarus, 233:179–193, 2014.
[3]  L.-I. Sorsa, M. Takala, P. Bambach, et al. Tomographic inversion of gravity gradient field for a synthetic Itokawa model. Icarus, 336:113425, 2020.
[4]  S. Le Maistre, A. Rivoldini, and P. Rosenblatt. Signature of Phobos’ interior structure in its gravity field and libration. Icarus, 321:272–290, 2019.
[5]  D. J. Scheeres, B. Khushalani, and R. A. Werner. Estimating asteroid density distributions from shape and gravity information. Planetary and Space Science, 48(10):965–971, 2000.
[6]  S. v. d. Walt, J. L. Schönberger, J. Nunez-Iglesias, et al. scikit-image: image processing in Python. PeerJ, 2:e453, 2014.

How to cite: Caldiero, A., Le Maistre, S., and Dehant, V.: Estimation of the interior density of a small body given its gravity field, Europlanet Science Congress 2021, online, 13–24 Sep 2021, EPSC2021-756, https://doi.org/10.5194/epsc2021-756, 2021.

The binary asteroid 65803 Didymos-Dimorphos is the target of NASA’s DART (Cheng et al., 2018) and ESA's Hera missions (Michel et al., 2018). Hera will arrive at the asteroid system several years after the DART impact. It will carry out a detailed characterisation of the asteroids’ overall properties and will measure the outcome of the DART impact on Dimorphos. 

Didymos, the primary body, is a fast spinning asteroid with a period of only 2.26 hours, which is very close to the critical spin limit of 2.2 hours for asteroids bigger than approx. 200 m (Walsh et al., 2018). The interior structure of Didymos is crucial for our understanding of its stability. Without cohesion and with an estimated bulk density of 2.1 g/cc, Didymos is not able to keep its shape stable (Holsapple, 2001; Zhang et al., 2017) and thus, cohesive forces might be present in its structure.

The interior strength properties of asteroids determine to a large degree their collisional evolution. For small bodies, even a small level of cohesion can significantly affect the outcome of an impact  (Raducan and Jutzi, 2021). 

Previous studies have applied N-body and soft-sphere discrete element (SSDEM) codes to study the structural stability of rubble pile asteroids (Sanchez and Scheeres, 2012; Zhang et al., 2017, 2021; Ferrari and Tanga, 2020). Recently, the stability and failure modes of such objects have been studied using finite element (FEM) (Hirabayashi et al., 2020). Here we use Bern’s Smooth Particle Hydrodynamics (SPH) code (Jutzi et al., 2008; Jutzi 2015) to simulate the rotating asteroid and to investigate its physical properties. To do this, we set up Didymos according to its recent shape model (provided by the Hera Didymos Reference Model) with a bulk density of 2.1 g/cc and we spin up the body to a rotation period of 2.26 hours.

For the preliminary simulations presented here, we assume that Didymos has a homogeneous structure. We perform simulations using a range of values for the cohesion and the internal friction coefficient of the asteroid, and then for each case we evaluate Didymos’ stability.

An example of such a simulation is shown in Figure 1. For a cohesion of 100 Pa and a friction coefficient of 0.4, we find that Didymos’ shape is stable and only very little deformation takes place, as indicated by the small values in the total accumulated strain. In Figure 2, we show the case with a lower cohesion of 0.1 Pa and a coefficient of friction of 0.4. In this case, Didymos experiences large-scale deformations. To quantify the degree of deformation, we compute the cumulative strain distributions for both cases (Figure 3). For the 100 Pa cohesion case, > 90 % of the total mass experiences a total strain of < 0.03 and therefore this case is considered to be stable. On the other hand, for the case with 0.1 Pa cohesion, > 50 % of the total mass experiences a total strain of > 0.5, consistent with the large scale deformations observed in Figure 2. 

Our initial results suggest that a cohesion of ~ 10-100 Pa is required to maintain structural stability, which is consistent with recent SSDEM modeling results (Zhang et al. 2021). 

However, so far only homogeneous structures have been explored. We plan to investigate also more complex interior structures that may lead to different outcomes in terms of required cohesion to obtain an overall structural stability.