Complex variables approach to the short-axis-mode rotation of a rigid body

The rotation of a triaxial rigid body in the absence of external torques is known to be integrable [1, 2]. In particular, the canonical transformation to Andoyer variables [3] reduces the free rigid body rotation to an integrable, one degree of freedom Hamiltonian, which immediately shows the preservation of the total angular momentum and allows for the representation of the possible solutions by contour plots of the reduced Hamiltonian [4]. However, because the solution to the torque-free motion depends on elliptical integrals and elliptic functions, the closed form solution is replaced in practical applications by corresponding expansions in trigonometric functions truncated to some order —as it is customary for integrable problems in which the solution depends on special functions and may further require series inversion [5, 6]. In particular, useful expansions by Kinoshita apply for the case in which the triaxiality of the rigid body is small [7] and for the case in which the rigid body rotation is close to the axis of either maximum or minimum momentum of inertia [8, 9] —the order of these expansions was later extended by other authors [10]. Alternatively, the expansions in trigonometric functions of the free rigid body solution can be directly constructed using perturbation theory [11, 12], an approach that systematizes the computation of higher orders of the expansions and eases the construction of perturbation solutions in the presence of external torques [13], as opossed to the typical first order approach [14, 15, 16].

Perturbation approaches to the torque-free motion of a rigid body start from the decomposition of the free rigid body Hamiltonian in Andoyer variables into a “main problem” and a perturbation term, the solution of the main problem being achievable in elementary functions, contrary to the special functions required in the solution of the full problem. When the triaxiality is small, the selection of either the axisymmetric case or the spherical rotor as the main problem results in a completely reduced zeroth order Hamiltonian that only depends on momenta of the canonical set of Andoyer variables, a fact that simplifies application of the perturbation method because, then, Andoyer variables result to be the action-angle variables of the main problem [11]. Practical application of the perturbative scheme based on the uniaxial model to the case of a non-rigid Earth (rigid mantle and liquid core) have been made in [17].

On the contrary, while the perturbative character of short-axis-mode (SAM) rotation can be shown directly in Andoyer variables, the consequent complete reduction of the main problem by solving the Hamilton-Jacobi equation is required in order to set up an efficient perturbative integration scheme [12]. Indeed, the main problem of SAM rotation is now a simplification of the free rigid body Hamiltonian, contrary to a reduction [18, 19]. That is, it involves the same Andoyer variables as the free rigid body problem, which include coordinates as well as their conjugate momenta, and, in consequence, Andoyer variables are not action-angle variables of the zeroth order Hamiltonian.

As far as the complete reduction of a Hamiltonian is unique [20], the canonical transformation used in finding the completely reduced Hamiltonian is irrelevant to some extent, and the standard approach is to find the action-angle variables of the integrable problem by solving the classical Hamilton-Jacobi equation or variations of it [21]. However, action-angle variables commonly involve singularities for particular configurations in which some angles may be not defined, and become ill-defined for values close to those causing the singularities. The classical instance comes from the orbital motion of a particle, where the action-angle reduction of the Kepler problem by solving the Hamilton-Jacobi equation is achieved in Delaunay variables, in which the argument of the periapsis is not defined for circular orbits, and the argument of the node is not defined in the case of equatorial orbits. Nevertheless, these singularities are of “virtual” nature [22] and can be avoided using non-singular variables. Poincaré canonical coordinates and conjugate momenta avoid the singularities in the orbital motion [23], yet non-canonical sets of variables are commonly preferred in different applications [24, 25, 26].

In the case of rotational motion, Andoyer canonical variables, which are the action angle variables of the uniaxial rotation, show singularities in different configurations. Alternative variables [27, 28] remove the singularities of Andoyer variables but at the cost of creating singularities at different locations [29].

Then, in spite of the standard approach used in [12] of finding the action-angle variables that completely reduce the main problem of SAM rotation is perfectly correct, it happens that the variables used share the same deficiencies of Andoyer variables. In particular, they may become ill-defined when the axis of instantaneous rotation evolves very close to the axis of maximum (or minimum) inertia. Therefore, in practice one must to take some care when applying the solution to this kind of motion, which is precisely the case in which the separation of the free rigid body rotation into the main problem of SAM rotation and a perturbation applies.

On the other hand, it will be shown in Section 2.3 that the perturbative arrangement of the free rigid body Hamiltonian in the case of SAM rotation is immediately disclosed when using non-singular variables of Poincaré type, cf. [30]. In these variables, the free rigid body Hamiltonian takes the form of the Hamiltonian of the simple harmonic oscillator disturbed by additional quartic polynomial terms. Hence, since it is well known that the use of complex variables renders very efficient the construction of higher order analytical solutions to perturbed harmonic oscillators (see, for instance, [31, 32, 33]), an additional transformation to complex variables is carried out in Section 3 that converts the integration of the free rigid body Hamiltonian by perturbations into a simple exercise of polynomial algebra.

The polynomial structure of the perturbation is not preserved, in general, when perturbation torques are taken into account. However, it is shown in Section 4 that, when the rotation is close enough to the axis of maximum inertia, the gravity-gradient perturbation can be easily tackled within the same perturbative scheme.

Finally, it is worth mentioning that how much fruitful expanding the special functions in which the solution of the torque-free motion inherently depends into trigonometric ones can be, there are specific cases in which approximate analytical solutions to perturbed attitude motion can be achieved directly in closed form [34, 35]. These kinds of solutions widen applicability by avoiding the constraint to particular physical characteristics or dynamical configuration, and have been recently proposed as an appealing alternative for practical application to actual problems of astrodynamics [36].

The free rigid body Hamiltonian is plainly stated in terms of Andoyer [3] canonical variables in the form, cf. Eq. (4) of [37],

H0=1Asin2⁡ν+1Bcos2⁡ν12M2−N2+12CN2,$$\begin{array}{} \displaystyle \mathscr{H}_0=\left(\frac{1}{A}\sin^2\nu+\frac{1}{B}\cos^2\nu\right)\frac{1}{2}\left(M^2-N^2\right)+\frac{1}{2C}N^2, \end{array} $$(1)

where the constants A ≤ B ≤ C define the principal moments of inertia, and the variables (λ,μ,ν,Λ,M,N) define the node of the invariable plane on the inertial x, y plane, the node of the equatorial plane of the rigid body over the invariable plane, the component of the rotation of the body around its axis of maximum inertia, the projection of the angular momentum vector along the inertial z axis, the total angular momentum, and the projection of the angular momentum vector along the body axis of maximum inertia, respectively. The fact that the variables λ and μ are cyclic in Hamiltonian (1) immediately shows that the total angular momentum M as well as its projection along the z axis of the inertial frame Λ are integrals. Besides, because Λ is also ignorable in Eq. (1), the node of the invariable plane remains fixed in the inertial x-y plane. Therefore, the free rigid body Hamiltonian in Andoyer variables is only of one degree of freedom in the canonical pair (ν,N).

It belongs to Sadov [38] the merit of finding the transformation to action-angle variables that completely reduces Hamiltonian (1). His approach resorts to a complete solution of the Hamilton-Jacobi equation and involves the elliptic integrals of the first and third kinds, both of them in their incomplete and complete versions. Because the transformation from Andoyer to action-angle variables cannot be obtained explicitly, the completely reduced Hamiltonian must remain also as an implicit function of the action variables. A modern, abridged re-derivation of Sadov’s solution, which does not need to deal explicitly with the time when solving the Hamilton-Jacobi equation, can be found in [35]. This solution is summarized in the Appendix for completeness.

Working independently of Sadov, Kinoshita [7] succeeded also in finding the transformation of the free rigid body Hamiltonian to action-angle variables, yet he expressed his solution in terms of the Heuman’s Lambda function rather than the elliptic integral of the third kind. However, it must be reminded that the characteristic parameter in which the elliptic integral of the third kind depends upon, is defined as a function of the principal moments of inertia of the free rigid body, and is always negative. Because of that, the validity of the elliptic integral of the third kind that appears in the complete reduction of the free rigid body Hamiltonian is constrained to one of the “circular” cases [39], and hence it is customarily evaluated in terms of the Heuman’s Lambda function. Therefore, both, Sadov’s and Kinoshita’s solutions of the torque-free motion must be considered equivalent.

To make the use of the free rigid body solution practical in the study of the rotation under external torques, the closed form solution is customarily expanded in terms of trigonometric functions. Both authors, Sadov and Kinoshita, provided the necessary expansions in powers of the Jacobi’s nome, which is a function of the elliptic modulus that improves convergence over the direct use of the elliptic modulus in the expansions. Alternatively, the expansion of the closed form solution can be avoided in the particular cases in which some small quantity is identified in the original Hamiltonian. Then, the free rigid body Hamiltonian can be rearranged as a perturbation Hamiltonian and the solution of the integrable problem is directly attacked by perturbation methods [6].

Regrettably, because the modulus of the elliptic functions that solve the torque free motion depends on the energy of each particular solution, which is in fact the Hamiltonian, neither the Jacobi’s nome nor the elliptic modulus are useful parameters to be identified in the original Hamiltonian in order to apply a perturbation approach. Still, as shown in the Appendix, it happens that this elliptic modulus can be decomposed into a product that splits the physical and dynamical characteristics of the motion, and, luckily, each kind of feature is easily recognized in Hamiltonian (1) after simple rearrangement.

For perturbed harmonic oscillators, it is well known that the use of complex variables makes the whole perturbation approach very efficient, cf. [31, 32, 33]. Therefore, the perturbation Hamiltonian given by Eqs. (15) and (19) is next reformulated in complex variables.

3.2

Free rigid body solution by perturbations

The perturbation procedure starts writing the free rigid body Hamiltonian in the form of the Taylor series expansion, viz.

H0=∑n≥0εnn!Hn,0(u,U),$$\begin{array}{} \displaystyle \mathscr{H}_0=\sum_{n\ge0}\frac{\varepsilon^n}{n!}H_{n,0}(u,U), \end{array} $$

where H_0,0 ≡ ℳ as given in Eq. (25), H_1,0 ≡ 𝒫 as given by Eq. (26), H_n,0 = 0 for n ≤ 2, and ε is a formal small parameter. Straightforward computations based on the previously discussed properties of the Lie operator, yield the normalized Hamiltonian, in new, prime variables

K=κ2CV21κ3/2+∑n≥0(−α)nκn/2pniu′U′Vn+1,$$\begin{array}{} \displaystyle \mathscr{K}=\frac{\sqrt{\kappa}}{2C}V^2\left[\frac{1}{\kappa^{3/2}}+ \sum_{n\ge0}\frac{(-\alpha)^n}{\kappa^{n/2}}p_n\left(\frac{iu'U'}{V}\right)^{n+1}\right], \end{array} $$(30)

where p_n are polynomials in the triaxiality coefficient β. The first few triaxiality polynomials are p₀ = 2, and

p1=1,p2=β2p3=58β2,p4=332β23β2+8,p5=732β25β2+4,p6=1128β245β4+354β2+128,p7=91024β2265β4+650β2+128,p8=58192β2953β6+14888β4+17120β2+2048,p9=118192β24075β6+20212β4+13104β2+1024,$$\begin{array}{} \displaystyle p_1 = 1, \\\displaystyle p_2 = \beta^2 \\\displaystyle p_3 = \frac{5}{8}\beta^2, \\\displaystyle p_4 = \frac{3}{32}\beta^2\left(3\beta^2+8\right), \\\displaystyle p_5 = \frac{7}{32}\beta^2\left(5\beta^2+4\right), \\\displaystyle p_6 = \frac{1}{128}\beta^2\left(45\beta^4+354\beta^2+128\right), \\\displaystyle p_7 = \frac{9}{1024}\beta^2\left(265\beta^4+650\beta^2+128\right), \\\displaystyle p_8 = \frac{5}{8192}\beta^2\left(953\beta^6+14888\beta^4+17120\beta^2+2048\right), \\\displaystyle p_9 = \frac{11}{8192}\beta^2\left(4075\beta^6+20212\beta^4+13104\beta^2+1024\right), \end{array} $$

which, as expected, are the same as those in Table 2 of [12] after adjusting subindices and scaling by β².

The transformation from prime to original variables

u=u′+∑n≥1εnn!u0,n(u′,U′),U=U′+∑n≥1εnn!U0,n(u′,U′),$$\begin{array}{} \displaystyle u=u'+\sum_{n\ge1}\frac{\varepsilon^n}{n!}u_{0,n}(u',U'), \qquad U=U'+\sum_{n\ge1}\frac{\varepsilon^n}{n!}U_{0,n}(u',U'), \end{array} $$

is obtained by successive evaluations of Deprit’s triangle

ξn,q=ξn+1,q−1+∑0≤m≤nnm{ξn−m,q−1;Sm+1},$$\begin{array}{} \displaystyle \xi_{n,q}=\xi_{n+1,q-1}+ \sum_{0\le{m}\le{n}}{n\choose{m}}\,\{\xi_{n-m,q-1};S_{m+1}\}, \end{array} $$(31)

using the generating function 𝒮 = ∑_n≥0(εⁿ/n!) S_n+1 where terms

Sn=αnκn/2Vnβu2+U2sn,$$\begin{array}{} \displaystyle S_n=\frac{\alpha^n}{\kappa^{n/2}V^n}\beta\left(u^2+U^2\right)s_n, \end{array} $$(32)

are computed using Eq. (29). The first few s_n are

s1=18uU,s2=−<!-- - -->i4u2U2,s3=i6424iβ<!-- ? -->2+2uU+5β<!-- ? -->(U2−<!-- - -->u2)u2U2,s4=3642i57β<!-- ? -->2+32uU+β<!-- ? -->9β<!-- ? -->2+20(U2−<!-- - -->u2)u3U3,s5=1642343β<!-- ? -->4+2024β<!-- ? -->2+480u2U2−<!-- - -->11β<!-- ? -->2β<!-- ? -->2+2u4+U4−<!-- - -->6iβ<!-- ? -->147β<!-- ? -->2+100uUU2−<!-- - -->u2u3U3,s6=1564i−<!-- - -->2909β<!-- ? -->4+1588β<!-- ? -->2+192u2U2+β<!-- ? -->263β<!-- ? -->2+44u4+U4+2iβ<!-- ? -->91β<!-- ? -->4+653β<!-- ? -->2+200uUU2−<!-- - -->u2u4U4,s7=45i2048β<!-- ? -->124067β<!-- ? -->4+272282β<!-- ? -->2+44800u2U2U2−<!-- - -->u2−<!-- - -->β<!-- ? -->3283β<!-- ? -->2+186U6−<!-- - -->u6−<!-- - -->32iβ<!-- ? -->2116β<!-- ? -->4+917β<!-- ? -->2+308uUu4+U4+64i720β<!-- ? -->6+8345β<!-- ? -->4+6496β<!-- ? -->2+448u3U3u4U4,s8=10520483β<!-- ? -->51619β<!-- ? -->6+725628β<!-- ? -->4+730320β<!-- ? -->2+71680u2U2U2−<!-- - -->u2−<!-- - -->3β<!-- ? -->3539β<!-- ? -->4+4396β<!-- ? -->2+1488U6−<!-- - -->u6−<!-- - -->4iβ<!-- ? -->242087β<!-- ? -->4+104446β<!-- ? -->2+19712uUu4+U4+4i391113β<!-- ? -->6+1385374β<!-- ? -->4+578816β<!-- ? -->2+24576u3U3u5U5.$$\begin{array}{} \displaystyle s_1 = \frac{1}{8}uU, \\\displaystyle s_2 =-\frac{i}{4}u^2U^2, \\\displaystyle s_3 = \frac{i}{64}\left[{24i\left(\beta ^2+2\right)uU+5\beta(U^2-u^2)}\right]u^2U^2, \\\displaystyle s_4 = \frac{3}{64}\left[{2i\left(57\beta^2+32\right) u U+\beta\left(9\beta^2+20\right)(U^2-u^2)}\right]u^3U^3, \\\displaystyle s_5 = \frac{1}{64} \left[ 2 \left(343\beta^4+2024\beta^2+480\right) u^2U^2 -11 \beta^2\left(\beta^2+2\right)\left(u^4+U^4\right) \right. \\\displaystyle \qquad\, \left.-6i\beta\left(147\beta^2+100\right)uU\left(U^2-u^2\right) \right]u^3U^3, \\\displaystyle s_6 = \frac{15}{64}i \left[ -2\left(909\beta^4+1588\beta^2+192\right)u^2U^2 +\beta^2\left(63\beta^2+44\right)\left(u^4+U^4\right) \right. \\\displaystyle \qquad\, \left. +2i\beta\left(91\beta^4+653\beta^2+200\right)uU\left(U^2-u^2\right) \right]u^4U^4, \\\displaystyle s_7 = \frac{45i}{2048} \left[ \beta\left(124067 \beta ^4+272282 \beta ^2+44800\right) u^2 U^2\left(U^2-u^2\right)-\beta ^3 \left(283 \beta ^2+186\right) \left(U^6-u^6\right)\right. \\\displaystyle \qquad\, \left. -32 i \beta ^2\left(116 \beta ^4+917 \beta ^2+308\right) u U \left(u^4+U^4\right) +64 i \left(720 \beta ^6+8345 \beta ^4+6496 \beta ^2+448\right) u^3 U^3 \right]u^4U^4, \\\displaystyle s_8 = \frac{105}{2048} \left[ 3\beta\left(51619\beta^6+725628\beta^4+730320\beta^2+71680\right)u^2U^2\left(U^2-u^2\right)\displaystyle \right. \\\displaystyle \qquad\, -3\beta^3\left(539\beta^4+4396\beta^2+1488\right)\left(U^6-u^6\right)-4i\beta^2\left(42087\beta^4+104446\beta^2+19712\right)uU\left(u^4+U^4\right) \\\displaystyle \qquad\, \left. +4i\left(391113\beta^6+1385374\beta^4+578816\beta^2+24576\right)u^3U^3 \right]u^5U^5. \end{array} $$

The construction of the perturbation solution by Lie transforms results extremely efficient in this way. As a demonstration of the performance, it has been extended to the order 35 in the small parameter ϵ. Thanks to the use of Eqs. (30) and (32), the most computationally expensive part of the procedure at each perturbation order is the evaluation of the Poisson brackets, which are needed for the computation of the known terms of the homological equation and also for filling Deprit’s triangle.

This fact is illustrated in Fig. 1 where the label H_n means the time spent in evaluating each Hamiltonian term n in Eq. (30), S_n the time spent in evaluating each term of the generating function using Eq. (29), “known” means the time spent in computation of the known terms of the homological equation at each step of the perturbation approach, “triangle” means the time spent in filling the Deprit’s triangle with the new terms computable at each step, and “total” is the total time spent at each step, which is roughly the sum of the previous times plus some minor operations required by the algorithm. The computations have been carried out with Wolfram Mathematica 9 in a MacBook Pro running under macOS Mojave, version 10.14, with a 2.8 GHz Intel Core i7 processor and 16 GB of RAM.

Fig. 1

To further emphasize the efficiency of the selection procedure used with the complex variables approach, as opposite to the usual, time consuming, averaging procedures, these times are shown in Fig. 2 relative to the time required in choosing the different terms of the transformed Hamiltonian in prime variables in Eq. (30). In this figure it is noted that the time spent in the computation of the generating function remains between 5 and 10 times higher than the time spent in the computation of the new Hamiltonian term except for the lower orders, where the ratio can be lower. On the other hand, the time spent in the computation of the known terms of the homological equation is similar to the time spent in filling Deprit’s triangle at each step, and it grows high with the order of the perturbation theory. In particular, we checked that the computation of H_n grows roughly at a quadratic rate with n, while the computation of S_n grows approximately at a cubic rate with n, and the other times grow close to a quintic rate with the order of the perturbation theory.

Figure 2

The perturbation approach based on the main problem of SAM rotation is also feasible for motion under external torques. In the particular case of gravity-gradient perturbations due to a distant body, one can make the simplifying assumption that the disturbing body moves with Keplerian motion. If, besides, the orbital and inertial planes are assumed to match, the gravity-gradient effect is derived from the potential

D=−n2a32r3(C−B)(1−3γ32)−(B−A)(1−3γ12),$$\begin{array}{} \displaystyle \mathscr{D}=-\frac{n^2a^3}{2r^3}\left[(C-B)(1-3\gamma_3^2)-(B-A)(1-3\gamma_1^2)\right], \end{array} $$(33)

where r, a, and n, are the radius, semi-major axis, and mean motion, respectively, of the Keplerian orbit in which the distant body evolves, and γ_i, i = 1,2, 3, are the direction cosines of the line joining the centers of mass of the rotating body and the distant body, which are expressed in Andoyer variables by simple rotations. Namely,

γ1γ2γ3=R3(ν)∘R1(J)∘R3(μ)∘R1(I)∘R3(ϑ)∘100,$$\begin{array}{} \displaystyle \left(\begin{array}{c} \gamma_1 \\\gamma_2 \\\gamma_3 \end{array}\right) =R_3(\nu)\circ{R}_1(J)\circ{R}_3(\mu)\circ{R}_1(I)\circ{R}_3(\vartheta)\circ \left(\begin{array}{c} 1 \\0 \\0 \end{array}\right), \end{array} $$(34)

where I = arccos(Λ/M), J = arccos(N/M), and ϑ = λ −ϕ, ϕ is the true anomaly of the distant body, and the simplifying assumption that the inertial plane is the same as the orbital plane has been made.

After carrying out the rotations involved in Eq. (34), it is found that the majority of perturbation terms that comprise Eq. (33) are factored by sin J or sin²J (see [13], for instance). Besides, typical values of the gravity-gradient perturbation for solar system bodies are small when compared with the torque-free rotation, say below 10⁻⁶. Therefore, in those cases in which the inclination angle J is small, terms O(sin⁡12J)$\begin{array}{} \mathscr{O}(\sin\frac{1}{2}J) \end{array} $ can be neglected. In this approximation, the only relevant terms of the gravity-gradient perturbation, in the new variables, are simply

D=−n24a3r3C−A+B22−3sin2⁡I+3sin2⁡Icos⁡2ϑ+34(B−A)(1−cosI)2cos⁡(2ϑ−2g)+2sin2⁡Icos⁡2g+(1+cosI)2cos⁡(2ϑ+2g).$$\begin{array}{} \displaystyle \mathscr{D}=-\frac{n^2}{4}\,\frac{a^3}{r^3}\left\{ \left(C-\frac{A+B}{2}\right)\left(2-3\sin^2I+3\sin^2I\cos2\vartheta\right) \right. \\\displaystyle \qquad\,\left. +\frac{3}{4}(B-A)\left[(1-\cos{I})^2\cos(2\vartheta-2g)+2\sin^2I\cos2g+(1+\cos{I})^2\cos(2\vartheta+2g)\right]\right\}. \end{array} $$(35)

Notably, to this order of approximation, Eq. (35) does not depend either on ν or J. Hence, it is free also from x and X, as follows from Eqs. (12)–(13), and, in consequence, it is free from the complex variables u and U, as follows from Eqs. (20) and (21).

Now, the full Lie derivative in Eq. (27) is involved in the solution of the homological equation because g, and, therefore v is present in the perturbation. Still, the first summand in the square brackets of Eq. (27) vanishes for terms of the form F(v,V,uU). Hence, because Eq. (35) is made precisely of terms of this type, and in view of the form of the second summand in the square brackets of Eq. (27), which only includes a factor uU, the homological equation is easily solved. Indeed, if we choose the new Hamiltonian term

〈D〉=12π∫02πDdg=−n24a3r3C−A+B22−3sin2⁡I+3sin2⁡Icos⁡2ϑ,$$\begin{array}{} \displaystyle \langle\mathscr{D}\rangle=\frac{1}{2\pi}\int_0^{2\pi}\mathscr{D}\,\mathrm{d}g= -\frac{n^2}{4}\,\frac{a^3}{r^3}\left(C-\frac{A+B}{2}\right)\left(2-3\sin^2I+3\sin^2I\cos2\vartheta\right), \end{array} $$

and assume that there is not coupling with the previous terms of the perturbation theory, a particular solution of the homological equation for the order n corresponding to this term, is

Sn=−9n2256a3r3(B−A)Cκiκ3/2uU−V(1+c)2sin⁡(2ϑ+2κv)+2s2sin⁡2κv−(1−c)2sin⁡(2ϑ−2κv).$$\begin{array}{} \displaystyle S_n=-\frac{9n^2}{256}\frac{a^3}{r^3}\frac{(B-A)C\sqrt{\kappa}}{i\kappa^{3/2}uU-V}\left[(1+c)^2\sin(2\vartheta+2\sqrt{\kappa}v)+2s^2\sin2\sqrt{\kappa}v-(1-c)^2\sin(2\vartheta-2\sqrt{\kappa}v)\right]. \end{array} $$(36)

Short-axis mode rotation of a free rigid body is naturally decomposed into a main problem and a perturbation, a fact that leads to the straightforward integration of the rotation by perturbation series. When using non-singular variables of the Poincaré type, the main problem has the form of a harmonic oscillator, whose frequency is related to the triaxiality of the rigid body, whereas the perturbation is a quartic polynomial. Then, the use of complex variables makes the construction of the perturbation solution trivial. The polynomial character of the perturbation does not persist, in general, when the motion is affected by external torques. However, when the rotation is close to the axis of maximum inertia, the gravity-gradient perturbation can also be approached in complex variables.