The Triaxiality Role in the Spin-Orbit Dynamics of a Rigid Body

The variables in which the problem is posed may have a significant impact on its treatment. Our choice is the use of the total angular momentum as the key object to define them, which application for the roto-translatory problem was first introduced in [17] as a result of the application of the elimination of the nodes in the n-body problem [15] to the roto-translatory model. Nevertheless, there is a saying in celestial mechanics that no set of coordinates is good enough. This claim highlights that in every choice of variables, a sacrifice must be done. More precisely, Cartesian variables have a simple formulation, but they do not take advance of the presence of symmetries. Conversely, by using variables referred to the total angular momentum, we incorporate the angles associated to the symmetries allowing for compact expressions and intuitive geometric insight of the relative equilibria. However, this is done at the expenses of having singularities, i. e. a global study of the system requires the use of several charts.

The complete set of canonical variables is (r, φ, ψ, θ, δ, ν, R, Φ, Ψ, Θ, Δ, N). We are not going to provide a complete derivation of them, which may be found in [17, 23]. Instead and with the aim of fixing notation, we provide the geometric meaning of the angles by means of Figure 1 and briefly recall the definition of the canonical angles by following [11]: Let us consider the reference frame S^∗ = (ℓ,n×ℓ, n), where ℓ is the unitary vector defining node of the total angular momentum plane with the horizontal spatial plane and n is the unitary vector pointing in the direction of the total angular momentum. In addition, S^E = {E₁, E₂, E₃} and S^b = {b₁, b₂, b₃} are the spatial and body frames respectively, where b_i corresponds with the principal moment of inertia of ℬ₁. The orientation and center of mass of the body are referred to the new frame by means of (r, φ, ψ, θ, δ, ν), see Figure 1. These angles are determined by the nodes; ℓ_μδ defined by the rotational angular momentum and the spatial plane, ℓ_r = –ℓ_o given by the intersection of the total, rotational and orbital angular momentum planes and ℓ_θ generated by the orbital and spatial planes intersection. Precisely, we have that ϕ=(E1,ℓ^),ψ=(ℓ,ℓI^),θ=(ℓo,r^),δ=(ℓr,ℓI^)andν=(ℓI,b1^).$\begin{array}{} \displaystyle \phi=(\widehat{\mathbf{E}_1,\mathbf{\ell}}), \,\psi=(\widehat{\mathbf{\ell},\mathbf{\ell}_I}), \theta=(\widehat{\mathbf{\ell}_{o},\mathbf{r}}),\, \delta=(\widehat{\mathbf{\ell}_{r},\mathbf{\ell}_I})\,\, {\rm and}\,\, \nu=(\widehat{\mathbf{\ell}_{I},\mathbf{b}_1}). \end{array}$ Moreover, there are three more auxiliary angles which are not among the canonical variables but we will use them later on; λ=(E1,ℓ^),μ=(ℓμδ,ℓI^),σ=(Πr,Πb^)andh=(E1,ℓθ^).$\begin{array}{} \displaystyle \lambda=(\widehat{\mathbf{E}_{1},\mathbf{\ell}}), \mu=(\widehat{\mathbf{\ell}_{\mu\delta},\mathbf{\ell}_I}), \sigma=(\widehat{\Pi_r,\Pi_b})\,\, {\rm and}\, h=(\widehat{\mathbf{E}_{1},\mathbf{\ell_{\theta}}}). \end{array}$

In addition, the conjugate momenta of the variables read as follows

Fig. 1

R,Φ=G⋅E3,Ψ=G⋅n=G,Θ=Go,Δ=Gr,N=Gr⋅b3,$$\begin{array}{} \displaystyle R, \quad \Phi =\mathbf{G} \cdot \mathbf{E}_3,\quad\Psi = \mathbf{G}\cdot \mathbf{n} = G ,\quad \Theta= G_o,\quad \Delta=G_r, \quad N = \mathbf{G}_r \cdot \mathbf{b}_3, \end{array}$$(4)

where G is the total angular momentum vector, G_r is the angular momentum of the secondary body in the body frame and G_o is the orbital angular momentum. Thus, we have the following interpretation of the momenta: (R) Radial velocity of the center of mass. (Φ) Third component of the total angular momentum in space frame. (Ψ) Magnitude of the total angular momentum. (Θ) Magnitude of the angular momentum of the center of mass. (Δ) Magnitude of the angular momentum of the rigid body. (N) Third component of the angular momentum of the rigid body in the body frame (principal axes of inertia).

Finally, we gather below the formulas for the relative positions of the different planes, as functions of the canonical momenta, which are given by

cos⁡(Io+Ir)=Ψ2−Δ2−Θ22ΘΔ,cos⁡Ir=Ψ2+Δ2−Θ22ΨΔ,$$\begin{array}{} \displaystyle \cos (I_o + I_r)=\frac{\Psi^2- \Delta^2- \Theta^2}{2 \Theta \Delta },\qquad \cos I_r=\frac{\Psi^2+ \Delta^2-\Theta^2}{2 \Psi \Delta}, \end{array}$$(5)

cos⁡Io=Ψ2+Θ2−Δ22ΨΘ,cos⁡σ=NΔ,cos⁡I=ΦΨ$$\begin{array}{} \displaystyle \cos I_o =\frac{\Psi^2 + \Theta^2 - \Delta^2}{2 \Psi \Theta},\qquad \cos \sigma=\frac{N}{\Delta},\quad \cos I=\frac{\Phi}{\Psi} \end{array}$$(6)

and the following relations between angles and momenta hold

Ψ=Δcos⁡Ir+Θcos⁡Io,Δsin⁡Io=Θsin⁡Ir=Ψsin⁡(ι),$$\begin{array}{} \displaystyle \Psi=\Delta\cos I_r+\Theta\cos I_o,\qquad \frac{\Delta}{\sin I_o}=\frac{\Theta}{\sin I_r}=\frac{\Psi}{\sin(\iota)}, \end{array}$$(7)

where I is the inclination of the total angular momentum plane Π, I_o and I_r are the inclinations of the orbital and rotational planes with respect to Π and ι = I_o + I_r.

The direction cosines appearing in (3) may be expressed in the body frame by means of the following composition of rotations:

γ=R3(ν)R1(σ)R3(δ)R1(ι)R3(π−θ)e1$$\begin{array}{} \displaystyle \mathbf{\gamma}=R_3(\nu)\,R_1(\sigma)\,R_3(\delta)\,R_1(\iota)\,R_3(\pi-\theta)\;\mathbf{e}_1 \end{array}$$(8)

where γ = (γ₁, γ₂, γ₃) and e₁ = (1, 0, 0). Finally, taking into account that γ12+γ22+γ32=1$\begin{array}{} \displaystyle \gamma_1^2+\gamma_2^2+\gamma_3^2=1 \end{array}$ and after some calculations, we are allowed to express the MacCullagh’s term (3) as follows

U=km32m1r3(2A3−A2−A1)V1+32(A2−A1)V2,$$\begin{array}{} \displaystyle U=\frac{ k m}{32m_1r^3}\left[(2A_3-A_2-A_1)V_1+\frac{3}{2}(A_2-A_1)V_2\right], \end{array}$$(9)

where

V1=−<!-- - -->2(1−<!-- - -->3cι<!-- ? -->2)(1−<!-- - -->3cσ<!-- s -->2)−<!-- - -->3sσ<!-- s -->2[(1−<!-- - -->cι<!-- ? -->)2C2,2,0+(1+cι<!-- ? -->)2C−<!-- - -->2,2,0]−<!-- - -->6sι<!-- ? -->2[sσ<!-- s -->2C0,2,0−<!-- - -->(1−<!-- - -->3cσ<!-- s -->2)C2,0,0]+12cσ<!-- s -->sι<!-- ? -->sσ<!-- s -->[(1−<!-- - -->cι<!-- ? -->)C2,1,0+2cι<!-- ? -->C0,1,0−<!-- - -->(1+cι<!-- ? -->)C−<!-- - -->2,1,0]$$\begin{array}{} \displaystyle V_1\!\! =\!\!\!\!\!& -2(1-3c_\iota^2)(1-3c_\sigma^2) \\[1ex] \displaystyle & -3s_\sigma^2\Big[(1-c_\iota)^2C_{2,2,0}+(1+c_\iota)^2C_{-2,2,0}\Big] \\[1ex] \displaystyle &-6s_\iota^2\Big[s_\sigma^2C_{0,2,0}-(1-3c_\sigma^2)C_{2,0,0}\Big] \\[1ex] \displaystyle & +12c_\sigma s_\iota s_\sigma\Big[(1-c_\iota)C_{2,1,0}+2c_\iota C_{0,1,0}-(1+c_\iota)C_{-2,1,0}\Big] \end{array}$$(10)

which is independent of ν, and V₂, the “triaxiality part” given by

V2=−<!-- - -->(1−<!-- - -->cσ<!-- s -->)2[(1−<!-- - -->cι<!-- ? -->)2C2,2,−<!-- - -->2+(1+cι<!-- ? -->)2C−<!-- - -->2,2,−<!-- - -->2+2sι<!-- ? -->2C0,2,−<!-- - -->2]−<!-- - -->(1+cσ<!-- s -->)2[(1−<!-- - -->cι<!-- ? -->)2C2,2,2+(1+cι<!-- ? -->)2C−<!-- - -->2,2,2+2sι<!-- ? -->2C0,2,2]−<!-- - -->6sι<!-- ? -->2sσ<!-- s -->2[C2,0,2+C2,0,−<!-- - -->2]+4sσ<!-- s -->2(1−<!-- - -->3cι<!-- ? -->2)C0,0,2+4sι<!-- ? -->sσ<!-- s -->(1−<!-- - -->cσ<!-- s -->)[(1−<!-- - -->cι<!-- ? -->)C2,1,−<!-- - -->2+2cι<!-- ? -->C0,1,−<!-- - -->2−<!-- - -->(1+cι<!-- ? -->)C−<!-- - -->2,1,−<!-- - -->2]+4sι<!-- ? -->sσ<!-- s -->(1+cσ<!-- s -->)[−<!-- - -->(1−<!-- - -->cι<!-- ? -->)C2,1,2−<!-- - -->2cι<!-- ? -->C0,1,2+(1+cι<!-- ? -->)C−<!-- - -->2,1,2],$$\begin{array}{} \displaystyle V_2 =\!\!\!\!& -(1-c_\sigma)^2\Big[(1-c_\iota)^2C_{2,2,-2}+(1+c_\iota)^2C_{-2,2,-2}+2s_\iota^2C_{0,2,-2}\Big] \\[1ex] \displaystyle & -(1+c_\sigma)^2\Big[(1-c_\iota)^2C_{2,2,2}+(1+c_\iota)^2C_{-2,2,2}+2s_\iota^2C_{0,2,2}\Big] \\[1ex] \displaystyle &-6s_\iota^2s_\sigma^2\Big[C_{2,0,2}+C_{2,0,-2}\Big]+4s_\sigma^2(1-3c_\iota^2)C_{0,0,2} \\[1ex] \displaystyle & +4s_\iota s_\sigma(1-c_\sigma)\Big[(1-c_\iota)C_{2,1,-2}+2c_\iota C_{0,1,-2}-(1+c_\iota)C_{-2,1,-2}\Big] \\[1ex] \displaystyle & +4s_\iota s_\sigma(1+c_\sigma)\Big[-(1-c_\iota)C_{2,1,2}-2c_\iota C_{0,1,2}+(1+c_\iota)C_{-2,1,2}\Big], \end{array}$$(11)

and the notation has been abbreviated by writing C_{i, j, k} ≡ cos(iθ + jδ + kν) and c_x ≡ cosx and s_x ≡ sinx.

Facing a non-integrable Hamiltonian system requires the development of a perturbation theory. A usual way to proceed is to expand the Hamiltonian function in power series and truncate it at a certain order. This procedure is in general an uncontrolled approximation since, for most cases, we are no longer sure that the solutions of the truncated model stay close to the full model. In the case of the FG2BP, it is customary to consider the orbital and rotational kinetic energies and choose the first term in the expansion of the gravitational potential, see (1). Thus, we end up with the free rigid body Hamiltonian (the rotational kinetic energy) and the Kepler system Hamiltonian. That is, the orbital kinetic energy plus the first non-zero term of the gravitational potential. However, there is no any theorem claiming that this is the only way it can be done. The novelty of an intermediary model is that it allows to consider truncation of partial order in the power series expansion. In other words, we pick some “entire” terms of the power series plus a “piece” of one term. The only requirement that must be fulfilled in this procedure is that the intermediary model should be searchable.

In [11], the authors proposed an axis-symmetric integrable intermediary model, whose accuracy was tested by comparing with the MacCullagh’s truncation and showing a good performance in the numerical experiments. Here, we complete this previous study by investigating the triaxial case. One of our aims is to analyze the physical-parametric bifurcations of relative equilibria due to the elimination of the axial symmetry of the body. Keeping this motivation in mind, we propose our intermediary following exactly the same procedure than in [11], except for the triaxial parameter. That is to say, we only take into account the first line of V₁ in (10). Then, let the perturbing potential of the intermediary model be given by

V=−km16m1r3(2A3−A2−A1)(1−3cι2)(1−3cσ2),$$\begin{array}{} \displaystyle \mathcal{V}=-\dfrac{km}{16m_1r^3}(2A_3-A_2-A_1)(1-3c_\iota^2)(1-3c_\sigma^2), \end{array}$$(12)

which leads us to the final expression of the intermediary Hamiltonian

H=12R2+Θ<!-- T -->2r2−<!-- - -->κ<!-- ? -->r+q2sin2⁡<!-- ? -->(ν<!-- ? -->)A1+cos2⁡<!-- ? -->(ν<!-- ? -->)A2(Δ<!-- ? -->2−<!-- - -->N2)+1CN2−<!-- - -->κ<!-- ? -->(2A3−<!-- - -->A2−<!-- - -->A1)16r31−<!-- - -->3cι<!-- ? -->21−<!-- - -->3cσ<!-- s -->2,$$\begin{array}{} \displaystyle \mathcal{H}\!\!\!\!&\displaystyle =\frac{1}{2}\left(R^2+\frac{\Theta^2}{r^2}\right)-\frac{\kappa}{r}+\frac{q}{2}\left[\left(\frac{ \sin ^2(\nu )}{A_1}+\frac{ \cos ^2(\nu )}{A_2}\right)(\Delta^2-N^2)+\frac{1}{C}N^2\right]\\ &\displaystyle -\frac{\kappa (2A_3-A_2-A_1)}{16\,r^3}\left(1-3 c^2_\iota\right)\left(1-3 c^2_\sigma\right), \end{array}$$(13)

where q = m/m₁. Furthermore, with the aim of alleviate formulas, we have considered the Hamiltonian per unit of mass by scaling the system and inertia momenta as follows:

H′=H/m;R′=R/m;Θ′=Θ/m;Δ′=Δ/m;N′=N/m;Ψ′=Ψ/m;Φ′=Φ/m;A1′=A1/m1;A2′=A2/m1;A3′=A3/m1.$$\begin{array}{} \displaystyle \mathcal{H}'\,\,=\mathcal{H}/m;\quad \,\, R'=R/m;\quad \Theta'=\Theta/m;\quad \Delta'=\Delta/m;\quad N'=N/m;\quad \Psi'=\Psi/m;\\ \displaystyle \,\,\Phi'\,=\Phi/m;\quad A_1'=A_1/m_1;\quad A_2'=A_2/m_1;\quad A_3'=A_3/m_1. \end{array}$$(14)

Nevertheless, for the sake of simplicity, we keep the original notation without primes on the variables. Then, the 2-DOF Hamiltonian system of differential equations associated with (13) is given by the following expressions:

r˙=R$$\begin{array}{} \displaystyle \dot{r}&=&R\\ \end{array}$$(15)

R˙=Θ2r3−κ1r2−3ακ16r41−3cι21−3cσ2$$\begin{array}{} \displaystyle \dot{R}&=&\frac{\Theta ^2}{r^3}-\kappa \frac{1}{r^2}-\dfrac{3 \alpha \kappa}{16 r^4}\left(1-3c_\iota^2\right) \left(1-3c_\sigma^2\right) \end{array}$$(16)

θ˙=Θr2−3ακ8Θr3cι1−3cσ2(cι+ΘΔ)$$\begin{array}{} \displaystyle \dot{\theta }&=&\frac{\Theta }{ r^2}-\dfrac{3 \alpha \kappa}{8 \Thetar^3}\, c_\iota \left(1-3c_\sigma^2\right) (c_\iota+\dfrac{\Theta}{\Delta} ) \end{array}$$(17)

ψ˙=3ακΨ8ΔΘr3cι1−3cσ2$$\begin{array}{} \displaystyle \dot{\psi }&=&\dfrac{3 \alpha \kappa \Psi}{8\Delta \Theta r^3}c_\iota \left(1-3c_\sigma^2\right) \end{array}$$(18)

δ˙=qsin2⁡(ν)A1+cos2⁡(ν)A2−3ακ8Δ2r31−3cι2cσ2+cι1−3cσ2(cι+ΔΘ)Δ$$\begin{array}{} \displaystyle \dot{\delta }&=& \left[q \left(\frac{ \sin ^2(\nu )}{A_1}+\frac{ \cos ^2(\nu )}{A_2}\right)-\dfrac{3 \alpha \kappa}{8 \Delta^2 r^3}\left(\left(1-3c_\iota^2\right)c_\sigma^2+ c_\iota \left(1-3c_\sigma^2\right)(c_\iota+\dfrac{\Delta}{\Theta} )\right)\right]\Delta \end{array}$$(19)

ν˙=qA3−qsin2⁡(ν)A1+cos2⁡(ν)A2+3ακ8Δ2r31−3cι2N$$\begin{array}{} \displaystyle \dot{\nu}&=& \left[ \frac{q}{A_3} -q\left(\frac{\sin ^2(\nu )}{A_1}+\frac{\cos ^2(\nu )}{A_2}\right) +\dfrac{3 \alpha\kappa}{8\Delta ^2 r^3}\left(1-3c_\iota^2\right)\right]N \end{array}$$(20)

N˙=q(A1−A2)2A1A2(Δ2−N2)sin⁡(2ν)$$\begin{array}{} \displaystyle \dot{N}&=&q\,\frac{(A_1-A_2) }{2\,A_1A_2}\,(\Delta^2 -N^2) \,\sin (2\nu ) \end{array}$$(21)

together with the integrals Θ̇ = φ̇ = Φ̇ = Ψ̇ = Δ̇ = 0.

Note that, in general, a 2-DOF system is not integrable. Thus, in the triaxial case, the analytical integration is not provided and the integrability of the system remains as an open question, which is not in the scope of the present paper.

In this part we analyze the impact of the triaxiality parameter ρ. For this purpose we fix eccentricity to zero in all the evaluations and consider four different triaxialities by modifying the secondary body principal axes. Precisely we study the cases ρ ∈ {0, 0.0005, 0.15, 0.35}, which measure the evolution of the model when ρ ranges from the axial-symmetric case to a moderate triaxiality. In Figure 3 it is shown how the performance becomes poorer as the triaxiality parameter grows. More in detail, we observe for the radial variable that differences duplicate in each step from ρ = 0 to ρ = 0.35. Analogously, slow angles ψ and θ reduce their accuracy. However, in this case there is a drastic change from the axial-symmetric case to the first step into the triaxial model, which increases differences by a factor of ten. Moreover, there is also a significant difference in the behavior of these angles, since the differences with the full model are now strict monotone functions instead of having an oscillatory nature as in the axial-symmetric case. Rather than the expected and moderate degradation in performance of the radial and the slow-angles, the main novelty in the triaxial case with respect to the axial symmetric one is that their approximation deteriorates more and much faster. For that reason and for the sake of brevity, fast angles will not be displayed in the following sections, since their behavior is quite similar.

Fig. 2

Fig. 3

Before we present our experiments in the general case, we study in detail the transition from ρ = 0 to a very low triaxiality. In this regard, Figure 2 shows how fast angles behaves for a nearly zero triaxiality. It is observed that the annunciated degradation in accuracy depends gradually on the triaxiality parameter ρ, though it is presented very early. However, our experiments are given for fast rotating bodies. As such, an analysis for the case of slow rotating bodies is required in order to find out how the rotational angular momentum modulus influences in accuracy.

Now we consider the case of higher triaxialities. Despite of what have been said, these experiments still show a very high accuracy for moderate values of ρ. Namely, in all the simulations provided in Figure 3, the error ratio (r_Interm–r_full)/r₀ is under 0.00064% after one orbital period and for the slow angles we obtain differences up to the order 10^–6 radians. Fast angles are no longer well approximated.

Initial conditions for simulations in Figure 2 and Figure 3.

r0=1550,ϕ0=2π9,ψ0=π4,θ0=0,σ0=π2,ν0=0,μ0=7π18,R0=0,Io=π18,Ir=π12,I=7π36,e=0,$$\begin{array}{} \displaystyle &r_0=1550,\; \phi_0=\dfrac{2\pi}{9},\,\; \psi_0=\dfrac{\pi}{4},\,\; \theta_0=0,\;\, \sigma_0=\dfrac{\pi}{2},\; \nu_0=0,\; \\ \displaystyle &\mu_0=\dfrac{7\pi}{18},\;R_0=0,\, I_o=\dfrac{\pi}{18},\; I_r=\dfrac{\pi}{12},\, I=\dfrac{7\pi}{36},\; e=0, \end{array}$$(23)

making use of the formulas in the Appendix, we are led to the following initial values for the canonical variables

r0=3.12⋅Rp,ϕ0=0.698131,ψ0=0.785398,θ0=0,δ0=0.630221,ν0=0,R0=0,Φ0=0.125460,Ψ0=0.153158,Θ0=0.093797,Δ0=0.062930,N0=0.$$\begin{array}{} \displaystyle &r_0=3.12\cdot R_p,\;\;\;\, \phi_0=0.698131,\;\;\;\, \psi_0=0.785398,\;\;\;\, \theta_0=0,\;\;\;\, \delta_0=0.630221,\;\;\;\, \nu_0=0,\; \\ \displaystyle &R_0=0,\;\; \Phi_0=0.125460,\;\; \Psi_0=0.153158,\;\; \Theta_0=0.093797,\;\; \Delta_0=0.062930,\;\;N_0=0. \end{array}$$(24)

In the formulation of the triaxial model in Section 2 we restrict the validity of the model to small eccentricities. That is to say, we consider eccentricity up to 0.5, since for higher values even the full model (MacCullagh’s approximation) does not provide a valuable estimation. In Figure 4e is set to 0.1 and 0.5, showing that accuracy decreases by factor two as eccentricity grows and leading to a similar scenario than in the triaxiality analysis.

Fig. 4

Next, we specify the initial conditions and parameters to assess the role of the eccentricity. Firstly, we fix the dimensions of the secondary body to be 60, 55 and 45 Km. Then, we obtain the following associated moments of inertia

A1=1.06⋅1021,A2=1.18⋅1021,A3=1.28⋅1021.$$\begin{array}{} \displaystyle A_1=1.06\cdot 10^{21}, \;A_2=1.18\cdot 10^{21}, \;A_3=1.28\cdot 10^{21}. \end{array}$$

The remaining variables, inclinations and momenta are also fixes with the exception of the eccectricity

r0=1550,ϕ0=2π9,ψ0=π4,θ0=0,σ0=π2,ν0=0,μ0=7π18,R0=0,Io=π18,Ir=π12,I=7π36,e,$$\begin{array}{} \displaystyle &r_0=1550,\; \phi_0=\dfrac{2\pi}{9},\; \psi_0=\dfrac{\pi}{4},\; \theta_0=0,\; \sigma_0=\dfrac{\pi}{2},\; \nu_0=0,\; \\ \displaystyle &\mu_0=\dfrac{7\pi}{18},\;R_0=0,\; I_o=\dfrac{\pi}{18},\; I_r=\dfrac{\pi}{12},\; I=\dfrac{7\pi}{36},\; e, \end{array}$$(25)

making use of the formulas in the Appendix, we are led to the corresponding initial values for each value of the eccentricity. Namely, for e = 0.1 we obtain

r0=3.12⋅Rp,ϕ0=0.698131,ψ0=0.785398,θ0=0,δ0=0.630221,ν0=0,R0=0,Φ0=0.124829,Ψ0=0.152388,Θ0=0.09332,Δ0=0.062614,N0=0,$$\begin{array}{} \displaystyle &r_0=3.12\cdot R_p,\;\;\;\, \phi_0=0.698131,\;\;\;\, \psi_0=0.785398,\;\;\;\, \theta_0=0,\;\;\;\, \delta_0=0.630221,\;\;\;\, \nu_0=0,\; \\ \displaystyle &R_0=0,\;\; \Phi_0=0.124829,\;\; \Psi_0=0.152388,\;\; \Theta_0=0.09332,\;\; \Delta_0=0.062614,\;\; N_0=0, \end{array}$$(26)

and for e = 0.5

r0=3.12⋅Rp,ϕ0=0.698131,ψ0=0.785398,θ0=0,δ0=0.630221,ν0=0,R0=0,Φ0=0.108650,Ψ0=0.132637,Θ0=0.081229,Δ0=0.054498,N0=0.$$\begin{array}{} \displaystyle r_0=3.12\cdot R_p,\;\;\;\, \phi_0=0.698131,\;\;\;\, \psi_0=0.785398,\;\;\;\, \theta_0=0,\;\;\;\, \delta_0=0.630221,\;\;\;\, \nu_0=0,\; \\ \displaystyle R_0=0,\;\; \Phi_0=0.108650,\;\; \Psi_0=0.132637,\;\; \Theta_0=0.081229,\;\; \Delta_0=0.054498,\;\; N_0=0. \end{array}$$(27)

This section present a representative selection for a low range of altitudes. In Figure 5 we compare the results obtained for 1000 and 2000 Km, with fixed e = 0.3, showing that a different picture emerges for the height analysis. Namely, accuracy is enhanced, at least duplicated, with height. However, a more exhaustive study considering a wider range of altitudes should be done to clarify at what extend is this claim true.

Fig. 5

Initial conditions and parameters are set as follows

r0=1550,ϕ0=2π9,ψ0=π4,θ0=0,σ0=π2,ν0=0,μ0=7π18,R0=0,Io=π18,Ir=π12,I=7π36,e=0.3,$$\begin{array}{} \displaystyle r_0=1550,\;\, \phi_0=\dfrac{2\pi}{9},\;\, \psi_0=\dfrac{\pi}{4},\;\, \theta_0=0,\; \sigma_0=\dfrac{\pi}{2},\;\, \nu_0=0,\; \\ \displaystyle \mu_0=\dfrac{7\pi}{18},\;R_0=0,\; I_o=\dfrac{\pi}{18},\; I_r=\dfrac{\pi}{12},\; I=\dfrac{7\pi}{36},\; e=0.3, \end{array}$$(28)

Then, we obtain the following initial conditions for heights H = 1000 and H = 2000 Km accordingly

r0=3.12⋅Rp,ϕ0=0.698131,ψ0=0.785398,θ0=0,δ0=0.630221,ν0=0,R0=0,Φ0=0.119679,Ψ0=0.146101,Θ0=0.0894754,Δ0=0.060031,N0=0,$$\begin{array}{} \displaystyle r_0=3.12\cdot R_p,\;\;\;\, \phi_0=0.698131,\;\;\;\, \psi_0=0.785398,\;\;\;\, \theta_0=0,\;\;\;\, \delta_0=0.630221,\;\;\;\, \nu_0=0,\; \\ \displaystyle R_0=0,\;\; \Phi_0=0.119679,\;\; \Psi_0=0.146101,\;\; \Theta_0=0.0894754,\;\; \Delta_0=0.060031,\;\; N_0=0, \end{array}$$(29)

for height = 1000Km and

r0=5.12⋅Rp,ϕ0=0.698131,ψ0=0.187160,θ0=0,δ0=0.630221,ν0=0,R0=0,Φ0=0.153312,Ψ0=0.146101,Θ0=0.114620,Δ0=0.076901,N0=0,$$\begin{array}{} \displaystyle r_0=5.12\cdot R_p,\;\;\;\, \phi_0=0.698131,\;\;\;\, \psi_0=0.187160,\;\;\;\, \theta_0=0,\;\;\;\, \delta_0=0.630221,\;\;\;\, \nu_0=0,\; \\ \displaystyle R_0=0,\;\; \Phi_0=0.153312,\;\; \Psi_0=0.146101,\;\; \Theta_0=0.114620,\;\; \Delta_0=0.076901,\;\; N_0=0, \end{array}$$(30)

for height = 2000Km.

Along this subsection we study the equilibria obtained under a fixed inclination of the dihedral angle ι which determines the relative position of the orbital and rotational planes. In particular we consider (1−3cι2)=0$\begin{array}{} \displaystyle (1-3c_\iota^2)=0 \end{array}$ , i.e. cι=±1/3,$\begin{array}{} \displaystyle c_\iota=\pm \sqrt{1/3}, \end{array}$ inclination dubbed as critical inclination. Observe that for this critical inclination the perturbative term in (16) and (20) vanishes and these two equations correspond with the decoupled motion of the Kepler and the free rigid body problems respectively.

Invariant 4-tori 𝕋⁴(θ, ψ, δ, ν).

Our first approach to tackle the problem of searching for orbits of constant radius within our intermediary is by examining the subsystem (ṙ, Ṙ), i.e. equations (15) and (16). Just by being (1−3cι2)=0$\begin{array}{} \displaystyle (1-3c_\iota^2)=0 \end{array}$ we get an equilibrium of (ṙ, Ṙ) which provide us the following expression for the constant radius:

r=Θ2κ.$$\begin{array}{} \displaystyle r=\frac{\Theta ^2}{ \kappa}. \end{array}$$(32)

The orbits obtained under this inclination are restricted to be into the invariant 𝕋⁴(θ, ψ, δ, ν).

Following the classification previously indicated, we study here the equilibria obtained under a fixed inclination of the σ angle. In particular, we called inclined equilibria to those for which c_σ ≠ 0, i.e. N ≠ 0. Due to the factorization obtained of (20) and in order to obtain orbits of constant radius for this family of equilibria we have:

qA3−qAi+3ακ81−3cι2Δ2r3=0,i∈{1,2}.$$\begin{array}{} \displaystyle \frac{q}{A_3} -\frac{q}{A_i} +\dfrac{3 \alpha\kappa}{8}\frac{\left(1-3c_\iota^2\right)}{ \Delta ^2 r^3}=0, \quad i\in \{1,2\}. \end{array}$$(33)

Note that the above formula implies (1−3cι2)>0$\begin{array}{} \displaystyle (1-3c_\iota^2)>0 \end{array}$ and the following restriction for Δ

Δ2≤3ακ8qr3A3AiA3−Ai.$$\begin{array}{} \displaystyle \Delta ^2 \leq\dfrac{3 \alpha\kappa}{8q r^3 }\left(\dfrac{A_3A_i}{A_3-A_i}\right). \end{array}$$(34)

As result of the appearance of the r cubed constrain in the denominator we have that these equilibria are just possible for very small values of Δ, which implies bodies in a slow rotational regime.

Fixed ν. Invariant 3-tori 𝕋³(θ, ψ, δ).

Considering the subsystems (ṙ, Ṙ) and (ν̇, Ṅ) under a inclined equilibria and 1−3cσ2≠0$\begin{array}{} \displaystyle 1-3c_\sigma^2\neq0 \end{array}$ along with some particular initial value problem conditions, i.e. R = 0, ν = ν_i, N ≠ 0, 1−3cσ2≠0$\begin{array}{} \displaystyle 1-3c_\sigma^2\neq0 \end{array}$ , from (16) and (20) we obtain an expression for the radius

r=Θ22κ+Θ44κ2−3α(1−3cι2)(1−3cσ2)16,$$\begin{array}{} \displaystyle r=\frac{\Theta ^2}{2 \kappa}+\sqrt{\frac{\Theta ^4}{4 \kappa^2}-\frac{3\alpha(1-3c_\iota^2) (1-3 c_\sigma ^2)}{16}}, \end{array}$$(35)

which leads to an implicit relation between c_ι and c_σ which provides us equilibrium that guarantees us the existence of orbits of constant radius within the invariant 3-tori 𝕋³(θ, ψ, δ)

3ακ1−3cι2AiA38Δ2q(A3−Ai)3=Θ44κ2−316α1−3cι21−3cσ2+Θ22κ.$$\begin{array}{} \displaystyle \sqrt[3]{\frac{3 \alpha \kappa \left( 1-3c_\iota^2\right) \text{A}_i \text{A}_3}{8 \Delta ^2 q (\text{A}_3-\text{A}_i)}}=\sqrt{\frac{\Theta ^4}{4 \kappa ^2}-\frac{3}{16} \alpha \left(1-3 c_\iota^2\right) \left(1-3 c_\sigma^2\right)}+\frac{\Theta ^2}{2 \kappa }. \end{array}$$(36)

Using Taylor series of the previous expression of c_ι about 0 we obtain the zero order approximation of c_σ

cσ=−16aΘ2+16a2κ+3ακ3ακ,$$\begin{array}{} \displaystyle c_\sigma=-\dfrac{\sqrt{16 a \Theta^2 +16 a^2 \kappa + 3 \alpha \kappa}}{3\sqrt{\alpha \kappa}}, \end{array}$$

where

a=3αA1A3κ8Δ2q(A1−A3)3.$$\begin{array}{} \displaystyle a=\sqrt[3]{\frac{3 \alpha \text{A}_1 \text{A}_3 \kappa }{8 \Delta ^2 q (\text{A}_1-\text{A}_3)}}. \end{array}$$

This expression shows us how the triaxiality of B₁ conects the relation between c_σ and c_ι.

Searching for conditions that lead us to find equilibria for δ and ψ within our 3-tori we notice analyzing (18) and (19) that there is no a simultaneous equilibrium for δ and ψ. Therefore, we continue studying the existence of 2-tori originated by fixing one of the remaining angles.

The variation of the remainder angles for these orbits are:

θ˙=nθ=Θr2−3ακ8cι(1−3cι2)(cι+ΘΔ)Θr3,$$\begin{array}{} \displaystyle \dot{\theta } & = & n_{\theta} = \frac{\Theta }{ r^2}-\dfrac{3 \alpha \kappa}{8}\frac{c_\iota(1-3c_\iota^2)(c_\iota+\dfrac{\Theta}{\Delta} )}{\Thetar^3}, \end{array}$$(37)

ψ˙=nψ=3ακ8cι(1−3cι2)ΨΔΘr3.$$\begin{array}{} \displaystyle \dot{\psi }&=&n_\psi=\dfrac{3 \alpha \kappa}{8}\frac{ c_\iota (1-3c_\iota^2) \Psi }{ \Delta \Theta r^3}. \end{array}$$(38)

Under this particular body-inclined equilibria the relations between inclinations are given by:

cσ12=1−(4−τ1)cι2−(1−τ1)ζcι(1−τ1)(1−6cι2−3ζcι) for νii∈{1,3},$$\begin{array}{} \displaystyle c_{\sigma 1}^2=\dfrac{1-(4-\tau_1)c_\iota^2-(1-\tau_1)\zeta c_\iota}{(1-\tau_1)(1-6c_\iota^2-3\zeta c_\iota)} \quad \text{ for } \nu_i \quad i \in \{ 1,3 \}, \end{array}$$

cσ22=1−(4−τ2)cι2−(1−τ2)ζcι(1−τ2)(1−6cι2−3ζcι) for νii∈{2,4},$$\begin{array}{} \displaystyle c_{\sigma_2}^2=\dfrac{1-(4-\tau_2)c_\iota^2-(1-\tau_2)\zeta c_\iota}{(1-\tau_2)(1-6c_\iota^2-3\zeta c_\iota)} \quad \text{ for } \nu_i \quad i \in \{ 2,4 \}, \end{array}$$(39)

where we introduce the following notation

τ1=A1A3,τ2=A2A3,ζ=ΔΘ.$$\begin{array}{} \displaystyle \tau_1=\dfrac{A_1}{A_3}, \quad \quad \tau_2=\dfrac{A_2}{A_3}, \quad \quad \zeta=\dfrac{\Delta}{\Theta}. \end{array}$$

These relationships lead us to the following family of curves

cσ1=cσ1(cι;τ1),cσ2=cσ2(cι;τ2)$$\begin{array}{} \displaystyle c_{\sigma_1}=c_{\sigma_1}(c_\iota;\tau_1), \quad \quad c_{\sigma_2}=c_{\sigma_2}(c_\iota;\tau_2) \end{array}$$

whose graph is studied in the Cartesian plane Oc_ιc_σ.

Notice that the triaxiality of B₁ is encapsulated in the ratios τ₁ and τ₂ respecting some conditions as 0 < τ₁ <τ₂ < 1 and 1 < τ₁ + τ₂. On figure 6 we distinguish curves of equilibria. A complete analysis of these curves is a work on progress in [6].

Fig. 6

A partial study of the triaxiality role is presented now in more detail due to the important effects produced on the inclinations of c_σ and c_ι for this particular equilibria. Some important restrictions for the sake of simplicity must be considered. We restrict the domain of c_ι to 0<cι<1/3$\begin{array}{} \displaystyle 0\lt c_\iota\lt 1/\sqrt{3} \end{array}$ and use the following nomenclature:

c_ι^∗ is the c_ι-intercept value of the curve c_σ1, i.e. 0 = c_σ1(c_ι^∗; τ₁) and c_ι^∗∗ is the c_ι-intercept value of the curve c_σ2 respectively, namely

cι∗=ζ(1−τ1)−ζ2(1−τ1)2−4(τ1−4)2(τ1−4),cι∗∗=ζ(1−τ2)−ζ2(1−τ2)2−4(τ2−4)2(τ2−4).$$\begin{array}{} \displaystyle {c_\iota}^*=\frac{\zeta(1-\tau_1) -\sqrt{\zeta^2(1 -\tau_1)^2-4 (\tau_1 -4 )} }{2 (\tau_1 -4)}, \quad {c_\iota}^{**}=\frac{\zeta(1-\tau_2) -\sqrt{\zeta^2(1 -\tau_2)^2-4 (\tau_2 -4 )} }{2 (\tau_2 -4)}. \end{array}$$

Within our restriction and fixing a certain inclination c̃_ι, we distinguish three regions:

–

c̃_ι ∈ (0, c_ι^∗): No equilibria is found for this fixed inclination.

–

c̃_ι ∈ (c_ι^∗, c_ι^∗∗): One equilibrium found. As we notice in figures 7 and 8 (central) for a fixed inclination on this region we find one c_σ value that guarantee us an equilibrium.

Fig. 7

Fig. 8

Moreover in figures 7 and 8 (left) we included a representation of the triaxiality region.

As we observe the triaxial shape of B₁ acquires a fundamental role in the amplitude of the different regions of equilibria. Notice that the mentioned amplitude is also affected by the ratio ζ. From a more detail study we conclude that this situation is robust enough even for values of ζ > 1.

Fixed ν and ψ. Invariant 2-tori 𝕋²(θ, δ) and periodic orbits.

Working on the subsystems (ṙ, Ṙ), (ν̇, Ṅ) and (ψ̇, Ψ̇) under an inclined equilibria and (1−3cσ2)=0$\begin{array}{} \displaystyle (1-3c_\sigma^2)=0 \end{array}$ along with the following initial value problem, we obtain an equilibrium which guarantee us the existence of an orbit of constant radius within the invariant 2-tori 𝕋²(θ, δ),

R=0,r=Θ2κ,ν=νi,1−3cσ2=0,cι2=131−qA3−AiA3Ai8Δ2Θ63ακ4.$$\begin{array}{} \displaystyle R=0,\quad r=\frac{\Theta ^2}{\kappa}, \quad \nu=\nu_i, \quad\\ \displaystyle 1-3c_\sigma^2=0, \quad c_\iota^2=\dfrac{1}{3}\left[1-q\left(\dfrac{A_3-A_i}{A_3A_i}\right)\dfrac{8\Delta^2\Theta^6 }{3\alpha\kappa^4}\right]. \end{array}$$(40)

This particular inclination of σ provides us a body-inclined equilibria and leads us to recover in (16) the Kepler problem formulation. Moreover, it also guarantees us the existence of orbits within 𝕋²(θ, δ).

Particularly the variation of the remainder angles are given by:

θ˙=nθ=Θr2,$$\begin{array}{} \displaystyle \dot{\theta}&=n_\theta=&\dfrac{\Theta}{r^2}, \end{array}$$(41)

δ˙=nδ=qΔ2A3+Ai3A3Ai.$$\begin{array}{} \displaystyle \dot{\delta}&=n_\delta=&q\Delta\left(\dfrac{2A_3+A_i}{3A_3 A_i }\right). \end{array}$$(42)

Notice that these mean motions for θ and δ provide us periodic orbits within 𝕋²(θ, δ).

Along this section we study the equilibria obtained under a fixed inclination of the σ angle which determine the relative position of the rotational and body planes. In particular we consider c_σ = 0 inclination dubbed as body-perpendicular equilibria. It is straightforward that N = 0 is equivalent to have a perpendicular inclination. Observe that some of the perpendicular equilibria study on this section present critical inclination. Note that not likewise the inclined equilibria, the body-perpendicular equilibria allow a free regime of Δ.

Fixed ν .Invariant 3-tori 𝕋³(θ, ψ, δ).

Examining the subsystems (ṙ, Ṙ) and (ν̇, Ṅ) under a perpendicular inclination c_σ = 0 along with the following initial value problem we obtain an equilibrium which guarantees us the existence of an orbit of constant radius within the invariant 3-tori 𝕋³(θ, ψ, δ)

R=0,r=Θ22κ+Θ44κ2−3α(1−3cι2)16,ν=νi,N=0.$$\begin{array}{} \displaystyle R=0,\quad r=\frac{\Theta ^2}{2 \kappa}+\sqrt{\frac{\Theta ^4}{4 \kappa^2}-\frac{3\alpha(1-3c_\iota^2)}{16}}, \\ \displaystyle \nu=\nu_i, \quad \quad N=0. \end{array}$$(43)

Observe that the radius obtained is according to [11] when the triaxial body tends to be axial-symmetric.

Worth noting is that when the orbital and the rotational planes present a critical inclination, i.e. 1−3cι2=0$\begin{array}{} \displaystyle 1-3c_\iota^2=0 \end{array}$ , the formula for the radius is simplified to r = Θ²/κ. Notice that by analyzing (18) and (19) it can be concluded that there is no a simultaneous equilibrium for δ and ψ. As such, we continue studying the existence of 2-tori originated by fixing each one of the remaining angles individually.

Consequently the mean motions of θ and δ provide us periodic orbits withing 𝕋²(θ, δ).