Analytical solution of the Colombo top problem

The Colombo top is a basic model in the rotation dynamics of a celestial body moving on a precessing orbit and perturbed by a gravitational torque. The paper presents a detailed study of analytical solution to this problem. By solving algebraic equations of degree 4, we provide the expressions for the extreme points of trajectories as functions of their energy. The location of stationary points (known as the Cassini states) is found as the function of the two parameters of the problem. Analytical solution in terms of the Weierstrass and the Jacobi elliptic functions is given for regular trajectories. Some trajectories are expressible through elementary functions: not only the homoclinic orbits, as expected, but also a special periodic solution whose energy is equal to that of the first Cassini state (unnoticed in previous studies).

precession cycle. Two centuries later, Tisserand (1891) rephrased these postulates as laws. The three Cassini laws of Tisserand state that: (1) rotation and orbital periods are equal, (2) the rotation axis has a constant inclination to the ecliptic, and (3) the three axes are coplanar. The second law differs from the original Cassini's statement, but in view of the third law, the difference is unimportant.
With the advent of the Newtonian dynamics, the question arose if the Cassini's model is consistent with equations of motion. This was a part of the problem issued by the French Académie Royale des Sciences for the Prize of 1764. In his prize dissertation and in a later work, Lagrange (1764Lagrange ( , 1780) demonstrated that the state described by Cassini is an equilibrium of the associated differential system and studied small librations in its vicinity.
The work of Colombo (1966) brought a new understanding of the Cassini laws and motion near the stationary configuration which they describe in the specific case of the Moon. In particular, Colombo demonstrated that the second and the third laws are conceptually independent from the first law and themselves serve as a basis of an interesting dynamical problem which describes the long-term evolution of the spin axis of an arbitrary rotating body. He thus dropped from his analysis the assumption of the direct spin-orbit resonance, but kept the assumption of general precession of the orbital plane due to perturbations (either caused by the oblate central body, or by other masses in the system). He showed that the long-term dynamics of the spin axis can be described by a simple, one-dimensional problem assuming the orbital node performs a uniform precession and the inclination remains constant. Its stationary points represent generalizations of the Cassini second and third laws.
It was soon understood that the Colombo problem is a very suitable starting point for analysis of the obliquity evolution of terrestrial planets, even when the orbital node and inclination undergo more complex evolution. A fascinating example are studies of Mars obliquity variations in relation to this planet's past paleoclimate, starting with Ward (1973Ward ( , 1974. Tides or internal process may additionally change some of the system's parameters, a situation relevant to all terrestrial planet studies, including the Moon-see Peale (1974), Ward (1975Ward ( , 1982, or Ward and de Campli (1979), to mention just few examples of early works. Later studies of Laskar and colleagues made a masterful use of detailed knowledge of planetary long-term dynamics and its implications on secular evolution of their spin axes (e.g., Correia and Laskar 2001, and many other with more technical details). Following earlier hints, mentioned already in Harris and Ward (1982), applications to giant planets were also developed in the past two decades (e.g., Ward and Hamilton 2004;Hamilton and Ward 2004;Ward and Canup 2006;Boué et al. 2009;Vokrouhlický and Nesvorný 2015;Brasser and Lee 2015;Rogoszinski and Hamilton 2020).
Beyond planets and satellites, studies of secular spin evolution of asteroids flourished recently, especially after Vokrouhlický et al. (2003) applied it to explain space parallelism of spin axes of large members in the Koronis family (see also Vokrouhlický et al. (2006)). Further applications include spin states of exoplanets (e.g., Atobe et al. 2004;Atobe and Ida 2007;Saillenfest et al. 2019), or artificial satellites and space debris (e.g., Efimov et al. 2018).
Taken altogether, we note that the backbone of all these studies is the basic Colombo model. Interestingly, a systematic mathematical treatment of this elegant Hamiltonian problem has not been significantly advanced beyond the state of art dating back to Henrard and Murigande (1987), and Henrard (1993). These classical works focused on the aspects directly related with the probability of capture into different phase space zones when a slow parameter evolution drives the system across the homoclinic orbits; hence, they paid no attention to regular orbits. One can also recognize the Colombo top in an anonymous quadratic Hamiltonian treated by Lanchares and Elipe (1995), who focused on the qualitative study of its parametric bifurcations.
The main reason to seek for the complete analytical solution of the Colombo top is its significance for further studies of more realistic, perturbed problems. Be it analytical perturbation techniques, or numerical splitting methods, the knowledge of explicit time dependence of the Colombo top motion is crucial. The present work is divided in two principal parts: Sect. 2 explores the problem using purely geometric and algebraic tools, whereas the integration of equations of motion is considered in Sect. 3. In other words, Sect. 2 concerns integral curves, that become time-dependent trajectories in Sect. 3.
The formulation the problem is given in Sect. 2.1, where we introduce two sets of variables: traditional x, y, z and shifted X , Y , Z . Throughout the text, we switch between the two sets, depending on convenience. In Sect. 2.2, the geometric construction of the integral curves is shown; intersections of the curves with the x = X = 0 meridians (their extremities in z) are found in Sect. 2.3 , expressed in terms of the energy constant and of the two parameters a, b. Some of these can be critical points (the Cassini states), which allows the expression in terms of the parameters only-given in Sect. 2.4. The information gained allows to partition the phase space and distinguish three types of the Colombo top problem.
In Sect. 3, we first provide a universal solution in terms of the Weierstrass elliptic function ℘ (Sect. 3.1), comparing various formulations of the same result. But since the Weierstrass function behaves differently in various domains of the phase space, we reformulate the solution in terms of the more common Jacobi elliptic functions (Sect. 3.2). Finally, the specific trajectories that admit solutions in terms of elementary functions are presented in Sect. 3.3. The closing Sect. 4 summarizes the results and their implications.

Equations of motion
The basic assumptions leading to the Colombo top problem involve a rigid body on an orbit around some distant primary. The orbital motion might be called Keplerian, but with one notable addition: the orbital plane rotates uniformly around some fixed axis in the inertial space with the angular rate μ. The body is assumed to rotate in the lowest energy state, namely about the shortest principal axis of its inertia tensor. Following the work of Colombo (1966), we assume the rotation period is not in resonance with the revolution period about the primary (see, e.g., Peale 1969 for generalizations to this situation). Considering the quadrupole torque due to the gravitational field of the primary, and averaging over both orbital revolution about the center and rotation cycle about the spin axis, one easily realizes that the total angular momentum of rotation (or, equivalently, the angular velocity ω of rotation) is conserved. The whole dynamical problem then reduces to the analysis of the motion of the unit vector r of rotation pole in space.
In order to introduce fundamental astronomical parameters of relevance, let us for a moment assume μ = 0, i.e., the orbital plane is fixed in the inertial space. In this simple case, r performs regular precession about the fixed direction K of the orbital angular momentum with a frequency μ r = −α cos ε, where α is the precession constant and ε is the obliquity, namely the angle between r and K (cos ε = r · K). Note the minus sign in Eq. (1) which indicates polar regression in the inertial space. In this model, ε stays constant and where n is the orbital mean motion about the primary, ω the angular rotational frequency, E b the dynamical ellipticity, and e the orbital eccentricity. The dynamical ellipticity expresses degree of non-sphericity of the body, and it is defined using principal moments A ≤ B ≤ C of the inertia tensor as Things become more interesting in the case where the orbital plane about the primary is not constant. As mentioned above, the Colombo top model describes the situation when it performs a uniform precession in the inertial space. In particular, K revolves uniformly on a cone about a fixed direction K in space, such that the orbital inclination I with respect to the reference plane normal to K is constant. The magnitude of the precession rate is μ, though most often the orbital plane performs again regression in the inertial space (assuming I < 90 • ). The interest and complexity of this model revolves about a possibility of a resonance between the two precession frequencies −μ and μ r . In order to describe it using a simple Hamiltonian model, Colombo (1966) observed it is useful to refer r to the reference frame following precession of the orbital plane, thus representing r T = (sin ε cos(h − π/2), sin ε sin(h − π/2), cos ε). Here, h − π/2 is a longitude reckoned from the ascending node of the orbital plane and ε is a colatitude measured from K as above. It is actually an advantage to introduce H = cos ε (and sin ε = √ 1 − H 2 ). This is because in terms of the symplectic variables (h, H ), the Colombo top is a one degree of freedom, conservative problem with a Hamiltonian function and two new nondimensional parameters are defined as As observed by Henrard and Murigande (1987), the discussion can be confined to nonnegative constants a and b thanks to the symmetries (a, h, H ) → (−a, h + π, H ), and Since the equations of motioṅ are singular at H 2 = 1, it is better to use the Cartesian coordinates of the unit momentum vector Then, similarly to Henrard and Murigande (1987), we obtain the Hamiltonian function that generates equations of motioṅ For the sake of minor simplification, let us propose the shifted variables Adding a constant a 2 to the Hamiltonian (8), we simplify it to with the new energy constant Equations of motion for the shifted variables arė Notably, when a = 0, the problem is simplified to the symmetric free top, with constant Z and uniform rotation of r around the third axis, with the frequency √ −2E = |Z |.
The symmetry plane of S2 is z = b, and its vertex line is parallel to the x-axis, passing through y = (E − a 2 )/a. Observe that combining S1 and S2 we can also obtain another surface (see Fig. 1): S3 a paraboloid of revolution Actually, the S3-based function is an alternative Hamiltonian of the Colombo top, leading to the same equations of motion (9) as the Hamiltonian (8).
The paraboloid S3 has the symmetry axis parallel to the z-axis, and passing through the points x = 0, y = −a. Its vertex is located at z = E /b. The advantage of S3 appears when discussing the limit of a = 0. Then, the paraboloid does not change the shape and the intersections of S1 and S3 are circles because of coincidence of the symmetry axes. Contrarily to this, setting a = 0 in S2 results in degeneracy: the cylinder breaks in two parallel planes. On the other hand, b = 0 turns the S3 paraboloid into a (circular) cylinder, whereas S2 retains its shape. Thus, S2 and S3 (or H and H a ) can be considered complementary, although the J. Haponiak et al. regularity at a = 0 seems to us more favorable (e.g., if some perturbation approach is based upon the small inclination assumption).
An integral curve can be represented as a parametric curve in a number of ways. Of course, the best is to solve Eq. (9) and find r(t). But before we accomplish it (and, actually, in order to do it), let us consider parameterizations, where two coordinates are expressed in terms of the third. The most straightforward is to solve the system of S1 and S2 equations, using z as a parameter variable, which leads to (15) All three invariant surfaces S1, S2, S3 intersect the plane x = 0, which is their common plane of symmetry. The integral curves, as the lines of intersection of S1, S2, and S3 also pass through x = 0 and are symmetric with respect to this plane. Moreover, both y and z coordinates of a given integral curve attain their local extremes at x = 0, according to equations of motion (13).

Intersections with the plane x = 0
To find the intersection points of an integral curve with the plane x = 0 for some specified energy E, it is enough to find either y or z coordinate. Knowing one of them, one might recover the other from the relation y 2 + z 2 = 1. But this involves the ambiguity of sign; thus, it is better to find z, and then use the parametric equation (15) for y(z). To benefit from minor simplifications, we first find Z and then Y (Z ).
According to Eq. (15), the relation between X and Z is where W (Z ) is a polynomial of degree 4 Note the absence of the cubic term thanks to the use of the shifted variable Z . Solving the quartic equation W (Z ) = 0 is a tedious task; the details can be found in "Appendix A". Briefly, the four roots Z j are given in terms of the three roots e i of the reduced cubic resolvent equation (additionally modified to match the standard cubic polynomial appearing in the Weierstrass elliptic integrals). Depending on the parameters a and b, selecting some energy value E, the number of meridian intersection points is 4, 3, 2, 1, or 0.
Four intersection points mean that there are two distinct integral curves with the same energy, each intersecting x = 0 in two points: for a lower curve, and for the upper curve, with e j defined in Eq. (141) as functions of E, a, and b, such that Three intersection points mean that either one of the two curves contracts to a single point, or two curves share the same intersection point. These situations are distinguished by the sign of invariant g 3 from Eq. (135).
-If g 3 > 0, then the upper curve shrinks into a point with The remaining intersection points of the lower curve are The expressions of e 1 and e 23 are given in Eq. (145). -If g 3 < 0, then the upper and the lower curves meet at with the remaining intersection points at If e 1 = e 2 = e 3 , then z 2 from Eq. (21) becomes equal to z 34 and only two intersection points remain with a 2 3 + b 2 3 = 1. This case requires a unique value of energy with a simple expression (154).
A more generic situation with two intersection points occurs when E defines only one trajectory, with z 1 < z 2 given by where e 1 and e c are defined in Eqs. (142) and (144). Two other roots of W (z − b) = 0, i.e., z 3 = Z 3 + b and z 4 = z 3 are complex, so they do not define the intersections-see Eq. (151). Finally, one intersection point means that for the energy E there is only one integral curve that contracted to a point with and e 1 given by Eq. (145).
In further discussion, the point of intersection of a regular curve with the x = 0 plane will be named a turning point, unless it is an equilibrium from the dynamical point of view.

Cassini states and homoclinic orbits
Finding the Cassini states can be approached from different points of view. Geometrically, they are the points of tangency of the surfaces S1, S2, and S3. Algebraically, they are the multiple roots of x(z, E) = 0, and x(y, E) = 0. Dynamically, they are the fixed points of equations of motion (9) or, equivalently, the local extremes and saddle points of the Hamiltonian on a sphere S1.
In principle, the multiple roots have been found in Sect. 2.3. But they are given in terms of E, which is an implicit function of a and b as a root of Δ(E, a, b) = 0, where Δ is the discriminant of W (Z ), defined in Eq. (136). Hence, one possible way is to find the roots of Δ = 0, which is a quartic equation in E, and substitute them into z 12 , z 23 , or z 34 from Sect. 2.3 (we already did it for z 234 ). The alternative is to find the stationary points (x, y, z) of the Hamiltonian and use them to evaluate the energy by the substitution into H (x, y, z) = E. The latter path is more convenient and was recently taken by Saillenfest et al. (2019).
Discussing the Cassini states as the multiple roots of W (Z ) = 0, we should include the information, that they are also the roots of its derivative W (Z ) = 0. Since two polynomials have a common root if and only if their resultant equals zero, we ask about the solutions of R(W (Z ), W (Z )) = 0. If we treat W (Z ) and W (Z ) as polynomials in Z , the result is simply the discriminant Δ from Eq. (136), multiplied by a constant (negative) factor. Yet, we can also treat W (Z ) and W (Z ) as polynomials of E with degrees 2 and 1, respectively, whose coefficients depend on Z . So, using we obtain the resultant and For easier reference to Saillenfest et al. (2019), we abandon the shifted variables and solve w(z − b) = 0, which is the same equation as the one they used. However, thanks to considering the Cartesian variables on a sphere, we do not have to warn about any loss of sign, because it does not occur neither in evaluating R(U 1 , U 2 ), nor in obtaining W (Z ). The value of y associated with a given z results directly from the conditionẋ = 0 in Eq. (9), being Solving w(z − b) = 0, we can recycle the procedure used for W (Z ) = 0. More details can be found in "Appendix B" so here we only summarize the final results. When considering the problem on a sphere, there are only two generic situations: either there are two Cassini states C 2 and C 3 , or there are four of them: C 1 , C 2 , C 3 and C 4 . Let us call the former case 'the Colombo top problem of type II,' and the latter-'type IV.' The special case when there are three Cassini states will be called 'the Colombo top problem of type III' (see Fig. 2).
In the following subsections, we briefly review each type, providing the coordinates of the Cassini states C j as the functions of a and b parameters. To avoid confusion with the z i expressions of Sect. 2.3, we label the coordinates of C j as z * j , and y * j (of course, x * j = 0 for all the points). We only provide z * j that allows to find y * j from Eq. (31), and then, the energy E j of the Cassini state is Substituting this energy value into an appropriate z 12 , z 23 , z 34 , or z 234 from Sect. 2.3, should result in returning to the expression of z * j in terms of a and b.

Type II
When a 2 3 + b 2 3 > 1, the dynamics on a sphere is relatively simple. There are two stable equilibria (elliptic points): C 3 at the lower (z < 0) hemisphere and C 2 at the upper hemisphere, located at where and according to Eqs. (164) and (163) We can identify the Cassini states coordinates as z * 3 = z 12 (E 3 ), and z * 2 = z 12 (E 2 ), where z 12 is given by Eq. (26).
The energy values in the type II problem are bounded by E 3 ≤ E ≤ E 2 . All trajectories with energy E 3 < E < E 2 are simple periodic curves (one for each value of E), oscillating in z between z 1 and z 2 , as given by Eq. (25). Extending and modifying the domains labeling of Henrard and Murigande (1987), let us label the entire sphere surface (with the two Cassini states excluded) as D 23 (Fig. 2, left).

Type IV
When a 2 3 + b 2 3 < 1, the flow is shaped by the presence of three elliptic fixed points C 1 , C 2 , C 3 , and a hyperbolic point C 4 (Fig. 2, right). The latter is accompanied by two homoclinic orbits-upper Γ 1 and lower Γ 2 . Thus, the sphere is first partitioned into three domains: D 1 between C 1 and Γ 1 , D 2 between C 2 , Γ 1 , Γ 2 , and the third area bounded by C 3 and Γ 2 . But the last of the three domains is further divided by the curve whose energy equals that of C 1 -the thick dashed curve Γ 3 in Fig. 2. The subdomains D 3 and D 4 may look similar from the geometrical point of view, but note the fact that each curve in D 4 has a companion in D 1 with the same energy, which is not the case in D 3 .
Equation (171) from "Appendix B" can be used as they are, but here we add an alternative form, based upon the transformation (129): where and The unstable equilibrium energy E 4 is of special importance, because it serves to determine the turning points of the separatrices from Eq. (23): z 4 (E 4 ) for Γ 1 , and z 1 (E 4 ) for Γ 2 .
For the energy values E 1 < E < E 4 , when Δ > 0, each E refers to two periodic curves: one in D 4 with the turning points (z 1 , z 2 ) given by (18), and one in D 1 with (z 3 , z 4 ) given by Eq. (19). An energy in E 3 < E < E 1 defines only one periodic trajectory in D 3 , and each E 4 < E < E 2 defines one periodic curve in D 2 . Since Δ < 0 in both cases, the turning points are given by (z 1 , z 2 ) from Eq. (25).

Type III
The case of a 2 3 + b 2 3 = 1 is specific, but it has to be included to understand the bifurcation between the two neighboring types. Increasing a and/or b from the type IV, we observe that the two Cassini states C 1 and C 4 merge into a single point C 14 of neutral stability, the homoclinic orbit Γ 1 contracts to a point, and the separatrix Γ 2 merges with the specific periodic orbit Γ 3 into Γ 23 . Thus, the domains D 1 and D 4 disappear. In course of transition from type III to type II, the division between D 2 and D 3 disappears; hence, we merge them into a single D 23 . Lanchares and Elipe (1995) refer to this transition under the name of teardrop bifurcation.
The coordinates of the Cassini states can be obtained from Eq. (174) and are expressible in terms of a or b alone, resulting in and where The expression of the energy at C 14 is so simple, that we write it explicitly It is useful also in finding the turning point of Γ 23 at x = 0, which is For regular trajectories with E 3 < E < E 2 , and E = E 14 , their turning points are given by z 1 and z 2 from Eq. (25), which also means z 1 < z 2 .

General framework
Equations of motion (9) or (13) admit an exact analytical solution in terms of the elliptic functions. Actually, the solution hinges upon the fact that the variable Z can be considered separately from the remaining two. To see it, we can differentiateŻ = −a X, substituteẊ from the first of Eq. (13), and use the energy integral (11) to eliminate Y , obtaining Equation (44) defines a 1 degree of freedom, conservative system with the potential and the energy integralŻ But the potential V z (Z ) and the energy E z are closely related to the polynomial W (Z ): Thus, whether we use Eq. (46), or the squared Eqs. (13) and (16), the outcome is 4Ż 2 = W (Z ), amenable to the separation of variables method.
Since W (Z ) is a quartic polynomial, finding Z (t) amounts to inverting the elliptic integral where σ z = sgnŻ = −sgnX , and the right-hand side evaluates to The integral to the left of (48) should be a monotonous function of Z to allow the inversion (solving for Z (t)). In the admissible range of Z between two turning points (or one turning point and an unstable equilibrium), the sign ofŻ is constant, which allows to establish its value from the initial condition and to pull σ z out of the integrand. The integral to the left of Eq. (48) is an elliptic integral, and it can be reduced to the Weierstrass normal form by an appropriate transformation Z → s, so that where S(s) = 4s 3 − g 2 s − g 3 is the cubic polynomial of the Weierstrass resolvent (133) for W (Z ) = 0, with the invariants g 2 , g 3 defined in "Appendix A.2"-Eqs. (134) and (135). Note that the initial value Z 0 is always mapped to s → ∞, regardless of the ordering of the integration limits Z 0 and Z , which explains the presence of σ z in the forthcoming transformations.
In all further instances, we will use the abbreviated notation for the Weierstrass ℘ function ℘ (u) = ℘ (u; g 2 , g 3 ), whenever the invariants g 2 and g 3 from "Appendix A.2" are used.
Only the invariants different than g 2 and g 3 will be added to the list of arguments if needed. The derivative of the Weierstrass ℘ function obeys and in the first half-period 0 as expected from Eq. (50). Once the solution for the first half-period is found, its continuation can be investigated by the substitution into Eq. (44) which is free from the restrictions imposed by the inversion procedure.

Initial conditions at turning point
Let us begin with the easiest situation, when the initial condition is Z 0 = Z j = z j − b, and t 0 = t j is the epoch of crossing the turning point X 0 = 0, i.e.,Ż = 0. Then, the integral to the left of (50) is reduced to the Weierstrass normal form in variable s through the rational transformation (Enneper 1890;Bianchi 1901) Replacing the subscript 0 with j in the formulae of Sect. 3.1.1, we find where Thanks to the symmetry of trajectory with respect to the turning point, the solution does not depend on σ z and remains valid for all values of τ j . Given the initial conditions X (t j ) = 0, Y (t j ) = Y j , Z (t j ) = Z j , and knowing Z (t), we may complete the solution for R(t). The two missing variables come from the energy integral (11), with and from the third of the equations of motion (13), i.e., X = −Ż /a, so The coordinates on the unit sphere are, as usually, x = X , y = Y − a, and z = Z + b.

Arbitrary initial conditions: Weierstrass-Biermann form
If the initial conditions are not at the turning point, i.e., Z 0 = Z j , the transformation Z → s is more cumbersome than (54). Whittaker and Watson (1927) quote two alternative forms of the final solution: one due to Weierstrass, published by Biermann (1865), and one by Mordell (1915). Actually, there is yet another, formally elegant form-the secondo metodo d'inversione of Bianchi (1901), but the relation of its constants to initial conditions is rather awkward, so we do not consider it here. The Weierstrass-Biermann form results from the substitution (Enneper 1890) where σ z is the sign of Z − Z 0 , and are the Taylor coefficients in After expressing E in terms of the initial conditions, they take the form Substituting s = ℘ (τ 0 ), and recalling that √ S(s) = −℘ (τ 0 ) over the first half-period, we obtain Then, accounting for Eq. (63) and the fact that √ W 0 = 2a|X 0 | = −2σ z a X 0 , The substitution into Eq. (44) with the initial conditionŻ 0 = −a X 0 shows that the formula (65) is actually valid for all values of τ 0 , so the initial restriction to the first half-period can be abolished. Like before, the remaining two variables Y , X are found from the energy integral and from the equations of motion, respectively. Thus, for Y , we have simply whereas X = −Ż /a, requires the differentiation and some manipulations leading to where the second derivative has been removed using the identity ℘ = 6℘ 2 − (g 2 /2).

Arbitrary initial conditions: Safford form
If X 0 = 0, so τ 0 = τ j , the above solution simplifies to Eq. (58) in a straightforward manner. On the other hand, as pointed out by Safford (1919), the general solution (64) can be derived from the particular solution (55) by assuming τ j = τ 0 + φ j , and making use of the addition theorem for the Weierstrass function ℘. Given the initial conditions X 0 , Y 0 , Z 0 at t = t 0 , we can find the turning point coordinates Z j , Y j for this trajectory using the formulae of Sect. 2.3. Then, according to Safford (1919), the phase φ j is defined through Eq. (60) at Z = Z j , giving s = s j = ℘ (φ j ). Thus, as an alternative to using Eqs. (65) and (66) for arbitrary initial conditions, one can first compute φ j = ℘ −1 (s j ), and then apply the particular solution (58) with τ j = τ 0 + φ j .
But Safford (1919) cared solely about the reduction of the integrand form, so he paid no attention to the problem that ℘ (φ j ) is uniquely invertible only in the domain 0 ≤ φ j ≤ ω 1 , i.e., within the first half-period. In order to properly place the phase in the full period range −ω 1 ≤ φ j ≤ ω 1 , one needs the information about the sign of the derivative ℘ (φ j ). To this end, we take a slightly different approach, that actually leads to a simpler expression for s j . Substituting τ j = φ j , and R = R 0 in Eq. (58), we can solve them to find and this set allows the unique determination of the phase. Then, we can apply Eq. (58) to obtain with Y j and Y derived from the energy integral. Compared to the Weierstrass-Biermann solution, we gain the simplicity at the expense of pre-computing the turning point coordinates and the phase for the given (arbitrary) initial conditions.

Arbitrary initial conditions: Mordell form
The inversion formula of Mordell (1915) was formulated in the language of homogeneous binary forms; in order to apply it to the Colombo top problem, let us translate it to the univariate polynomials framework. The link is simple: the quartic polynomial W (Z ) from Eq. (122), and the quartic binary form can be matched by Thus, after the substitution Z = ξη −1 , Eq. (50) is equivalent to where the definite integral in Z has been replaced by a path-independent line integral. Skipping the intermediate steps described in Mordell (1914Mordell ( , 1915, the inversion of (73) results in where V 0 = V (ξ 0 , η 0 ), andh 0 is the Hessian covariant 1 of V . From the correspondence rules (72), we derive where W n are defined in Eq. (63). By letting Z = ξ/η, and Z 0 = ξ 0 /η 0 , we find that (74) and (75) amount to Proceeding like in Sect. 3.1.3, we obtain the final form looking different from the Weierstrass-Biermann solution (65), yet providing the same values of Z . Notably, the Mordell solution involves only the first power of ℘.
In order to demonstrate the equivalence of the Weierstrass-Biermann and the Mordell solutions, one can multiply the numerator and the denominator in Eq. (65) by the factor τ 0 )), and substitute the expressions of the invariants in terms of W n by analogy with the definitions (134) and (135). The result of this procedure is the Mordell solution (77). Obviously, both (65) and (77) admit the same limit expression (58) when X 0 = 0 and τ 0 = τ j . As usually, the solution for X and Y can be derived from the equations of motion and the energy integral, like in Sect. 3.1.3.

Weierstrass functions in terms of the Jacobi functions
Although the solution in terms of the Weierstrass ℘ function presented in Sect. 3.1 may look universal, its qualitative properties depend on the values of the invariants through the sign of the discriminant Δ. Indeed, the sign plays the central role in expressing the solution in terms of the Jacobian elliptic functions. In this section, only the generic, Δ = 0 cases are to be discussed.
The basic relation between the Weierstrass and Jacobi functions is formally universal (Byrd and Friedman 1971) where But if we restrict considerations to the real arguments and moduli of the elliptic functions, the above expressions are valid only if the discriminant Δ from Eq. (136) is positive.
When Δ < 0, which means complex e 2 and e 3 , the appropriate form is (Abramowitz and Stegun 1972) ℘ (τ ) = = e 1 + γ n cn(u n , k n ) sn(u n , k n )dn(u n , k n ) 2 = e 1 + γ n 1 + cn(2u n , k n ) 1 − cn(2u n , k n ) , where The two cases are linked by the complex modulus transformation-see Byrd and Friedman (1971, formula 165.07). Let the initial conditions at the epoch t 0 be (X 0 , Y 0 , Z 0 ). The equations relating the Jacobi and Weierstrass functions can be substituted in to any of the solution forms provided in Sect. 3.1. For the Weierstrass-Biermann or the Mordell form, it is enough to compute the energy E = E(Y 0 , Z 0 ) and the discriminant Δ to choose the appropriate set (79,80) or (82,83). Then, from the invariants g 2 , g 3 the roots e 1 , e 2 , e 3 are found, which allows the computation of R(t) or r(t) for any epoch t.
Below, we discuss the Safford form from Sect. 3.1.4, which allows to use the initial conditions at any t 0 , but requires the turning point coordinates (X j = 0, Y j = y j + a, Z j = z j − b) as supplementary parameters. Having computed two appropriate turning points, such that either Z 1 < Z 0 < Z 2 , or Z 3 < Z 0 < Z 4 , we pick one of them as the reference point Z j to be used in Eq. (70). Then, after determining phase φ j with respect to the turning point Z j , the motion can be computed for any epoch t, using τ j = τ 0 + φ j .

Motion in D 1 and D 4
The case Δ > 0 occurs only in type IV, when the energy is bounded by E 4 < E < E 1 : either in the domain D 1 , where turning points Z 3 , Z 4 are given by Eq. (19) with Z * 4 < Z 3 < Z 4 , or in the domain D 4 , where the turning points are Z 1 , Z 2 given by Eq. (18) and Z 1 < Z 2 < Z * 4 . Then, selecting any appropriate Z j as the reference point, we can use the expressions where and the phase φ j in can be computed from the incomplete elliptic function of the first kind F: The motion is periodic, and the period P t (with respect to time t) is given by the complete elliptic integral of the first kind K

Motion in D 2 , D 3 , and D 23
The common feature of trajectories in domains D 2 , D 3 , and D 23 is the negative discriminant Δ. Thus, given the initial conditions and resulting energy, we pick Z 1 or Z 2 computed from Eq. (25) as the reference point Z j and compute, for any t where A n = 2(ab + Y j Z j ) The phase in results from The related period is The above formulation is universal, i.e., appropriate in types II, III, and IV, provided Δ < 0.

Reduction rules
By special cases, we mean the ones where elliptic functions reduce to the elementary ones. They could be studied by taking limits of the Jacobi functions at k n or k p tending to 0 or 1. But it is more direct to observe that each of the special cases, be it the Cassini states, separatrices, or the special curve Γ 3 , results from the reduction of the Weierstrass function ℘ (u, g 2 , g 3 ) to the special case Since Δ = 0 implies (g 2 /3) 3 = g 2 3 , reduction to the above form is always possible by the homogeneity relations (Abramowitz and Stegun 1972) For the derivative ℘ , respective equations are obtained by straightforward differentiation. Depending on the sign of g 3 , two procedures are available. For g 3 > 0, the substitution 3 leads straight to ℘ (u; g 2 , g 3 ) = e 1 ℘ ( √ e 1 u; 3, 1), where g 3 = e 3 1 , according to Eq. (145). If g 3 < 0, two steps are taken. First, letting λ = i, we convert ℘ (u; g 2 , g 3 ) = −℘ (i u; g 2 , −g 3 ). (105) Then, with λ = (−g 3 ) − 1 6 , we recover ℘ (u; g 2 , g 3 ) = −2e 12 ℘ (i 2e 12 u; 3, 1), according to Eq. (146).

Cassini states C 2 and C 3
When discussing the Cassini states, we can use the simplest form (58) for Z . Although the time dependence at the Cassini states does vanish due to ab + Y 0 Z 0 = −W (Z 0 )/(8a) = 0, and X 0 = 0, it remains of interest to inspect ℘ (τ ) in the numerator, because its period is the period of small oscillations around the stable equilibrium. Given a and b, one should first establish the problem type in order to compute the appropriate coordinate Z * 2 = z * 2 − b, or Z * 3 = z * 3 − b. These are Eqs. (33), (40), and (36) for the types II, III, and IV, respectively. Starting from this point, the procedure is common: Z * j gives the energy of the Cassini state which substituted in Eq. (134) or (135) gives the invariants g 2 (E j ) or g 3 (E j )-both positive. Resorting to Eqs. (104) and (101), we find that Z in solution (58) depends on the squared tangent of √ 3e 1 /8 t which implies the period where e 1 (E j ) = 3 √ g 3 = √ g 2 /2. The same result can be obtained by taking the limit of (100) at k n → 0. C 4 and homoclinic orbits 0 1 , 0 2 The energy of the unstable Cassini state C 4 can be evaluated from Eq. (107) with j = 4, and Z * 4 given by (36). However, we need it not for the Cassini state itself, but rather to describe the motion on homoclinic orbits having the energy E 4 , which are Γ 1 and Γ 2 . To this end, we will use the Safford form (70) with the reference points given by Eq. (23), namely Z 4 for Γ 1 , and Z 1 for Γ 2 .

Cassini state
Since g 3 < 0, the reduction (106) leads to with the coefficientsÃ where e 12 is given by Eq. (146). The argument whose phase with respect to Z j at the initial epoch t 0 is given by Note that the time rate of argumentũ is the same at both the separatrices, and as expected the solution tends to the Cassini state C 4 asymptotically at t → ±∞.

Cassini state C 1 and special orbit 0 3
The Cassini state C 1 , being a stable equilibrium with g 3 > 0, is characterized by the period of small oscillations given directly by the formula (108) with the energy E 1 evaluated at Y * 1 , Z * 1 deduced from Eq. (36). What makes the difference, compared to C 2 or C 3 , is the presence of another trajectory having the energy E 1 -the special curve Γ 3 .
The motion along Γ 3 can be described using the Safford form solution (70) with respect to the turning points Z 1 or Z 2 , given by Eq. (21). Then, performing the reduction (104), we obtain the solution in terms of trigonometric functions X =Ā γ 3 sin 2ū 1 +B sin 2ū 2 , The argumentū isū where γ 3 = √ 3e 1 /2, as in Eq. (81) for e 3 = −e 1 /2; hence, the time period of the solution is the same as that of small oscillations around C 1 , i.e., P 1 given by Eq. (108). The phase φ j can be computed from The above solution can be obtained either from (86) with k p = 0, and e 3 = −e 1 /2, or-in a different form-from (94) with k n = 0, and g 2 = 3e 2 1 .
The time offset τ 23 is given by

Conclusions
It is common to describe the integral curves of the Colombo top problem as an intersection of a parabolic cylinder and a unit sphere, the latter being centered at the origin of the x, y, z coordinate system. We have proposed a new point of view: the curves can be tracked along the intersection of two out of the three invariant surfaces: a parabolic cylinder, a sphere, and a paraboloid of revolution. This thread has been merely signaled, but it can be of possible interest when designing geometric integrators for the numerical treatment of the Colombo top motion. Moreover, using the shifted coordinates X , Y , Z , one introduces the symmetry to the parabolic cylinder on the paraboloid (at the expense of having an off-centered sphere) which does simplify a number of expressions given in this work. When partitioning the phase space of the Colombo top problem, we have completed the landscape, well known from earlier works, with an interesting but hitherto overlooked feature: the trajectory Γ 3 which is unique by being periodic, yet expressible in terms of elementary functions of time. Its presence calls for the distinction of D 3 and D 4 domains even though qualitatively they look similar. It also adds to a better understanding of the parametric bifurcation associated with the passage from type II, to type IV.
The analytical expressions for the turning points of the Colombo top trajectories as functions of energy, given in Sect. 2.3, had not been reported so far. The expressions for the location of the Cassini states from Sect. 2.4 depend only on parameters a, b. Up to some rearrangement of terms, they are similar to those of Saillenfest et al. (2019) in type II or III. For type IV, when a 2 3 + b 2 3 < 1, the Cardano form provided by Saillenfest et al. (2019) is formally correct, but it gives real values only as the sums of two complex conjugates (casus irreducibilis of the resolvent cubic). In the present work, we have preferred to use expressions based on the purely real trigonometric form whenever the quartic has a positive discriminant.
The differential equation forŻ , with its right-hand side proportional to the square root of the degree 4 polynomial, is not a novelty in celestial mechanics. The same form pops up while discussing the second fundamental model of resonance (Henrard and Lemaitre 1983). Its solution in terms of the Weierstrass elliptic function has always been given either in the simplified form of Eq. (55), as in Ferraz-Mello (2007), or in the Biermann-Weierstrass form (Nesvorný and Vokrouhlický 2016). We have taken an opportunity to recall other possibilities (Safford and Mordell forms) than can be of use in other applications as well.
We hope that the results of the present work will facilitate the study of perturbed Colombo top problems. They should be useful either as the kernel of analytical perturbation procedures, or as a building block of numerical integrators based upon composition methods.
where 3ρ = 1 − a 2 − b 2 , according to Eq. (30). The Descartes rule guarantees that (for nonzero a and b) the resolvent has at least one positive real root to be used in factorization (124).
Before we proceed to solving the resolvent, let us make some important remarks. The four roots of the quartic equation W (Z ) = 0 come in two pairs of the roots of W + (Z ) = 0 and W − (Z ) = 0, i.e., Let ξ 1 be the only, or the greatest positive root of the resolvent (127). Then, from the Vieta's formulas, we find for the remaining two roots ξ 2 + ξ 3 = −2E − ξ 1 , and ξ 2 ξ 3 = a 4 b 2 ξ −1 1 , which allows to see that This leads to the Euler form of the solution Assuming for the real roots 0 < ξ 3 ≤ ξ 2 ≤ ξ 1 , we guarantee a number of properties like the ordering Z 1 ≤ Z 2 ≤ Z 3 ≤ Z 4 , the fact that a given trajectory contains only a pair (Z 1 , Z 2 ) or (Z 3 , Z 4 ), and that if ξ 2 , ξ 3 are complex conjugates, then Z 1 and Z 2 remain real, whereas Z 3 and Z 4 become complex.

A.2 Weierstrass resolvent and its roots
The cubic resolvent equation (127) can be brought to a reduced form without the square term in a number of ways. We choose the substitution based upon the seminvariant (Janson 2011) with Applying it to Eq. (127), and multiplying both sides by 4, we obtain S(s) = 4s 3 − g 2 s − g 3 = 0.
The ordering of roots in (146) is exceptional (e 3 is the greatest), but helps to maintain a coherent notation in further applications. The triple root e 123 = 0 may appear only for g 2 = g 3 = 0, which is possible only when a 2 3 + b 2 3 = 1.

A.3 The roots Z j
Although the roots of W (Z ) are to be expressed in terms of the roots of S(s), we need to include in the discussion not only the invariants Δ, g 2 , and g 3 , but also seminvariants P 2 and Q 2 = 2a 2 0 a 3 − 6a 0 a 1 a 2 + 4a 3 1 = −4a 2 b.
This is due to the fact that although formally it is enough to substitute into (130), the signs of ξ j play a significant role in determining which of the roots are real and which are complex, and there are various ways to create multiple roots. In the following discussion, we will refer to the Theorem 9.3 of Janson (2011), adjusted to the different scaling of our invariants and seminvariant (namely, his P = 48P 2 , Q = 16Q 2 , J = 432g 3 , and I = 12g 2 ).

A.3.1 Four simple real roots (1 > 0 and all j ≥ 0)
If the four real roots exist, they take the form (130) with ξ j defined in Eq. (148) and e j as in (141). This requires not only Δ > 0 to have three simple real roots e j , but also that ξ j ≥ 0 for each j ∈ {1, 2, 3}. The latter is secured by P 2 < 0 and 12P 2 2 − a 2 0 g 2 ≥ 0 (Janson 2011). Substituting (131), (134), and (123), we obtain Δ > 0, and − ρ ≤ 2E 3 < 0, as the condition for the real quadruple Z 1 < Z 2 < Z 3 < Z 4 . If Δ is positive, but the second condition in (149) is not fulfilled, there are no real roots, and Z j form two distinct pairs of complex conjugate numbers.

A.3.3 Multiple roots (1 = 0)
The statement Δ = 0 means only that at least one of the roots is at least a double root. Further distinction is based upon the signs and values of g 2 , g 3 and P 2 . Let us inspect five possibilities involving multiple real roots from the Theorem 9.3 of Janson (2011). A quadruple real root is not possible, because it requires P 2 = g 2 = g 3 = 0, whereas substituting E = 0 we obtain g 3 = 4a 4 b 2 = 0. Two real double roots are also impossible, because they require (among other conditions) that Q 2 = 0, which is not the case. The remaining three cases are the following. 1. A triple real root and one single real root appear when g 2 = g 3 = 0, and P 2 < 0. Taking the resultant of g 2 and g 3 considered as the polynomials in E, one finds that both the invariants admit the common root if the relation well known from Henrard and Murigande (1987). With this constraint, g 2 = 0 can be solved to give a unique negative root According to the statement below Eq. (146), g 2 = g 3 = 0 refers to the triple root e 123 = 0; hence, with ξ 1 = ξ 2 = ξ 3 = − 2 3 E, we obtain where Z 1 is the single, and Z 234 is the triple root. 2. Two simple real roots Z 1 , Z 2 and double real root Z 34 appear when g 2 > 0, P 2 < 0, and 12P 2 2 − g 2 > 0. So, if E is a real root of Δ = 0 in the interval then either Z 1 = − ξ 1 − 2 ξ 23 , Z 2 = − ξ 1 + 2 ξ 23 , Z 34 = ξ 1 , for g 3 > 0,