---

A&A 507, 1073-1081 (2009)

The YORP effect on 25 143 Itokawa

S. Breiter¹ - P. Bartczak¹ - M. Czekaj^2,1 - B. Oczujda¹ - D. Vokrouhlický³

1 - Astronomical Observatory, Adam Mickiewicz University, Soneczna 36, 60-286 Poznan, Poland
2 - Departament d'Astronomia i Meteorologia, Universitat de Barcelona, Av. Diagonal 647, 08028 Barcelona, Spain
3 - Institute of Astronomy, Charles University, V Holesovickách 2, 18000 Prague 8, Czech Republic

Received 20 May 2009 / Accepted 21 August 2009

Abstract
Context. The asteroid 25143 Itokawa is one of the candidates for the detection of the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect in the rotation period. Previous studies were carried out up to the 196 608 facets triangulation model and were not able to provide a good theoretical estimate of this effect, raising questions about the influence of the mesh resolution and the centre of mass location on the evolution the rotation period.
Aims. The YORP effect on Itokawa is computed for different topography models up to the highest resolution Gaskell mesh of 3 145 728 triangular faces in an attempt to find the best possible YORP estimate. Other, lower resolution models are also studied and the question of the dependence of the rotation period drift on the density distribution inhomogeneities is reexamined. A comparison is made with 433 Eros models possessing a similar resolution.
Methods. The Rubincam approximation (zero conductivity) is assumed in the numerical simulation of the YORP effect in rotation period. The mean thermal radiation torques are summed over triangular facets assuming Keplerian heliocentric motion and uniform rotation around a body-fixed axis.
Results. There is no evidence of YORP convergence in Gaskell model family. Differently simplified meshes may converge quickly to their parent models, but this does not prove the quality of YORP computed from the latter. We confirm the high sensitivity of the YORP effect to the fine details of the surface for 25 143 Itokawa and 433 Eros. The sensitivity of the Itokawa YORP to the centre of mass shift is weaker than in earlier works, but instead the results prove to be sensitive to the spin axis orientation in the body frame.
Conclusions. Either the sensitivity of the YORP effect is a physical phenomenon and all present predictions are questionable, or the present thermal models are too simplified.

Key words: celestial mechanics - minor planets, asteroids - methods: numerical

1 Introduction

The success of the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect after the article of Rubincam (2000) is by no means surprising. After the epoch of considering weak tidal damping torques (Burns & Safronov 1973) or planetary flybys (e.g. Scheeres et al. 2000) as the only effects that systematically alter rotation states of asteroids, astronomers now have an attractive mechanism capable of shaping the spin statistics surprisingly quickly. Indeed, the timescale of the YORP effect, resulting from nonisotropic thermal radiation, proved to be short enough to allow its direct observational confirmation.

The candidate asteroids should be small enough, sufficiently irregular in shape and have well determined rotational elements over at least a few years. In 2007, three objects were reported to reveal a detectable change of the rotation period coherent with the YORP estimates: , subsequently renamed 54 509 YORP (Lowry et al. 2007; Taylor et al. 2007), 1862 Apollo (Kaasalainen et al. 2007), and 25 143 Itokawa (Kitazato et al. 2007). Although Itokawa's YORP detection was withdrawn soon after its announcement, Durech et al. (2008a) reported the observational confirmation of the YORP effect on 1620 Geographos.

Curiously, Itokawa (formerly ) had looked like a perfect candidate for the YORP detection, with its small (0.3 km) diameter, observationally convenient orbital semi-axis of , irregular shape and well determined pole position and rotation period. Vokrouhlický et al. (2004) published an optimistic prediction that the change of rotation rate  should induce the quadratic phase drift detectible in future oppositions of 2004 and 2006, as well as from the expected high precision observations from the Hayabusa spacecraft expected in 2005. But the next years observational efforts indicated that the asteroid's rotation does not agree with the YORP predictions; not only with the original Vokrouhlický et al. (2004) results, but also with more recent ones of Scheeres et al. (2007). According to Durech et al. (2008b) if there is any sign of the YORP-induced rotation rate evolution for Itokawa, it must be hidden within the current observational error bounds. More precisely, a detectible should be at least    - the value that raises the  of the present observations by at most . Theoretically predicted values are significantly higher than that, but they have not been detected in reality.

Whether fulfilled or not, expectations of detecting the YORP on Itokawa required apparently routine computations of the theoretically predicted  . Vokrouhlický et al. (2004) used early shape models derived from radar observations (Ostro et al. 2004) and lightcurve inversion (Kaasalainen et al. 2003). Both models were quite small (4036 and 2040 facets respectively) and convex, leading to positive  (the negative value of Vokrouhlický et al. 2004; was due to a sign error, as shown by Scheeres et al. 2007). But then the rich flux of data from Hayabusa has flooded the existing YORP modeling software with topography models of increasing complexity. Scheeres et al. (2007) were capable of handling the 50 696 triangles of the preliminary Gaskell et al. (2006b) model. Comparing the results with subsequently improved radar shapes of Ostro et al. (2005), they noticed that the appearance of concavities in the Itokawa shape models leads to negative  with increasing magnitude and the values of YORP derived from increasingly finer topography models do not seem to stabilize. Durech et al. (2008b) went a step further, computing the YORP for a larger model, reaching 196 608 facets resolution and including 1D thermophysical modelling at the same time. In agreement with the trend observed by Scheeres et al. (2007), the resulting was negative and greater in magnitude than its 49 152 facets counterpart. Hoping that there might be no reason to use more detailed models, they derived a sequence of 21 triangular meshes by simplifying Gaskell's 196 608 mesh to 1000 facets. But the resulting sequence of the YORP effect values revealed no sign of ``saturation'' while approaching the finest mesh, indicating that the last value of   is not the ultimate. Of course, there was a slight chance that at higher resolution the trend might be reversed, lifting  back to the interval consistent with the observational bounds, but higher resolution models have been beyond the capability of the authors' computing power. Within a computationally cheap, but rather unrealistic ``pseudo-convex'' algorithm that ignores the shadowing effects, Scheeres et al. (2007) computed the effect due to the 196 608 facet model, and Scheeres & Gaskell (2008) went up to the 3 million facets mesh. According to their conclusions, there was no sign of convergence even in this simple case.

An interesting path to lower the computed YORP strength for Itokawa to the observational constraints was proposed by Scheeres & Gaskell (2008), who pointed out that shifting the centre of mass by  may cancel the component of the mean YORP torque responsible for the spin rate evolution. Such a translation may be justified by inhomogeneous density - quite probable in a rubble-pile body like Itokawa. The influence of the centre of mass location for the YORP effect was first signaled in the analytical solution of Nesvorný & Vokrouhlický (2008), but this feature did not look like an important factor for the mildly irregular objects discussed in their model.

The present paper continues the computations of the theoretical YORP effect on Itokawa. Using our optimized code run in parallel on about 12 computers, we have been able to compute the YORP effect on Itokawa up to the most detailed Hayabusa model of over 3 million triangular faces. The results are summarized in Sect. 2. We also recomputed the chain of models discussed by Durech et al. (2008b) and compare it with alternative families of smaller models derived using various mesh simplification (decimation) strategies. For comparison, we compute the YORP effect on 433 Eros - a second object with a mega-facet topography model.

The hypothesis of Scheeres & Gaskell (2008) is discussed in Sect. 3. We account for the translation of the centre of mass (COM) using a less approximate approach. Moreover, we consider the influence of the rotation of principal axes rejected a priori by Scheeres & Gaskell (2008).

  2 YORP effect values from different topography models

  2.1 Computational model

All known recipes for computing the YORP torques can be roughly divided in two main groups: most of them use a discrete mesh, summing the contributions from each triangular face (Scheeres 2007; Vokrouhlický & Capek 2002; Scheeres & Mirrahimi 2008; Capek & Vokrouhlický 2004; Statler 2009), whereas others use a continuous shape model described in terms of spherical harmonics (Nesvorný & Vokrouhlický 2008; Breiter & Michalska 2008; Mysen 2008). The latter group cannot be applied to Itokawa, because its irregular shape, not being starlike, cannot be expressed in terms of spherical harmonics. In these circumstances, to compute  as a function of time t, we use the standard ``per triangle'' summation

where c is the speed of light,  is the Stefan-Boltzmann constant, and I₃ is the maximum principal moment of inertia along the axis determined by the unit vector . The same vector coincides with the spin vector  , because we assume the principal axis rotation mode without wobbling. The summation index j runs over all triangular faces with centroids , outward oriented surface elements  and temperature . Assuming a homogeneous object we let all triangles share the common emissivity  .

The temperature distribution on the body surface should be modeled by solving the heat diffusion equation with boundary conditions at each triangular face

resulting from the energy conservation principle. The left hand side is a sum of thermally re-radiated power flux and of the heat convected into an asteroid's interior; the right-hand side  is the irradiation of the jth triangle by the sunlight. None of the published YORP models uses the exact form of Eq. (2) as the boundary conditions. Rubincam (2000), Vokrouhlický & Capek (2002) and Statler (2009) neglect the conductivity K obtaining the ``Rubincam approximation''

Note that in this case the surface temperature is directly given by the boundary conditions and the heat diffusion equation itself can be discarded. Capek & Vokrouhlický (2004), Scheeres (2007), Scheeres & Mirrahimi (2008), Mysen (2008), Nesvorný & Vokrouhlický (2008), and Breiter & Michalska (2008) use various approximations to account for the nonzero conductivity. Their thermal models use a ``plane parallel'' approximation reducing the problem to a 1D heat diffusion model. The asteroid is treated like a union of disjoint, noninteracting, one end heated semi-infinite rods.

All the enumerated models indicate that - in contrast to the YORP effect in obliquity - the secular evolution of  does not depend on conductivity in the plane parallel thermal model. In these circumstances, being interested only in the YORP effect on the rotation period, we have assumed K=0 - a choice that avoids having to solve the heat diffusion equation numerically. Using Rubincam's approximation, one has to face only one serious problem: computing the irradiation. If the direction to the Sun in the body-fixed frame is given by the unit vector  , then for a convex body with the uniform surface albedo A

where is the solar constant and R_j(t) is the heliocentric distance in astronomical units. But asteroids are rarely (if ever) convex, so we have to add a visibility function to the final formula

The value of v_j is either 1, if the Sun is visible from the centroid , or 0 if the Sun is occluded by any other triangle of the surface mesh or its zenith distance exceeds  . In our algorithm, the occlusion tests are based on the exact ``ray-triangle stabbing'' queries, optimized by appropriate sorting of triangles projected onto local celestial spheres.

Assuming uniform rotation and a Keplerian orbit, the values of  are computed in nested loops (the rotation phase loop inside the mean anomaly loop). Then the composite Simpson rule is used to compute the mean value. To our surprise, sampling the mean anomaly instead of adhering to the usual routine of uniformly sampling the true anomaly (followed by a correcting factor multiplication) gives a better accuracy when the sampling density is low (in our case 100  100 points).

Unlike Scheeres (2007) or Statler (2009), we use the simple Lambertian emission/diffusion model. Hence, setting A=0 in all computations we concatenate the contributions due to thermal radiation and diffuse reflection. To some degree, our algorithm can be seen as the extension of the old Vokrouhlický & Capek (2002) approach, with most of the efficiency gained by appropriate programming. The models involving up to 200 000 facets were treated on a single CPU (taking at most 2 days of computations), but more challenging cases were run in parallel on about 12 computers (a mesh of 3  10⁶ faces required about a week).

In all computations with Itokawa, we adopted the orbital semi-axis , eccenricity e = 0.280126, the moment of inertia I₃ = 7.77  , and the obliquity . The tests performed indicated that the formal accuracy of our results is about .

Except for explicitly mentioned instances, all the results reported in this paper were computed with the shadowing effects included. No reduction to the centre of mass and principal axes was made (except for the transformation explicitly described in Sect. 3) and all shape models were considered in their original reference systems.

  2.2 Will it converge?

Starting our work, we used the preliminary Gaskell et al. (2006a) topography models of Itokawa downloaded from Abe (2007). This collection consists of four meshes build of quadrilateral patches with additional triangulated variants obtained by simple diagonal cuts of the original facets. Referring to these triangulations, we use the following abbreviations: 6G6 for the mesh of 49 152 triangles, 8G6 for the 196 608 facet set, 10G6 for the 786 432 triangles, and 12G6 for the mesh of 3 145 728 triangular patches. The coding rule we use for the models is the following: kM stands for a model of 768  2^k triangles, belonging to a family M (for example, M is ``G6'' for the aforementioned Gaskell et al. 2006a, collection).

Spanning all obliquities with a step of , we obtained the dotted curves in Fig. 1. All curves have similar shapes, practically flat in the regions of and . This property indicates that the YORP effect on Itokawa is not sensitive to errors in obliquity. To our surprise, Gaskell models with different resolutions result in curves differing mostly by a vertical shift. The vertical distance between 6G6 and 8G6 is larger than, for example, between 10G6 and 12G6, but not to the extent that might suggest approaching some limit with the 12G6 mesh.

Later, we abandoned the G6 family in favor of more recent G8 Hayabusa-based topography models available on the NASA PDS server (Gaskell et al. 2008). The meshes had the same origin and structure as their G6 predecessors. The resulting YORP curves are plotted in Fig. 1 as continuous lines. Once again, the pattern is repeated: all curves are similar, but with a vertical translation depending on the number of facets. Moreover, comparing the distance between G6 and G8 models having the same resolution, we see that the distance increases for more complicated models, so that the results of 12G6 are quite similar to 10G8 for the actual obliquity of Itokawa, . The G8 Gaskell models show even less evidence of approaching a limit, because the vertical distance from 10G8 to 12G8 is almost as large as the one from 6G8 to 8G8 .

Thus we confirm the sensitivity of YORP on Itokawa to the resolution of the Gaskell model reported in Durech et al. (2008b) or Scheeres et al. (2008). Moreover, we find Itokawa highly sensitive to the version of the Gaskell model if we compare the meshes with the same high number of triangular faces. We found no sign of saturation in the YORP effect when approaching 12G8 - the finest available mesh.

The question of shadowing is important for any nonconvex object and ignoring occlusions (i.e. using the zenith distance of the Sun as the only visibility criterion) cannot be a viable approximation unless concavities are very shallow and smooth. But for Itokawa we see many different level shadows: from the large scale body-head occlusions to local shadow lines drawn by small boulders and crater rims. Our computations confirm - in principle - the sensitivity of the computed YORP effect to the occlusions (shadowing) announced by Scheeres et al. (2008,2007). The black dashed line in Fig. 1 was computed according to the 12G8 mesh without the occlusion tests and its shape agrees with the top panel of Fig. 2 of Scheeres et al. (2008). The remaining G8 meshes also give results matching their plots and that are significantly different from the values with the shadowing turned on. However, we find the influence of shadowing in a form different to the bottom panel of the figure from Scheeres et al. (2008). The shadowing effect effect is illustrated in Fig. 2, where we present the differences for each surface triangle, where the pseudo-convex values of were computed without the occlusion tests. The differences reached the maximum of 1.09  and the minimum - 1.78  , but we clipped the data to the range from - 1  (black) to 1  (white) in order to enhance details in Fig. 2. The error introduced by the pseudo-convex model is best visible in the ``head'' part occluded by the ``body'', and in the ``body''-facing slopes of major boulders.

Why do subsequent G6 or G8 models result in simulated YORP curves that differ mostly (although not exclusively) by a vertical translation, hence, by a constant, obliquity-independent torque? There is no doubt that the phenomenon is related to shadowing effects, although one might also suspect the accumulation of discretization or roundoff errors in our software. Such errors should systematically depend on the number of triangles considered. However, a simple test excluded this possibility: we took the 6G8 mesh of 49 152 facets and performed a planar subdivision of each triangle into 4 and 16 subtriangles using edge midpoints as new vertices. In this way we increased the number of triangles without affecting the shape itself. The results for these subdivided meshes agreed with the original 6G8 curve up to  in the worst case of and up to  for the actual obliquity of Itokawa.

If the phenomenon of gradual shifting of curves in Fig. 1 is not a numerical artifact of the YORP simulating algorithm (and after numerous tests we believe it is not), it may be related to a particular algorithm deriving the lower resolution models from the basic 12G8 mesh. To shed some light on this possibility, we decided to create our own downsampled versions of the Itokawa shape using different mesh simplification methods. Unexpectedly, the results raised doubts concerning the use of convergence in a family of different resolution models as an argument for the YORP accuracy.

2.3 Simplified meshes

According to the description accompanying the Gaskell et al. (2008) data files at the NASA PDS repository, lower resolution models were derived from the basic, highest resolution mesh by repeated removal of vertices merging four incident quadrilaterals into one. In terms of triangular variants that we dubbed G8, it also implies that the number of facets drops by factor of 4 after each coarsening operation. Due to the technical limitations, Durech et al. (2008b) did not go beyond 6G6 and 8G6 models, and not much can be deduced from a sequence of two YORP values. So, the 8G6 model was taken as a departure point and 21 meshes were produced with the number of triangular faces decreasing first by 2  10⁴, then by 10⁴, and finally by 10³ (Durech et al. 2008b). A simple algorithm of edge contraction was used, consisting in replacing the shortest edge of the mesh with a single vertex located at the midpoint and repeating the procedure until the requested number of faces is reached (Durech, private communication). The fact that edge contraction for a triangular mesh leads to a different outcome than the merging of quadrilaterals is seen in the results of Durech et al. (2008b): the results for 6G6 and the alternative mesh of 49 152 triangles differ by . The results reached the final   with no sign of approaching a well defined limit.

Comparing the value for the 8G6 Gaskell model from Durech et al. (2008b) with our results shown in Fig. 1 (   ) we found a  difference. But the comparison with the remaining values from Durech et al. (2008b, Table 3) was found to be complicated, because some of the models provided by the author contained a number of nonmanifold faces that had to be removed and patched before use.

Comparing our computations based upon the patched Durech meshes (Fig. 3, red diamonds) with the original data from Durech et al. (2008b) (Fig. 3, magenta triangles), we find differences that are almost systematic. This offset of about 5  is almost independent of the mesh density and reveals systematic errors in the results of Durech et al. (2008b), but it does not contradict the main conclusion, that the sequence of derived meshes YORP values converges to the final 8G6 very slowly; as slowly as the sequence of G8 models approaches 12G8 (Fig. 3, blue circles).

Observing that two meshes with 49 152 triangles derived by simplifying the same 8G6 model using two different methods result in the relative difference of  reaching 0.6  i.e. almost of the computed value, we tried other methods of mesh simplification, but this time taking the best 12G8 model as the primary source. Mesh manipulations were performed using two open-source programs: ReMESH 2.0 (Attene & Falcidieno 2006) and MeshLab 1.2 (Cignoni 2008).

We studied three basic types of simplified meshes:

• E was generated by the principle of the shortest edge contraction, similarly to the Durech meshes. E meshes were obtained directly by the ``Simplify'' filter of the ReMESH package with the option ``priority on the edge length''.

• UR provided an almost uniform meshing with almost equilateral triangles and was generated by the ``Uniform remesh'' filter of ReMESH. The interpolation in this algorithm typically produces noisy patterns that we removed by adjusting the uniform mesh to the original 12G8 model by the RIMLS method (Robust Implicit Moving Least Squares projection of Öztireli et al. 2009, implemented in MeshLab).

• Q was substantially different from the previous ones: decisions about the choice of an edge to be collapsed were made not only according to its length, but also with a goal to preserve sharp features of the surface for as long as possible. The classical Garland & Heckbert (1997) method of quadric error metrics implemented in MeshLab was used for this purpose.

Each mesh family consisted of 11 members. In contrast to Durech et al. (2008b), we assumed a logarithmic spacing with the 768  2^k rule. Thus, starting from the 11E mesh of 1 572 864 triangles, we stepped down to 1E with 1536 faces and similarly so for  and .

The results of the YORP computations are presented in Fig. 4. Compared with Fig. 3 it shows that our E family converges to Gaskell 12G8 roughly the same way as the Durech meshes did to 8G6, although more smoothly and slightly more horizontally in the 9E-11E zone. The YORP values based upon the uniform triangulations UR show a more nonlinear trend on the semi-logarithmic plot and their convergence is better than for E, but neither of these two approach that of the Q family. The differences between 12G8 and  remain below 5  until k=3. The 3Q mesh of 6144 triangles still gives the YORP closer to the 12G8 based reference value than Gaskell's 10G8 with its 786 432 facets.

Figure 5 presents the dependence of the YORP effect on obliquity for models with the same number of triangular faces (49 152). It confirms our statement that vertical translation of the YORP curves is a shadowing-based phenomenon. It does not depend on the number of triangles, but rather on the presence or absence of sharp features of the surface. It is also an example of how misleading the pseudo-convex approximation can be.

The superiority of the Garland-Heckbert decimation for the YORP simulation deserves discussion. In order to understand why it works so well for Itokawa, recall:

1.

With the obliquity of Itokawa close , most of the YORP torque should be generated in the equatorial regions of the asteroid.

2.

Two major smooth regions that exist on Itokawa are located close to the poles (Barnouin-Jha et al. 2008; Demura et al. 2006).

3.

The Garland-Heckbert method spares small triangles at sharp features, creating large triangles in smooth areas (Garland & Heckbert 1997).

For these reasons Itokawa is particularly well suited for the application of the Garland-Heckbert method and even a mesh with a small number of triangles may provide good YORP values. Figure 6 clearly illustrates the superiority of the Garland-Heckbert decimation and its feature-preserving properties. The 6Q mesh of 49 152 triangles looks almost like the original 12G8, whereas the 6G8 of Gaskell et al. (2008), using a different mesh simplification method, presents a much smoother surface in spite of having the same number of facets as 6Q.

The notion of ``mesh resolution'' is justified to some extent for an almost uniform UR mesh, but it would be misleading in the context of Q models which are based upon the Garland-Heckbert decimation. In terms of the edge lengths of a given model, the meshes with the same number of facets gives quite different triples (min, max, mean) expressed in meters. Taking the 49 152 facets as an example, we obtain (1.6, 15, 4.7) for 6G8, (1, 13, 5.2) for 6E, (2.6, 6.8, 4.4) for 6UR, and (0.8, 22, 4.6) for 6Q. Judging by the number of triangles or by the mean value of edge lengths, one would qualify all models as ``4-5 m resolution'' ones. Yet their influence on the computed YORP values is quite different.

The rapid convergence of Q meshes to 12G8 or relatively good behaviour of UR leads to a question: can we use the rapid convergence of a sequence of simplified meshes to 12G8 as an argument that the model is sufficiently accurate for the YORP computations?

In our opinion, the answer to this question is negative according to an argument from another family of meshes: the 8G6 mesh was taken as the reference model and we reduced its number of triangles using the Garland-Hecbkert edge contraction method. Figure 7 compares our ad hoc models with the results based on Durech meshes. Looking at the almost flat curve approaching the 8G6 based YORP value and ignoring all the remaining points in Fig. 7 one might conclude that the Gaskell mesh 8G6 with its 196 806 facets provides an ultimate YORP value and there is no need for higher resolution models.

Obviously, this conclusion is false. But we can reject it only because we know the results for 10G6 and 12G6. So, even if Scheeres et al. (2007) or Durech et al. (2008b) observed a good convergence to the 8G6 values, this would not be evidence of the YORP value accuracy, but only a proof that lower resolution meshes had been appropriately sampled.

2.4 More than Itokawa

Is the sensitivity to surface fine details a special property of Itokawa, or is it a more common problem? One cannot rule out the possibility that only a special combination of a large scale shape and small scale roughness impedes the YORP modeling. Although Statler (2009) addressed a similar question, but he used relatively simple models of synthetic objects. Our attention was drawn to 433 Eros for two reasons. First, Gaskell (2008) provides space-based shape models of Eros with the same number of facets as for Itokawa. On the other hand, the topography features of Eros make it the negative of Itokawa: the former is mostly cratered, whereas the latter is mostly bouldered. Assuming an obliquity of  , the maximum moment of inertia I₃ = 4.9781  , and the orbital elements , e = 0.223, we found the YORP values of  ranging from -6.46  for the 49 152 triangles of 6GE, to -9.47  for the complete 12GE mesh of Gaskell (2008). The relative difference is smaller than for Itokawa, but the convergence is equally poor, as shown in Fig. 8.

Curiously, the last value, obtained from over 3  10⁶ facets, is quite close to what we computed according to a much older 7790 triangle shape model by Thomas et al. (2002): -1.07  . This time we are in good agreement with the -1.12  computed by Capek & Vokrouhlický (2004) using the same model, as well as with the observational estimate of -1.0  suggested by Durech (2005).

  3 Centre of mass shift versus the spin axis tilt

Scheeres & Gaskell (2008) were the first to observe that the strength of the mean YORP torque on Itokawa can significantly decrease if the centre of mass (COM) is shifted with respect to the present origin of the topographic model reference frame. In their approach, the insolation and temperature distribution (hence the forces and reference torques on each facet) were computed once, assuming the ``centre of figure'' origin and rotation around the original model Oz axis. Then, corrections to the torques were computed to account for different locations of the current rotation axis, say O'z', passing through a new origin in the Oxy plane and parallel to Oz. The corrections were linearized according to the principle that the translation of COM is a small fraction of the body size. The possibility of skewing O'z' with respect to Oz was discarded by Scheeres & Gaskell (2008) as a small quantity of the second order. But the most important simplification introduced was to neglect the shadowing, i.e. using a pseudo-convex approximation. It was not stated explicitly in the article, but may be deduced from the comparison of Scheeres et al. (2008, Fig. 2) and Scheeres & Gaskell (2008, Fig. 1).

The sensitivity of YORP on Itokawa to the shape models, so clearly manifested in Sect. 2, incited us to check the influence of the Scheeres & Gaskell (2008) simplifications on the computed  . To this end we first sampled various locations of O' in the body frame, maintaining the rotation axis O'z' parallel to Oz. Using the 6QR model (a Garland-Heckbert decimated mesh of 49 152 facets, derived from 12G8), we were able to scan a wide range of presumed centres of mass, each time performing a full simulation with obliquity . For O' located at in the Oxyz frame, we subtracted and  from the coordinates of the mesh vertices. Then we computed the mean YORP torque projection on Oz' assuming uniform rotation around this axis. The computations were performed as described in Sect. 2.1, without any additional simplifications, and repeated independently for each and . Interested in the location of  points on the  plane (with the linear interpolation on the grid), we did not recompute the moments of inertia for a new COM, because their values are irrelevant when the numerator in Eq. (5) is zero. The axes Ox' and Oy' were maintained parallel to the original Ox and Oy, because their rotation around Oz' has no effect after averaging with respect to the daily rotation.

The obtained results confirm those of Scheeres & Gaskell (2008) qualitatively, but differ in numerical values. While Scheeres and Gaskell found the ``zero deceleration line'' as a straight line passing through the points and , our simulation locates the points at and . Thus our line of possible zero YORP locations of COM has a slightly different direction and - what is more important - passes further from the origin, requiring a stronger inhomogeneity to justify the COM translation. The minimum distance OO' according to our results is about , compared to  obtained by Scheeres & Gaskell (2008) - the difference that should be attributed indirectly to the influence of shadowing. Directly, the greater OO' is a consequence of three the times stronger YORP effect to be suppressed as compared to the pseudo-convex approximation of Scheeres & Gaskell (2008) (see the green lines in Fig. 5 - without shadowing, the effect is three times weaker).

Recalling that the Ox axis points to the ``head'' of Itokawa, the translation in the positive x direction agrees with a hypothesis that the ``head'' part has a higher density compared to the ``body'' part (Scheeres & Gaskell 2008). Another hypothesis of Scheeres & Gaskell (2008), that the ``neck'' part has a higher density, is also possible, although it requires a greater inhomogeneity.

Investigating the parallel translation of COM, we identified the line of the steepest descent of on the  plane

the line passing through the origin , and perpendicular to the line of . Our next step consisted of scanning the line (6) in the direction of increasing and . For each new COM, located at the distance

from the original centre, we computed at various orientations of the spin axis Oz' with respect to Oz, applying a sequence of two 3-1 rotations: first by angle  around Oz, and then by the ``tilt angle''  around the new Ox axis. More precisely, we took the original mesh, translated its vertices to the O' centered frame by subtracting and , and then rotated the vertices to the final O'x'y'z' frame using the angles  and . These angles are, respectively, the longitude and colatitude of the Oz' pole in the Oxyz frame. The moments of inertia were not recomputed, because: a) assuming uniform rotation around Oz', we implicitly postulate that the third row and column of the inertia matrix are free from off-diagonal terms; b) the value of I₃ is irrelevant for the  condition; and c) contribution of remaining non-diagonal terms disappears after averaging with respect to the daily revolution.

So, for each pair of a grid of was filled with the associated values of computed from the full simulation with , as usual. The line of was identified and the smallest value of colatitude  was registered together with its associated . Thus we collected a set of quadruples with the property that for a given COM distance R they produce the null YORP effect by a minimum deviation of the rotation pole from the original Oz direction. Figure 9 presents the minimum values of found at the COM located on line (6) at a given distance from the Gaskell's model centre. The points location is accurate up to about due to the interpolation error. Both angles, expressed in degrees, can be roughly approximated as

(solid lines in Fig. 9), and

Scheeres & Gaskell (2008) neglected the effect of axis rotation; their reasoning was fairly justified on the ground of transforming a torque while keeping the force on each facet unchanged. However, the rotation of the spin axis changes the mean energy acquired and re-radiated by different points on an asteroid. Thus, it relocates polar circles and the equator on a highly irregular object like Itokawa. Interestingly, the sensitivity to the spin axis orientation in the body frame (which should not be confused with a weak sensitivity to the obliquity at  close to  ) can be noticed even in the pseudo-convex approximation, being related to the global shape irregularity.

4 Conclusions

The main conclusion is that the value of the YORP effect on Itokawa within the Rubincam's zero conductivity approximation strongly depends on the adopted surface model mesh. In this respect we have confirmed the general, qualitative conclusions of Scheeres et al. (2007) or Durech et al. (2008b), although we differ in quantitative terms. Our simulations were pushed up to the best available mesh resolution, showing that the sensitivity does not decrease when a few centimeter level is reached. Moreover, the YORP value at the highest resolution 12G8 (-5.53  ) is the farthest from the observational constraint provided by Durech et al. (2008b). In the agreement with Statler (2009), we find the YORP surface detail sensitivity is a generic phenomenon that exists also for Eros or for artificial Gaussian random spheres with additional surface features.

As far as the sensitivity of the Itokawa YORP to the inhomogeneities of mass distribution is concerned, we used a more accurate modelling approach than Scheeres & Gaskell (2008). According to our results, the YORP-cancelling parallel translation of axes should be at least twice as large as the one found by Scheeres & Gaskell (2008). This might suggest a weaker dependence on the mass distribution, but combining the translation with the rotation of axes we restore the high sensitivity.

In our opinion, the implications of the attempts to model the YORP effect on Itokawa are far reaching, and the present paper adds more weight to pessimistic conclusions from the papers of Scheeres et al. (2007), Durech et al. (2008b), Scheeres & Gaskell (2008), and Statler (2009). After the lesson of Itokawa, the state of art in YORP modelling has to be judged according to an answer to a fundamental question: What are the error bounds of the simulated YORP values? The answer can no longer be limited to the numerical properties of an applied algorithm.

If the assumptions of the present modelling tools are justified, we cannot trust even the sign of the YORP effect computed from low resolution models. And once we obtain a detailed topography model, like a mega-faceted Itokawa or Eros mesh, we have no way to decide whether its resolution is fine enough. Exploring the sequence of coarsened models, it is difficult to obtain a satisfactory answer. As we tried to demonstrate in the present paper, this kind of test tells us much about the quality of the simplification algorithm, and nothing about the quality of the departure mesh as a source of the YORP effect.

In the worst case, if the sensitivity of YORP to tiny features on a surface is a physical fact, simulations of this effect will have to be qualified similarly to those of chaotic motion: uncertain, although not meaningless. But it is not evident as yet whether the hypersensitivity is a physical phenomenon or a product of incorrect assumptions (c.f. the remarks of Scheeres et al. (2007) concerning the secondary illumination). Recalling that isothermal bodies reveal no YORP torque regardless of their shapes (Breiter & Michalska 2008), we should observe a key role of the temperature surface gradient. All classical assumptions, like zero conductivity, a plane-parallel model, or ignoring indirect lighting, tend to increase the gradient. It looks as if all phenomena that make the surface temperature distribution more even and lower the gradient are neglected in current YORP models. The assumptions were justified for large, locally smooth and convex bodies, but there are no reasons to expect that they are suitable for irregularly shaped asteroids with craters and boulders. Regardless of computational cost, thermal models for the YORP effect have to be improved until a realistic interplay of conduction and radiation heat exchange is reached.

Acknowledgements

The authors appreciate the help and advice of Bartosz Oczujda concerning computer graphics rendering issues. We thank Dr. Robert Gaskell for information about details of Itokawa models. The work was supported by the Polish Ministry of Science and Higher Education grant N N203 302535. The work of D.V. was supported by the Czech Grant Agency (grant 205/08/0064) and the Research Program MSM0021620860 of the Czech Ministry of Education. The MeshLab tool used in our research has been supported by the European Networks of Excellence EPOCH and (partially) by Aim@Shape.

References

• Abe, M. 2007, Hayabusa Project Science Data Archive, http://hayabusa.sci.isas.jaxa.jp/shape.pl
• Attene, M., & Falcidieno, B. 2006, in Proceedings of Shape Modeling International (SMI'06) (IEEE C.S. Press), 271
• Barnouin-Jha, O. S., Cheng, A. F., Mukai, T., et al. 2008, Icarus, 198, 108 [CrossRef] [NASA ADS]
• Bottke, Jr., W. F., Vokrouhlický, D., Rubincam, D. P., & Nesvorný, D. 2006, Annu. Rev. Earth Planet. Sci., 34, 157 [CrossRef] [NASA ADS]
• Breiter, S., & Michalska, H. 2008, MNRAS, 388, 927 [CrossRef] [NASA ADS]
• Burns, J. A., & Safronov, V. S. 1973, MNRAS, 165, 403 [NASA ADS]
• Capek, D., & Vokrouhlický, D. 2004, Icarus, 172, 526 [CrossRef] [NASA ADS]
• Cignoni, P. 2008, MeshLab, http://meshlab.sourceforge.net/
• Demura, H., Kobayashi, S., Nemoto, E., et al. 2006, Science, 312, 1347 [CrossRef] [NASA ADS]
• Durech, J. 2005, A&A, 431, 381 [EDP Sciences] [CrossRef] [NASA ADS]
• Durech, J., Vokrouhlický, D., Kaasalainen, M., et al. 2008a, A&A, 489, L25 [EDP Sciences] [CrossRef] [NASA ADS]
• Durech, J., Vokrouhlický, D., Kaasalainen, M., et al. 2008b, A&A, 488, 345 [EDP Sciences] [CrossRef] [NASA ADS]
• Garland, M., & Heckbert, P. S. 1997, in SIGGRAPH 97 Conference Proceedings, Aug., ed. T. Whitted, Annu. Conf. Ser. ACM SIGGRAPH (Addison Wesley), 209
• Gaskell, R., Barnouin-Jha, O., Scheeres, D., et al. 2006a, in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 21-24 August, Keystone, Colorado
• Gaskell, R., Saito, J., Ishiguro, M., et al. 2006b, in Lunar and Planetary Inst. Technical Report, 37th Annu. Lunar Planet. Sci. Conf., ed. S. Mackwell & E. Stansbery, 37, 1876
• Gaskell, R., Saito, J., Ishiguro, M., et al. 2008, Gaskell Itokawa Shape Model V1.0, HAY-A-AMICA-5-ITOKAWASHAPE-V1.0, NASA Planetary Data System
• Gaskell, R. W. 2008, Gaskell Eros Shape Model V1.0, NEAR-A-MSI-5-EROSSHAPE-V1.0, NASA Planetary Data System
• Kaasalainen, M., Kwiatkowski, T., Abe, M., et al. 2003, A&A, 405, L29 [EDP Sciences] [CrossRef] [NASA ADS]
• Kaasalainen, M., Durech, J., Warner, B. D., Krugly, Y. N., & Gaftonyuk, N. M. 2007, Nature, 446, 420 [CrossRef] [NASA ADS]
• Kitazato, K., Abe, M., Ishiguro, M., & Ip, W.-H. 2007, A&A, 472, L5 [EDP Sciences] [CrossRef] [NASA ADS]
• Lowry, S. C., Fitzsimmons, A., Pravec, P., et al. 2007, Science, 316, 272 [CrossRef] [NASA ADS]
• Mysen, E. 2008, A&A, 484, 563 [EDP Sciences] [CrossRef] [NASA ADS]
• Nesvorný, D., & Vokrouhlický, D. 2008, AJ, 136, 291 [CrossRef] [NASA ADS]
• Ostro, S. J., Benner, L. A. M., Nolan, M. C., et al. 2004, Meteoritics and Planetary Science, 39, 407 [NASA ADS]
• Ostro, S. J., Benner, L. A. M., Magri, C., et al. 2005, Meteoritics and Planetary Science, 40, 1563 [NASA ADS]
• Öztireli, A. C., Guennebaud, G., & Gross, M. 2009, EUROGRAPHICS 2009, 28, accepted
• Rubincam, D. P. 2000, Icarus, 148, 2 [CrossRef] [NASA ADS]
• Scheeres, D. J. 2007, Icarus, 188, 430 [CrossRef] [NASA ADS]
• Scheeres, D. J., & Gaskell, R. W. 2008, Icarus, 198, 125 [CrossRef] [NASA ADS]
• Scheeres, D. J., & Mirrahimi, S. 2008, Celest. Mech. Dyn. Astron., 101, 69 [CrossRef] [NASA ADS]
• Scheeres, D. J., Ostro, S. J., Werner, R. A., Asphaug, E., & Hudson, R. S. 2000, Icarus, 147, 106 [CrossRef] [NASA ADS]
• Scheeres, D. J., Abe, M., Yoshikawa, M., et al. 2007, Icarus, 188, 425 [CrossRef] [NASA ADS]
• Scheeres, D. J., Mirrahimi, S., & Gaskell, R. W. 2008, in Lunar and Planetary Inst. Technical Report, Lunar and Planetary Institute Science Conference Abstracts, 39, 2348
• Statler, T. S. 2009, Icarus, 202, 502 [CrossRef] [NASA ADS]
• Taylor, P. A., Margot, J.-L., Vokrouhlický, D., et al. 2007, Science, 316, 274 [CrossRef] [NASA ADS]
• Thomas, P. C., Joseph, J., Carcich, B., et al. 2002, Icarus, 155, 18 [CrossRef] [NASA ADS]
• Vokrouhlický, D., & Capek, D. 2002, Icarus, 159, 449 [CrossRef] [NASA ADS]
• Vokrouhlický, D., Capek, D., Kaasalainen, M., & Ostro, S. J. 2004, A&A, 414, L21 [EDP Sciences] [CrossRef] [NASA ADS]

Footnotes

... strategies

All simplified shape models derived for the purpose of the present study are available from the first author upon request in the form of Wavefront OBJ data files.

... model

Although Fig. 2 of Capek & Vokrouhlický (2004), reproduced also as Fig. 6 in Bottke et al. (2006), might suggest a weak dependence of the spin YORP on conductivity for 6489 Golevka, this effect may merely reflect numerical errors. Also, the spin related part of these figures has an incorrect shape due to the interpolation of undersampled data. For a more reliable version, see Fig. 9 in Vokrouhlický & Capek (2002).

... derived

The term ``dumbed down'' used by Gaskell et al. (2008) is more descriptive.

... weaker)

We verified that without shadowing our results agree with Scheeres & Gaskell (2008) exactly in and with a 2 m difference in .

Copyright ESO 2009

All Figures

---

Copyright ESO 2009