A&A 376, 302-309 (2001)
DOI: 10.1051/0004-6361:20010935

Interpretation of lightcurves of precessing asteroids

M. Kaasalainen

Observatory, PO Box 14, 00014 University of Helsinki, Finland

Received 23 April 2001 / Accepted 26 June 2001

Abstract
I show that the lightcurves of freely precessing asteroids can be analyzed using basically the same methods as with objects in relaxed rotation. At least two different observation geometries (with respect to the angular momentum vector) are needed for a unique full solution; an informative basic model can be constructed at only one geometry, if the observations span several precession cycles. Even noisy data yield a good estimate of the dynamical parameters, whereas accurate data are required for a good shape solution.

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

1 Introduction

Most asteroid lightcurves are single-periodic and thus imply that the rotational state of the object is relaxed, i.e., that the asteroid rotates about the axis corresponding to its maximum moment of inertia, and this axis is aligned with the angular momentum vector. Several lightcurves, however, exhibit double-periodic behaviour. In most of these cases the object has been shown to be a binary asteroid either by lightcurve analysis (Pravec & Hahn 1997; Mottola & Lahulla 2000) or by direct imaging. So far, few asteroids have been suspected or found to be in an excited state of rotation, i.e., precessing. This is a direct consequence of the fact that no real object is completely rigid: tumbling motion causes the asteroid to lose rotational energy and thus reach the relaxed state in a relatively short time interval (Harris 1994). Indeed, to date, the only complete analysis of a tumbling asteroid is that of 4179 Toutatis based on radar images (Hudson & Ostro 1998). However, collisions and other fast perturbative events ensure that the supply of precessing asteroids is never exhausted. With the growing number of observations, lightcurves of potential precessing targets are certain to be measured every now and then. Small Near-Earth Objects are natural candidates for such observations: as the body size decreases, relaxation time increases while the time intervals between significant collisions decrease. The major bonus in observing a precessing asteroid is that only a short observation span is needed to construct a model since, due to the complicated motion, the dynamical parameters leave strong fingerprints on the lightcurve, and all parts of the target are usually seen and illuminated within a short period of time.

The lightcurves of asteroids in relaxed rotation contain a wealth of information on the objects' shapes, scattering properties, and rotational states, well extractable with suitable optimization techniques (Kaasalainen et al. 2001a,b). In principle, precessing objects are no exception: the parameter space only contains a few more dynamical parameters. However, they have some idiosyncratic properties, so care must be taken in extending the scheme for principal-axis rotators to include excited states of rotation. I describe the general aspects of the dynamical parameters in Sect. 2 (a detailed account of some necessary concepts is given in the Appendix). The inversion procedure is discussed in Sect. 3, and examples are presented in Sect. 4. Section 5 presents a summary.

2 Dynamical parameters

A total of eight parameters are needed for a complete description of the motion of the inertia ellipsoid in force-free precession. The body shape is described in the coordinate frame defined by the inertia ellipsoid, so the shape/scattering parameter part of the problem is exactly the same as in relaxed rotation (including the possible straightforward use of the computed inertia tensor in a regularizing function). Many eight-parameter sets are possible; I use the set ${\cal D}=(\vec{L},\phi_0,\theta_0,\psi_0,I_1,I_2)$, where $\vec{L}$ denotes the constant angular momentum vector (its size given by L and its direction by the ecliptic latitude $\beta$ and longitude $\lambda$), while $\phi_0,\theta_0,\psi_0$ are the standard Euler angles of, respectively, precession, nutation (or tilt), and rotation at some epoch t₀. I give these three angles in the so-called x-convention, used in, e.g., Goldstein (1980), Samarasinha & A'Hearn (1991), and Kryszczynska et al. (1999).

	Figure 1: Synthetic lightcurve of an asteroid in short-axis mode.
Open with DEXTER

Hudson & Ostro (1998) also use the same convention, but their angles are different from the standard set, describing the orientation in a given inertial frame rather than precession around the vector of angular momentum (the geometric transformation between the two is straightforward). The last two parameters $I_1,\ I_2$ of $\cal D$ are the lengths of the inertia ellipsoid axes perpendicular to the one around which the ellipsoid can be seen as rotating (this third axis is always normalized to unity). When $I_1<1,\ I_2<1$, the motion can thus be called short-axis mode (SAM), and when $I_1>1,\ I_2>1$, the motion is in long-axis mode (LAM). Note that these modes are mostly geometric descriptions - from the dynamical point of view, the ellipsoid always rotates about one extremal axis and this axis precesses about $\vec{L}$ (the case of I_i on different sides of unity is redundant; it can always be described as a SAM or a LAM and its limiting principal-axis state is unstable). Whether the ellipsoid is prolate or oblate with respect to this axis is mostly inconsequential. Cigar-shaped objects are more likely to be in LAM rather than SAM, since an axisymmetric elongated shape has only long-axis modes (except for one short-axis mode when it is in exact principal-axis rotation about the short axis). The reverse holds for flattened, disk-shaped objects.

The time evolution of the ellipsoid's orientation in space is given by the equations of motion for $\phi$, $\theta$, and $\psi$ (see the Appendix for this and other details of the concepts discussed in the main text). If $\vec{v}_{\rm ecl}$ denotes a vector in the ecliptic coordinate system, the same vector seen from the asteroid's coordinate system, $\vec{v}_{\rm ast}$, is

$$\vec{v}_{\rm ast}=R_z(\psi)\,R_x(\theta)\,R_z(\phi)\, R_y(\tilde\beta)\,R_z(\lambda)\vec{v}_{\rm ecl},$$ (1)

where the polar angle $\tilde\beta=90^{\circ}-\beta$, and $R_i(\alpha)$ is the rotation matrix corresponding to a rotation through angle $\alpha$ about the current i-axis in the positive direction. If a given $\cal D$ turns out to describe rotation about an extremal axis other than the one set to unity, the correct coordinate system is obtained by interchanging and renormalizing the extremal axes. In this way one can always use the same equations for both LAM and SAM, which simplifies the analysis considerably - the other option of staying in a fixed coordinate system would result in different equations and types of motion for the two modes (as in, e.g., Samarasinha & A'Hearn 1991).

The tumbling motion is governed by two periods: one, $P_{\psi}$, for the rotation about the extremal axis, and the other, $P_{\phi}$, for the precession about $\vec{L}$. The nutation period is exactly one half of the rotation period. While $P_{\psi}$ is a real period and a constant of the motion, $P_{\phi}$ is not an actual period but rather defined through another constant  $\Delta\phi$:

$$P_{\phi}= 2\pi {P_{\psi}\over \Delta\phi}\cdot$$ (2)

The precession angle $\phi$ always changes by the amount $\Delta\phi$ during the time interval $P_{\psi}$ (measured from any epoch). Unless $P_{\psi}$ and $P_{\phi}$ are commensurable, $P_{\phi}$ can only be interpreted as an average period of $\phi$.

The two characteristic periods are just as important as the single period in relaxed rotation, since they define the power spectrum of a lightcurve. Constraining these periods is of great value in constraining the parameter space of the inverse problem. The (main) peaks of the spectrum are located at frequencies that are linear combinations of $f_{\psi}$ and $f_{\phi}$, where f=1/P. Prominent peaks are usually found at $2f_{\phi}$ and $2(f_{\phi}\pm f_{\psi})$, where the plus sign holds for LAMs and the minus sign for SAMs; other low harmonics and combinations are typically seen as well. The factor two is the same shape-related phenomenon as in double-sinusoidal ordinary lightcurves. The combination of the frequencies is perhaps best understood when one considers the apparent frequency in the limiting case of principal-axis rotation (i.e., when the nutation angle is $\theta=0$): the two frequencies are still present, but they are seen only as the single frequency $f_{\phi}\pm f_{\psi}$. When the body is slightly tilted, an additional spectrum peak corresponding to the precession frequency will make its appearance. As an example, Fig. 1 shows a simulated lightcurve of an object (the same as in Fig. 6) in SAM (the nutation angle wobbles between $10^\circ$ and $40^\circ$), and the power spectrum is plotted in Fig. 2.

	Figure 2: The power spectrum of the lightcurve of Fig. 1.
Open with DEXTER

The scattering law is a combination of Lommel-Seeliger and Lambert laws such as in Kaasalainen & Torppa (2001). From left to right, the visible peaks correspond to $2f_{\phi}+2f_{\psi}$, $2f_{\phi}$, $2f_{\phi}-2f_{\psi}$, and $f_{\phi}$, with $P_{\phi}=4.70,\ P_{\psi}=19.9$. Though not necessary, a Fourier series representation is sometimes useful for checking the quality and the underlying patterns of the data. In general, the brightness B of any bona fide doubly periodic lightcurve can be expanded as

$$B=\sum_{ij} a_{ij}\cos(i\omega_1 t+j\omega_2 t)+ b_{ij}\sin(i\omega_1 t+j\omega_2 t),$$ (3)

where $\omega_i=2\pi/P_i$, i=0,1,2,..., and j=...,-1,0,1,... (only positive integers when i=0). Note that several i,j-combinations in coefficients are present in principle even for low orders. Thus the sometimes used approximative expansion with two one-dimensional Fourier series (Pravec 1997) is more suitable for binaries than for complex rotators: the lightcurves of the former have an additive two-component physical origin. Also, the two most prominent peaks in the power spectrum are not directly the two characteristic frequencies. Strictly speaking, a proper Fourier series can only be obtained if the observing geometry is fixed with respect to the angular momentum vector $\vec{L}$ so that the observed brightness is truly cyclic in both $\phi$ and $\psi$. If the geometry changes, the expansion may still be valid, at least if the rate of change is slow compared to the tumbling motion, and the compensation due to the changing solar phase angle can be well approximated.

3 Inverse problem

A convenient parameter set for the inverse problem is ${\cal D}_{\rm inv}=(\tilde\beta,\lambda,\phi_0,\psi_0,P_{\psi},P_{\phi}, I_{\rm s},I_{\rm a})$, where the symmetric $I_{\rm s}=(I_1+I_2)/2$ and the antisymmetric $I_{\rm a}=(I_1-I_2)/2$. Using the search region $I_{\rm s}>0.5$, $\vert I_{\rm a}\vert<0.5$ automatically takes into account the triangle inequalities fulfilled by the moments of inertia of all real bodies. This diagonal strip in the (I₁,I₂)-plane is further reduced by using only one sign for $I_{\rm a}$ (the other sign describing the same situation with indices 1 and 2 interchanged) and dropping out the redundant region where I₁ and I₂ are on different sides of unity value. The resulting sections - one for LAMs and one for SAMs - in which I₁ is always an extremal axis and I₂ the intermediate one are shown in Fig. 3.

	Figure 3: The boundaries of the allowed regions in $I_1,\ I_2$-plane.
Open with DEXTER

Also shown are the hyperbolae bounding the regions where it is possible to have a given ratio $P_{\psi}/P_{\phi}$. These hyperbolae (on which rotation is principal-axis for the given $P_{\psi}/P_{\phi}$) are given by

$$I_2={1-I_1\over 1-{c-1\over c}I_1},$$ (4)

where $c=(P_{\psi}/P_{\phi}\pm 1)^2$, the plus sign being for LAMs and the minus one for SAMs. Thus,

$$I^{\rm min/max}_{\rm s}={P_{\psi}/P_{\phi}\pm 1\over P_{\psi}/P_{\phi}},$$ (5)

min/max and plus/minus being for LAMs/SAMs. Note that for LAMs $P_{\psi}/P_{\phi}>0$, whereas for SAMs $P_{\psi}/P_{\phi}>2$. The hyperbolae shown in Fig. 3 are those of c=9 (dashed line) and c=4 (dash-dot). The allowed regions are on the concave sides of the hyperbolae. The SAM region is the triangle on the left, the LAM region the strip on the right, continuing to infinity. A sphere would be at (1, 1), a flat disk at (2, 1) and (0.5, 0.5), and a thin needle at (0, 1) and $(\infty,\ \infty)$.

Some further limiting factors apply as well. One is the lightcurve amplitude, another the expected maximum elongations of the body dimensions. If, for example, one deems an elongation ratio of more than three between any two axes unlikely, $I^{\rm max}_{\rm s}\approx 5$ for LAMs, while $I^{\rm min}_{\rm s}\approx 0.55$ for SAMs (i.e., any state with $P_{\psi}/P_{\phi}<2.2$ would be declared LAM outright), and $I^{\rm max}_{\rm a}\approx 0.4$. Amplitude is not an important constraint as it is geometry-dependent, and the other constraints already bound the search region well. For example, a rough guess for a SAM near opposition would be $I^{\rm min}_{\rm a}\approx {1\over 2}{A^2-1\over A^2+1}$, where A is the relative intensity difference between maxima and minima. Since the tilt angle $\theta$ of an axisymmetric body is given by

$$\cos\theta=\left\vert {P_{\psi}\over P_{\phi}}(I_{\rm s}-1)\right\vert^{-1},$$ (6)

the amplitude of a LAM, as a basic rule, increases with $I_{\rm s}$.

It is useful first to scan the dynamical parameter space by using an ellipsoid whose dimensions are consistent with the current moments of inertia so that, the scattering law aside, eight parameters describe the whole problem. The scanning be done on a relatively loose grid, choosing a starting point from each grid cell and homing in on the local minimum with a gradient algorithm (such as Levenberg-Marquardt - see Kaasalainen et al. 2001a,b). Only a few sampling values are needed for each angle ( $\tilde\beta,\lambda,\phi_0,\psi_0$). Sparse sampling is usually sufficient for $I_{\rm s}$, $I_{\rm a}$ as well, though there may be raggedness in $\chi^2$ along the $I_{\rm s}$ - parameter axis in places. Genetic algorithms, though slower than gradient ones, usually clear the possible unsmooth bits well. The two periods behave just as the single period of a principal-axis rotator, i.e., the sampling interval for the period P should be smaller than

$$\Delta P\approx {1\over 2} {P^2\over T},$$ (7)

where T is the time span of the lightcurve set. In this way, all possible peak shifts between lightcurves will be covered. The upper and lower estimates for the trial periods are obtained from the peak widths of the power spectrum. If more than two peaks are seen, one can often infer the correct periods directly by examining the possible frequency combinations; in any case, wrongly chosen periods usually give an obviously incorrect fit. Once the ellipsoid fit has been done, one can pick the best results and, using their parameter values as initial guesses, continue by fitting a general shape to the lightcurves as in Kaasalainen et al. (2001a,b).

	Figure 4: Synthetic lightcurve of an asteroid in long-axis mode (crosses) and the corresponding curve of the obtained convex model (solid line).
Open with DEXTER

This shape can be either convex or nonconvex; the former has better stability properties, while the latter facilitates an efficient use of the inertia tensor in a regularizing function. Local minima tend to be denser in dynamical parameters for the general shape than for the consistent ellipsoid: the optimization procedure may find adjusting the shape to be locally more profitable than adjusting the dynamics, especially if the data set is small.

The main bonus of observing a precessing target is that only one observing geometry (with respect to the angular momentum vector $\vec{L}$) is already sufficient for a preliminary model if the lightcurve covers several cycles. Two geometries are needed for a unique solution. If the solar phase angle $\alpha$ at the one observing geometry is nonzero, there are always two possible mirror solutions for $\vec{L}$, with corresponding mirror-image shape solutions. One solution is the correct $\vec{L}$, while the other one is the reflection of $-\vec{L}$ with respect to the Sun-Earth-asteroid plane. If the one geometry is at opposition (or close to it), there is an infinite number of $\vec{L}$-solutions, all on a cone around the line of sight, its surface defined by the angle between the viewer and $\vec{L}$. The shape solutions and other dynamical parameters can still be correct (within a possible reflection). If at least one more geometry is available (sufficiently different from the first one, and the asteroid is not moving in the plane of the ecliptic), the symmetry breaks and there is only one solution.

In the case of short-period observations (e.g., only two lightcurves a few days apart), we are attempting a feat that would be strictly impossible for principal-axis rotators. This is why the quality of the data now plays an especially important role. As shown in Kaasalainen et al. (2001a,b), noise is not a significant factor in the case of relaxed rotators, because there always must be lightcurves from a relatively long time span and various geometries. Inversion of a few short lightcurves of a precessing object is obviously more vulnerable to noise and especially outliers. Also, the difference between models obtained with different scattering laws may now be more noticeable.

4 Examples

Inversion results of long-spanned observations of precessing asteroids are similar to those of principal-axis rotators, so I concentrate here on the case of minimal observations. A synthetic noisy lightcurve of a LAM asteroid, produced by the same object (shown in Fig. 6) and scattering law as before, and observed at the solar phase angle $\alpha=40^\circ$, is shown in Fig. 4, and the corresponding power spectrum in Fig. 5.

	Figure 5: The power spectrum of the lightcurve of Fig. 4.
Open with DEXTER

The nutation angle $\theta$ wobbles around $50^\circ$. The main peaks correspond to, from left to right, $4f_{\phi}+2f_{\psi}$, $2f_{\phi}+2f_{\psi}$, and $2f_{\phi}$, with $P_{\phi}=13.8,\ P_{\psi}=17.5$. Though the lightcurve contains only of order ten characteristic periods and the noise is about 0.04 mag, we can get a decent convex representation of the nonconvex object, as shown in Fig. 7 (Fig. 6 shown SAM-like, 7 LAM-like).

	Figure 6: The object used in generating the synthetic lightcurves.
Open with DEXTER

	Figure 7: The convex shape model obtained from the lightcurve of Fig. 4.
Open with DEXTER

The lightcurve fit is superposed in Fig. 4. A mirror shape (and angular momentum vector) would be valid as well; even a short additional lightcurve is sufficient to single out the correct solution. The differences between the correct and the obtained dynamical parameters are typically $5^\circ$-$10^\circ$ for the angles, and $3\%$-$7\%$ for the moments of inertia. A slightly incorrect scattering law does not alter the figures much because of the already considerable noise. The inversion result also remains the same, if a couple of basic periods are missed in the middle of the observation span. The general principle is to prevent aliasing effects by covering several characteristic periods densely, either in one go or by combining separate observing runs. Adding more gaps, noise or occasional outliers to the lightcurve of Fig. 4 distorts the shape until it can no longer be accepted due to the clear discrepancy between it and the obtained moments of inertia. It is possible to use the inertia tensor as a regulating factor, but generally it seems that either the data are accurate enough for a detailed analysis, or the noise quickly prevents anything but a rough solution such as a consistent ellipsoid. The dynamical parameters obtained are good even in the latter case; thus lightcurves indeed contain robust information on the object's rotational state.

An obvious example of real observations is 4179 Toutatis. The main problem with the data set of Spencer et al. (1995) is poor calibration combined with the slowness of the precession/rotation. Fitting relative brightnesses is not possible since the observing geometry changes considerably during the basic cycles. Significant new parts of rotational cycles cannot be obtained during one night, so the absolute magnitude should be calibrated very accurately. As the lightcurve fit of Hudson & Ostro (1998) (based on radar observations) shows, this is not the case, and we should expect the corresponding noise to be of order 0.1 mag. Lightcurve analysis without prior information confirms this.

Figure 8 presents two fits to the data set of Dec. 1992-Jan. 1993, obtained with the Hapke scattering model.

	Figure 8: Two model fits to the observed magnitudes of 4179 Toutatis (crosses): the solid line is produced by a consistent ellipsoid, and the dashed line by a general convex shape.
Open with DEXTER

Even though the data are lumpy and scattered, suitably filtered power spectrum analysis discloses the basic periods, as shown by Mueller et al. (2000). I simply adopted the values $P_{\psi}=7.35,\ P_{\phi}=5.41$ days as the data do not lend themselves to an accurate analysis; many nearby values would have given equally good fits. Only the observations at solar phase angles less than $90^{\circ}$ were included to minimize the uncertainties due to the scattering law. The solid line is that of a consistent ellipsoid (axial ratios $a/c=0.40,\ b/c=0.48$), and the dashed line that of a general convex shape. The key dynamical parameters of the former are $\beta=-60^\circ,\ \lambda=182^\circ,\ I_1=3.15,\ I_2=2.97$, and those of the latter $\beta=-58^\circ,\ \lambda=175^\circ,\ I_1=3.21,\ I_2=2.86$, both well consistent with the $\beta=-52^\circ,\ \lambda=180^\circ,\ I_1=3.19,\ I_2=3.01$of Hudson & Ostro (1998). The respective initial Euler angles of the two types of fit are within a few degrees from each other, except for $\psi_0$ that varies some twenty degrees as its role is not very significant, the inertia ellipsoid being almost axisymmetric.

Gaps and outliers cause the general shape result to drift far away from the correct one. This is quite evident even without the information from the radar images, since the discrepancy between the shape and the obtained moments of inertia is considerable. The lightcurve fit is much better (0.08 mag if albedo variegation allowed, 0.09 mag if not) than it ought to be (0.12 mag for the radar model). A consistent ellipsoid fits the data as well as the radar solution, and it is evident that no more details can be obtained from the lightcurve; myriads of different shape solutions fit the data within 0.12 mag, and many better than that. The dynamical parameters, however, are good, so the Toutatis data set confirms the implications of numerical simulations.

5 Conclusion

Lightcurves of precessing bodies (spanning a number of both characteristic periods) have quite distinct appearances and should usually be distinguishable - after either a preliminary or a full analysis - from other two-period lightcurves, viz. those of binary asteroids in non-synchronous rotation. The fingerprints of the rotational state are quite strong, so the dynamical parameters can usually be well inferred even from noisy data, if the observations cover at least two sufficiently different observing geometries. Some basic precession parameters can be obtained already at one geometry. If the observation span is short, the shape model is more sensitive to noise; thus, a detailed shape model requires accurate measurements and/or at least a few (long) lightcurves.

We may expect the observed tumbling bodies to be more often slow (such as Toutatis) rather than fast rotators since the lifetime of slow tumbling motion is longer. In such cases, accurate calibration of absolute brightnesses is of primary importance.

Acknowledgements

I would like to thank Stefano Mottola for valuable comments and discussions.

Appendix A: Force-free precession

Free precession has been well discussed by many authors (see, e.g., Samarasinha & A'Hearn 1991 and references therein). However, some subjects needed in the inverse problem require special attention, so I give a concise account of the dynamical problem here.

The time evolution of the Euler angles $\phi,\ \theta,\ \psi$ is easy to derive from the basic kinematic equations (see, e.g., Goldstein 1980) that can be written as

$${L_1\over I_1} \dot\phi\,\sin\theta\sin\psi+\dot\theta\cos\psi,$$ =

$${L_2\over I_2} \dot\phi\,\sin\theta\cos\psi-\dot\theta\sin\psi,$$ = (A.1)

$${L_3\over I_3} \dot\phi\,\cos\theta+\dot\psi,$$ =

where L_i are the projections of the angular momentum vector $\vec{L}$ onto the principal axes of the tumbling inertia ellipsoid; $\vec{L}$ is constant in an inertial frame in torque-free precession. The components L_i are, by the definition of the Euler angles,

$$(L_1,L_2,L_3)=(L\sin\theta\sin\psi,L\sin\theta\cos\psi,L\cos\theta).$$ (A.2)

Combining the two sets of equations, we get

 

$$\dot\phi L(I_+-I_-\cos 2\psi),$$ =

$$\dot\theta L\,I_-\sin\theta\sin 2\psi,$$ = (A.3)

$$\dot\psi \cos\theta\left({L\over I_3}-\dot\phi\right)\!,$$ =

where $I_+={1\over 2}(I_1^{-1}+I_2^{-1})$ and $I_-={1\over 2}(I_1^{-1}-I_2^{-1})$. The basic properties of the time evolution of the Euler angles are quite apparent from this set of equations: when one considers $I_-<\!<I_+$: the precession angle $\phi$ grows almost linearly, with slight undulations synchronized with the rotation angle $\psi$, while the nutation angle $\theta$ is of sinusoidal type, giving rise to fluctuations (larger than those of $\phi$) in the otherwise linear development of $\psi$. Note that for I₁,I₂<I₃ (SAM), the factor $L/I_3-\dot\phi$ in $\dot\psi$ is negative; for I₁,I₂>I₃ (LAM), it is positive. Thus the rotation and precession directions of SAMs are opposite, while LAMs both rotate and precess in the same direction. Using $\theta>\pi/2$ does not really change this, since such values of $\theta$ only correspond to interchanging the "North'' and the "South'' of the body and are thus redundant.

The semiaxes a,b,c and the corresponding principal moments of inertia of an ellipsoid are related by (setting c=1, I_c=1)

$$I_a={b^2+1\over a^2+b^2},\qquad I_b={a^2+1\over a^2+b^2},$$ (A.4)

and

$$a^2={I_b-I_a+1\over I_a+I_b-1},\qquad b^2={I_a-I_b+1\over I_a+I_b-1}\cdot$$ (A.5)

The latter equations reflect the fact following directly from the definition of the moments of inertia of a rigid body: for all real bodies, the three moments must satisfy the triangle inequalities $I_i+I_j\ge 1$, $\vert I_i-I_j\vert\le 1$.

Using the angular momentum $\vec{L}$ and the energy E,

$$E={1\over 2} \vec{L}^2 \left\lbrack\sin^2 \theta \left({\sin^... ...\psi\over I_2}\right)+{\cos^2 \theta\over I_3}\right\rbrack\!,$$ (A.6)

as integrals of motion, one can readily employ the equations of motion to express the tumbling time evolution of the inertia ellipsoid in terms of elliptic integrals and Jacobian elliptic functions (e.g., Samarasinha & A'Hearn 1991; note that their description of short-axis mode differs from the one in this paper because of the different coordinate system and period definitions). Such expressions, however, are not very practical; it is easier just to integrate Eq. (A.3) numerically in one go, tabulating the results for each desired epoch. In this way one avoids expensive function calls to each point separately. It is useful to evaluate one quantity via an elliptic integral; this is the rotation period $P_{\psi}$, given by

$$P_{\psi}=4\sqrt{I_1 I_2 I_3\over 2 E\left(I_1-I_3\right)\left(I_2-{\vec{L}^2\over 2E}\right)} \,K(k),$$ (A.7)

where K(k) is the complete elliptic integral $K(k)=\int_0^{\pi/2} (\sqrt{1-k^2\sin^2 x})^{-1}\,{\rm d}x$, and

$$k^2={\left(I_2-I_1\right)\left({\vec{L}^2\over 2E}-I_3\right)... ...\left(I_1-I_3\right)\left(I_2-{\vec{L}^2\over 2E}\right)}\cdot$$ (A.8)

The nutation period is $P_{\theta}=P_{\psi}/2$. Note that the real argument of K is actually k²; k is never needed, but I use the notation K(k) to conform with the historical custom. Also, k² can by all means be negative: the allowed values are $k^2\le 1$(the sign of k² depends on the chosen ordering of the moments of inertia labelled $I_1,\ I_2$ - the resulting $P_{\psi}$ is naturally independent of this ordering). If k²>1 (or the square root term of Eq. (A.7) is imaginary), one has chosen the wrong extremal axis for rotation and should switch the z-axis of the coordinate system to the other extremal axis. (This never happens in iteration, since the formulae given here enable one to stay precisely in the allowed region for SAMs or LAMs.) If needed, the new Euler angles of the correct coordinate system are

$$\phi' \varphi_z+{\pi\over 2},$$ =

$$\theta' \vartheta_z,$$ = (A.9)

$$\psi' {\rm sign}(\cos\vartheta_x){\rm acos} (\sin\vartheta_x\cos(\varphi_x-\phi')),$$ =

where $(\varphi_z,\vartheta_z)$ are the longitude and polar angle of $\hat{\vec{z}}=M^{-1}(0,0,1)^T$, $M=R_x(\pi/2)R_z(\psi)R_x(\theta)R_z(\phi)$. Similarly, $(\varphi_x,\vartheta_x)$ apply to $\hat{\vec{x}}=M^{-1}(1,0,0)^T$. The above formulae for $P_{\psi}$ are correct for all configurations, regardless of the ordering of I_i, as long as I₃ is either the smallest or the largest one of the three (and is the correct extremal axis of rotation).

Using the constant

$$\Delta\phi=\int_0^{P_{\psi}} \dot\phi \,{\rm d}t =2\,\int_0^{P_{\psi}/2} \dot\phi \,{\rm d}t$$ (A.10)

(note that $\Delta\phi$ is independent of L, whose occurrences in $P_{\psi}$ and $\dot\phi$ cancel each other out), we can define

$$P_{\phi}= 2\pi {P_{\psi}\over \Delta\phi}\cdot$$ (A.11)

Thus, for example, the initial nutation angle $\theta_0$corresponding to an observed ratio R of $P_{\psi}/P_{\phi}$ and a given $\psi_0$ is found by solving numerically the equation

$$\int_0^{P_{\psi}(\theta_0)} \dot\phi(\theta_0) \,{\rm d}t=2\pi R$$ (A.12)

for $\theta_0$.

$P_{\psi}\rightarrow\infty$ when either k²=1, i.e., ${L^2\over 2E}=I_1$, or ${L^2\over 2E}=I_2$. In terms of the nutation angle $\theta$, this happens when

$$\cos^2\theta_{\rm max}= {\rm max}\left\lbrack{I_i^{-1}-I_0\over I_3^{-1}-I_0}\right\rbrack,$$ (A.13)

where $I_0=\sin^2\psi/I_1+\cos^2\psi/I_2$ and i=1,2. On the other hand, $\Delta\phi$ (or $P_{\psi}/P_{\phi}$) attains its minimum when $\theta=0$, i.e., in principal-axis rotation. In this case,

$$\Delta\phi_{\rm min}=\left\lbrack\sqrt{I_1I_2\over (I_1-I_3)(I_2-I_3)}\,\pm 1 \right\rbrack2\pi,$$ (A.14)

where the plus sign is for $I_1,\ I_2<I_3$ and the minus sign for $I_1,\ I_2>I_3$. Note that the square root is $\ge$1 for real bodies, so $P_{\psi}/P_{\phi}>0$ for $I_1,\ I_2>I_3$ and $P_{\psi}/P_{\phi}>2$for $I_1,\ I_2<I_3$. From the inequality (A.14) we directly obtain the limiting hyperbolae in $I_1,\ I_2$-plane discussed in Sect. 3.

Knowing $P_{\psi}$ enables one to jump easily over the gaps between lightcurves, thus keeping the integration time short (for maximal saving of integration time, all epochs can be projected within one time interval of $P_{\psi}$). Thus

 

$$\phi(t) \phi(t_0)+n\Delta\phi,$$ =

$$\theta(t) \theta(t_0),$$ = (A.15)

$$\psi(t) \psi(t_0)\pm n 2\pi,$$ =

when t is some suitable epoch at $t=t_0+n P_{\psi}$ with n an integer. When using gradient-based optimization algorithms, one needs the derivatives of the Euler angles with respect to the dynamical parameters; these are obtained by integrating the corresponding derivatives of $\dot\phi,\ \dot\theta,\ \dot\psi$ simultaneously with Eq. (A.3). The derivatives contain a component growing linearly with time, so using Eq. (A.15) to derive the jumps for the derivatives also serves to keep the integration result stable and accurate over long periods of time.

References

 
• Goldstein, H. 1980, Classical Mechanics (Second Edition) (Addison-Wesley, Reading, Massachusetts) In the text  
• Harris, A. W. 1994, Icarus, 107, 209 In the text NASA ADS  
• Hudson, R. S., & Ostro, S. J. 1998, Icarus, 135, 451 In the text NASA ADS  
• Kaasalainen, M., & Torppa, J. 2001a, Icarus, in press In the text  
• Kaasalainen, M., Torppa, J., & Muinonen, K. 2001b, Icarus, in press  
• Kryszczynska, A., Kwiatkowski, T., Breiter, S., & Michalowski, T. 1999, A&A, 345, 643 In the text NASA ADS  
• Mottola, S., & Lahulla, F. 2000, Icarus, 146, 556 In the text NASA ADS  
• Mueller, B. E. A., Belton, M. J. S., & Samarasinha, N. H. 2000, BAAS, 32, 1003 In the text  
• Pravec, P., & Hahn, G. 1997, Icarus, 127, 431 In the text NASA ADS  
• Samarasinha, N. H., & A'Hearn, M. F. 1991, Icarus, 93, 194 NASA ADS  
• Spencer, J. R., et al. 1995, Icarus, 117, 71 In the text NASA ADS