A&A 507, 10731081 (2009)
The YORP effect on 25 143 Itokawa
S. Breiter^{1}  P. Bartczak^{1}  M. Czekaj^{2,1}  B. Oczujda^{1}  D. Vokrouhlický^{3}
1  Astronomical Observatory, Adam Mickiewicz University, Soneczna 36, 60286 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 YarkovskyO'KeefeRadzievskiiPaddack (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 bodyfixed 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 YarkovskyO'KeefeRadzievskiiPaddack (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 semiaxis 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 YORPinduced 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 ``pseudoconvex'' 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 rubblepile 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 megafacet 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 StefanBoltzmann constant, and I_{3} 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 reradiated power flux and of the heat convected into an asteroid's interior; the righthand 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 semiinfinite 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 bodyfixed 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 ``raytriangle 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^{6} faces required about a week).
In all computations with Itokawa, we adopted the orbital semiaxis , eccenricity e = 0.280126, the moment of inertia I_{3} = 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).
Figure 1: YORP effect for Gaskell models of Itokawa. Dots ( from top to bottom): 6G6 (magenta), 8G6 (red), 10G6 (green), 12G6 (blue). Lines ( from top to bottom, same colours as for G6): 6G8, 8G8, 10G8, 12G8. Dashed line: pseudoconvex 12G6. 

Open with DEXTER 
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 Hayabusabased 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.
Figure 2: Differences between the complete shadowing and pseudoconvex YORP values for the Itokawa 8G8 model and obliquity . 

Open with DEXTER 
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 bodyhead 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 pseudoconvex 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 pseudoconvex 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, obliquityindependent 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^{4}, then by 10^{4}, and finally by 10^{3} (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.
Figure 3: YORP effect for Itokawa Gaskell models G6 (green squares), G8 (blue circles) and Durech meshes. Triangles represent the data quoted from Table 3 of Durech et al. (2008b); red diamonds label the results of our computation using patched Durech meshes. 

Open with DEXTER 
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 opensource 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.
Figure 4: YORP effect for Itokawa Gaskell models G8 (blue circles) and derived meshes: E (magenta triangles), UR (red diamonds), Q (green squares). 

Open with DEXTER 
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 9E11E zone. The YORP values based upon the uniform triangulations UR show a more nonlinear trend on the semilogarithmic 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: YORP effect as a function of obliquity for selected Itokawa models: 12G8 (blue dots), 6E (magenta, dashed), 6UR (red, dashdot), and 6Q (green, solid). Left  pseudoconvex approximation, right  with shadowing effects. The remaining part of each curve is symmetric with respect to the line. 

Open with DEXTER 
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 shadowingbased 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 pseudoconvex approximation can be.
Figure 6: Itokawa shape according to three models: 12G8 (3 145 728 triangles)  top, 6Q (49 152 triangles)  middle, and 6G8 (49 152 triangles)  bottom. 

Open with DEXTER 
The superiority of the GarlandHeckbert 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 (BarnouinJha et al. 2008; Demura et al. 2006).
 3.
 The GarlandHeckbert method spares small triangles at sharp features, creating large triangles in smooth areas (Garland & Heckbert 1997).
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 GarlandHeckbert 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 ``45 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 GarlandHecbkert 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.
Figure 7: YORP effect for Itokawa Gaskell models G6 (blue circles) and two meshes derived by simplification of 8G6: Durech set (green squares) and GarlandHeckbert (red diamonds). 

Open with DEXTER 
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 spacebased 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_{3} = 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.
Figure 8: YORP effect for Eros Gaskell (2008) models GE (blue circles) and derived meshes: E (green squares) and Q (red diamonds). 

Open with DEXTER 
Curiously, the last value, obtained from over 3 10^{6} 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 pseudoconvex 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 GarlandHeckbert 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 pseudoconvex 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 31 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 offdiagonal terms; b) the value of I_{3} is irrelevant for the condition; and c) contribution of remaining nondiagonal 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 reradiated 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 pseudoconvex approximation, being related to the global shape irregularity.
Figure 9: Null YORP combinations of centre of mass (COM) shift R and Oz axis tilt for Itokawa. 

Open with DEXTER 
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 YORPcancelling 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 megafaceted 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 planeparallel 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.
AcknowledgementsThe 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
 BarnouinJha, 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., BarnouinJha, O., Scheeres, D., et al. 2006a, in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2124 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, HAYAAMICA5ITOKAWASHAPEV1.0, NASA Planetary Data System
 Gaskell, R. W. 2008, Gaskell Eros Shape Model V1.0, NEARAMSI5EROSSHAPEV1.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 .
All Figures
Figure 1: YORP effect for Gaskell models of Itokawa. Dots ( from top to bottom): 6G6 (magenta), 8G6 (red), 10G6 (green), 12G6 (blue). Lines ( from top to bottom, same colours as for G6): 6G8, 8G8, 10G8, 12G8. Dashed line: pseudoconvex 12G6. 

Open with DEXTER  
In the text 
Figure 2: Differences between the complete shadowing and pseudoconvex YORP values for the Itokawa 8G8 model and obliquity . 

Open with DEXTER  
In the text 
Figure 3: YORP effect for Itokawa Gaskell models G6 (green squares), G8 (blue circles) and Durech meshes. Triangles represent the data quoted from Table 3 of Durech et al. (2008b); red diamonds label the results of our computation using patched Durech meshes. 

Open with DEXTER  
In the text 
Figure 4: YORP effect for Itokawa Gaskell models G8 (blue circles) and derived meshes: E (magenta triangles), UR (red diamonds), Q (green squares). 

Open with DEXTER  
In the text 
Figure 5: YORP effect as a function of obliquity for selected Itokawa models: 12G8 (blue dots), 6E (magenta, dashed), 6UR (red, dashdot), and 6Q (green, solid). Left  pseudoconvex approximation, right  with shadowing effects. The remaining part of each curve is symmetric with respect to the line. 

Open with DEXTER  
In the text 
Figure 6: Itokawa shape according to three models: 12G8 (3 145 728 triangles)  top, 6Q (49 152 triangles)  middle, and 6G8 (49 152 triangles)  bottom. 

Open with DEXTER  
In the text 
Figure 7: YORP effect for Itokawa Gaskell models G6 (blue circles) and two meshes derived by simplification of 8G6: Durech set (green squares) and GarlandHeckbert (red diamonds). 

Open with DEXTER  
In the text 
Figure 8: YORP effect for Eros Gaskell (2008) models GE (blue circles) and derived meshes: E (green squares) and Q (red diamonds). 

Open with DEXTER  
In the text 
Figure 9: Null YORP combinations of centre of mass (COM) shift R and Oz axis tilt for Itokawa. 

Open with DEXTER  
In the text 
Copyright ESO 2009