Analyzing the spatial motion of a rigid body subjected to constant body-fixed torques and gyrostatic moment

This paper aims to explore the rotatory spatial motion of an asymmetric rigid body (RB) under constant body-fixed torques and a nonzero first component gyrostatic moment vector (GM). Euler's equations of motion are used to derive a set of dimensionless equations of motion, which are then proposed for the stability analysis of equilibrium points. Specifically, this study develops 3D phase space trajectories for three distinct scenarios; two of them are applied constant torques that are directed on the minor and major axes, while the third one is the action of applied constant torque on the body’s middle axis. Novel analytical and simulation results for both scenarios of constant torque applied along the minor and middle axes are provided in the context of separatrix surfaces, equilibrium manifolds, periodic or non-periodic solutions, and periodic solutions’ extreme. Concerning the scenario of a directed torque on the major axis, a numerical solution for the problem is presented in addition to a simulation of the graphed results for the angular velocities' trajectories in various regions. Moreover, the influence of GM is examined for each case and a full modeling for the body's stability has been present. The exceptional impact of these results is evident in the development and assessment of systems involving asymmetric RBs, such as satellites and spacecraft. It may serve as a motivating factor to explore different angles within the GM in similar cases, thereby influencing various industries, including engineering and astrophysics applications.

attitude motion is provided.The effectiveness of integrated controls for attitude, taking into account a distributed delay, is demonstrated through the use of numerical modeling.A previous study conducted by 11 examined a scenario in which the rotational axis of the RB was affected by the GM and another moment around the same axis.The authors were able to achieve analytical solutions for the RB's motion that closely aligned with the body's physical properties, thus establishing the uniqueness of these solutions.In 12 , the issue for the RB's movement is examined when it is subjected to a constant GM that is due to potential and gyroscopic forces.The authors successfully derived three new solutions for the EOM, which are governed by three linear invariant relationships for the angular velocity vector components.For the scenario where the RB is considered to be heavy with mass distribution, they obtained a solution that aligns with the Kovalevskaya and Goryachev-Chaplygin generalized conditions.
In 13 , the authors presented new precise solutions for the rotary movement of an RB analogous to Lagrange's conditions.These solutions pertained to cases where the RB is subjected to a constant external torque of magnitude.Specifically, the solutions are derived for the following scenarios; firstly when the torque is parallel to the axis of symmetry and for arbitrary initial angular velocity; secondly for an orthogonal torque on this axis with stationary rotation around that axis, besides the assumption of arbitrary initial angular velocity; and finally when both the torque and initial angular velocity are perpendicular to the axis of symmetry, with the torque being fixed to the body.The kinematic solutions are represented using the rotation matrix.The obtained exact solutions are applicable to any duration of motion and rotation amplitude.In 14 , the behavior of RBs undergoing perturbed rotations near regular precession according to Lagrange's case is investigated.The influence of a restoring moment and a slowly varying perturbing one are taken into account during the processes of the gained solutions.In the absence of resonance, an approximated system of EOM is derived for this nonlinear two-frequency system.
The required solutions are obtained in 15 for the overall rotary movement of a nearly symmetrical RB when the action of variable torques is considered.The authors specifically focused on the scenarios of acted constant torque along the rotation's axis, as well as variable transverse torques.In the case of RB with axial symmetry and consistent axial torque, the solutions of Euler's EOM are completely accurate.However, the solutions of the Euler's angles are given in approximate form.In order to consider the scenario of a rotating RB that experiences a time-varying torque in the axial direction 16 , it is necessary to expand upon the approach outlined in 15 .The resulting analytical solutions described the overall attitude movement of a near-symmetric RB that was subjected to time-varying torques around all three spatial axes relative to the body.In 17 , the analysis focused on studying the movement of an axisymmetric gyrostat satellite in a circular orbit within the impact of a Newtonian force field (NFF) is presented.It identifies and examines all the stable positions of the satellite within the orbital coordinate system, while also investigating the factors that determine their existence.Furthermore, the authors identified the specific values of the system parameters that trigger changes in the number of EPs.
In 18 , a study was conducted on the rolling of an asymmetrical RB on a horizontal plane when it is acted by a periodic GM.The authors approached the problem using a rubber body model, which assumed no slipping or spinning at the contact point.The results showed that under specific values of the system's parameters and time-dependent of the GM, the system exhibited acceleration which resulted in an unbounded growth of energy.Further investigations were carried out to analyze how the acceleration depends on the system's parameters and initial conditions.It should be noted that the small parameter approach has been widely utilized in 19 as one of the approximate methods to obtain analytical solutions for the RB problem.However, the obtained periodic solutions using this approach, whether in a uniform gravity field or in an NFF; contained singular points.These singular points presented a significant challenge because the solutions aren't defined as whole numbers or their negative counterparts.Consequently, it was crucial to address these singularities for all values of these frequencies.As a result, a significant amount of scientific research is required to bridge this gap, making it impossible to find a solution that is completely free of these singularities.In 20 , this problem was addressed by incorporating the effect of the third component of the GM, which led to the discovery of a new frequency, known as Amer's frequency.This achievement was confirmed when considering the complete impact of the GM, regardless of whether the motion of a symmetric or asymmetric RB.It was determined that the solutions obtained were free from any irregularities and were applicable for all values of this frequency.
In 21 , the author explored the analytic solution of free rotary movement of an RB that is powered by a lowpower motor.Through the application of asymptotic methods, it has been shown that the motion of the carrier body is closely related to the rotation around a stationary axis, which depends on the problem's parameters and initial conditions.The analysis in 22 focuses on investigating the equilibrium attitude and stability of a rigid spacecraft in a stationary orbit around a uniformly rotating asteroid.The linearized EOM governing attitude motion are derived based on the assumption of small motions.Subsequently, the equilibrium attitude is established for both a general and a symmetrical spacecraft.Owing to the presence of higher-order inertia integrals, the equilibrium attitude deviates slightly from zero Euler angles.In 23 , the study extends the inquiry into attitude stability to encompass a rigid spacecraft in a stationary orbit around a uniformly rotating asteroid.The authors observe that, owing to the markedly non-spherical shape and swift rotation of the asteroid, the attitude stability domain undergoes significant alterations compared to the classical stability domain predicted by the Beletskii-DeBra-Delp method for a circular orbit in a central gravity field.Notably, when the spacecraft is positioned along the intermediate-moment principal axis of the asteroid, the stability domain may exhibit a complete divergence from the classical stability expectations.In 24 , the investigation employs a differential geometric methodology to derive the Poisson tensor, Casimir functions, and equations of motion governing the phase flow and phase space structures inherent to the studied system.The equilibrium attitude of the spacecraft, serving as a stationary point for the equations of motion, is determined from a holistic perspective by considering the Hamiltonian constrained by Casimir functions.Subsequently, nonlinear stability conditions for the identified equilibrium attitude are formulated utilizing an adapted energy-Casimir method.The research further delves into the examination of nonlinear attitude stability concerning three key asteroid parameters, specifically the ratio of mean radius to orbital radius and the harmonic coefficients.In 25 , a method for depicting the turn-tensor of an axisymmetric RB using the angular momentum vector was suggested.The author proved that when specific external moments are applied, the movement of an axisymmetric RB is essentially the same as that of a spherical RB, except for the extra rotation around its axis of symmetry.Additionally, an accurate solution to the problem of unrestricted rotation of an axisymmetric RB was constructed when the impact of linear viscous friction was considered.Under the influence of the GM and NFF, an earlier investigation of the RB's problem with a zero-value assigned to the first component of the GM vector was found in 26 .To address scenarios with irrational frequencies, EOM was solved using the Poincaré method of small parameters.The influence of the GM, CBFTs, and resistive forces on a charged RB has been examined in 27 .A suitable governing system for EOM was approached using the averaging method.To reach the required results, Taylor's method was used along with some initial conditions to solve the averaged system of the EOM.
In this paper, the rotary spatial motion of an asymmetric RB with CBFTs, influenced by a first component of the GM, is investigated.To eliminate their reliance on both the inertia properties and the magnitude of torque, the controlling EOM in the case of CBFT is carried out in a dimensionless form.The determination of the dimensionless system's EPs, the derivation of the linearized EOM, the characteristic equation, and the stability features are also presented.The analytic solution for the scenario of applied CBFTs along the minor and middle axes is addressed by a comprehensive schematic simulation for the SS, trajectories, stability areas, and extreme values for the torques in the 3D phase plane.The numerical solution in the scenario of an applied CBFT along major axis is also presented, along with 3D and 2D histograms of the dimensionless angular velocities, resulting in a typical spin-up maneuver as expected in some of the analyzed regions.The impact of various GM values on the locomotion and stability trajectories is supplied because it provides a useful resource for outcomes in such a case.Finally, each scenario includes a detailed examination of the minimum and maximum values of different dimensionless angular velocity components.

Problem's formulation
In order to improve the description of this problem, this section seeks to provide us with more information.Therefore, the spatial rotation of an asymmetric RB about a fixed origin O of two coordinate systems is consid- ered.The first coordinate system Ox 1 y 1 z 1 is the inertial and the second system Ox 2 y 2 z 2 rotates with the body.The RB's motion is influenced by both GM and CBFTs M vectors, as graphed in Fig. 1.
The governing EOM for the body 1,6,27 , is given by Here, D = (D 1 , D 2 , D 3 ) is the tensor components of the principal inertia's moment such as D 1 > D 2 > D 3 , ω = (ω 1 , ω 2 , ω 3 ) represents the body's angular velocity, and (M 1 , M 2 , M 3 ) are the components of the CBFTs M , while 1 is the first component of the GM (where 2 = 3 = 0 ) along Ox 2 , and t represents the time.Taking (1) www.nature.com/scientificreports/into consideration the following new parameters to eliminate the dependence of the above system on the inertial properties of the RB as Therefore, one can rewrite the equation of system (1) as follows where the scaled quantities x j , µ j , and τ represent the angular velocities, CBFTs, and time, respectively.
The system in (1) can no longer rely on the magnitude of the torque vector k is introduced.Making good use of u to redefine x j , µ j , and τ as follows Substituting (4) into system (1) to yield where the adjusted quantities x j , µ j , and τ represent, respectively, the dimensionless forms of the angular veloci- ties, CBFTs, and time.According to system (4), one gets Therefore, the EOM at a steady state becomes where X j are the components of the EPs.Therefore, one obtains In other words, if and only if µ 1 µ 2 µ 3 ≥ 0 , the EPs exist.Eight EPs could be determined by knowing the value of the constant torque µ j as follows Let us define the state perturbation concepts of angular velocity as 3 Then taking into account the consistency of µ j , to derive the below linearized EOM Then, the characteristic equation for (11) can thus be expressed as ( 2) Vol.:(0123456789)The prerequisites for the existence and stability of the roots of Eq. ( 12) are given in Table 1 along with an extensive list of all achievable and possible combinations of CBFTs and angular velocities.

CBFT along the minor axis
In this section, an analytical approach that determines the angular velocities of the RB when undergoing CBFT along its minor axis is investigated.Therefore, at the value (µ 1 , µ 2 , µ 3 ) = (0, 0, 1) along this axis, system (5) takes the form These equations depict the EPs of the system, which form a characterized hyperbola according to the equation 1 , a new variable α 1 can be inserted as which also can be rewritten as dα 1 dτ = x 3 (τ ); α 1 (0) = 0 , and then one can transform system (13) into By taking the derivative of the second equation from (15) with respect to α 1 , and subsequently employing the first equation to get which is a simple harmonic motion equation.Therefore, the following solutions hold for the first two equations of system (15)   where Table 1.Explores the linear stability of different EPs.17) into the third equation in system (15), yields Let's assume that As a consequence, Eq. ( 18) undergoes reformulation as follows At the equilibrium's angle α * with the absence of the GM, sin α * = −2 F 2 is satisfied.Evidently, there is no EP for F < √ 2 .To investigate the system's stability properties in the scenario where F ≥ √ 2 will be presented.Therefore, the following new parameter is introduced in the form This means that Eq. ( 21) has been transformed in the (α, x) phase plane as Based on Eqs. ( 22) and ( 23), one can write Likewise noticing that for a certain trajectory, the motion's first integral F is constant.Therefore, Eq. ( 24) can be integrated to get where the second integral of the motion corresponds to the dimensionless energy constant G 1 .It may be contour- plotted on the (α, x) phase plane using Eq. ( 25) and provided a specific value for the dimensionless first integral of the motion F .This type of contour-plot is shown in Fig. 2 for the value F = 4.72 , where each dashed or solid line contour corresponds to a G 1 trajectory with constant values.The phase plane axis is periodically positioned along unlimited stable centers ( • ) and unstable saddles ( • ) when F ≥ 1.41.

Case EP Characteristic equation Stability Manner CBFT Notes
Therefore, (α * , x * ) = [sin −1 (−2 F 2 ), 0] , where sin −1 ϕ 1 ≤ π 2 for all ϕ 1 , is the first EP inside (α, x) con- cerning left of the origin.Center points (CPs) are situated at π − α * ± 2mπ ; (m = 1, 2, . ..), while the saddle points (SPs) are at α * ± 2mπ .For an SP α s , the energy constant G 1s takes the form as Moreover, the separatrix's equation travels via α s is Substituting ( 26) into (27) yields Bringing back that α(0) = 2ψ 1 and |ψ 1 | ≤ π .Regarding stability analysis, there are just two separatrices that are relevant: the left separatrix (LS), going via the left SP at α SL and wrapping around the left CP at α CL ; the right separatrix (RS), going via the right SP at α SR and wrapping around the right CP at α CR .The LS and RS at F = 4.72 are displayed in Fig. 3.The various CPs and SPs connected to the LS and RS are also represented.Additionally, it illustrates the cyclical, spin-up, and vertical crossing (VC) trajectories throughout the action.
The graphed curves in Fig. 4 expresses the impact of distinct GM's values 1 (= 100, 300, 500) kg m 2 s −1 that acted on the body's main axes of inertia.The distinction between every case has been shown in the graphical simulation as in (a) of Fig. 4 which represents the LS and RS at 1 = 100 kg m 2 s −1 , noticing that the cyclical trajectories at the RS make closed oval trajectories bigger than the LS with a united direction for the VC and the spin-up trajectories.The left SP α SL and right one α SR are located, respectively, at (0, 0) and (6.2, 0) , while the left CP α CL and the right one α CR are found, respectively, at (−3, 0) and (3.3, 0) on the (α, x) phase plane.At 1 = 300 kg m 2 s −1 , as seen in Fig. 4b, the cyclical trajectories at the RS make distinct closed ovals other than the aforementioned at 1 = 100 , which are now smaller than the LS, with a non-united direction for the VC and the spin-up trajectories.In this case, the left SP is located at (0.19, 0) and the position of the right SP is found at (6.15, 0) , while the left CP is found at (−3.1, 0) and the location of the right CP is found at (3.29, 0) .When the value of the GM becomes 1 = 500 kg m 2 s −1 , the difference between the two SS becomes so obvious as the LS is becoming much smaller than the RS with a distribution in the direction of the trajectories, see Fig. 4c.The left SP is now located at (0.32, 0) and the right one SP is found at (5.8, 0) , while the left CP and the right CP are located at (−3.13, 0) and (3.28, 0) , respectively.Considering The substitution from (29) into (27), yields It must be noted that cos α s = (F 4 − 4) A 2 and sin(α s 2) = ± (1 − cos α s ) 2 .Therefore, the substitution from ( 26) into (30), yields Parts of Fig. 5 show 3D representations of the LS and RS surfaces at D(= 2000, 1500, 1000) kg m 2 when 1 = 100 kg m 2 s −1 in the space (x 1 , x 2 , x 3 ) .It is noted that the mentioned space is divided into three infinite areas by two surfaces.Two of them are stable, while the third one is an unstable region.One of these areas will remain associated with any motion that has starting conditions there.If the body starts its motion within one of the stable domains, it will follow a closed and limited path within that region, revolving around its respective center.However, if the initial conditions lie within the realm of instability, an upward spin will be executed.When these conditions lie in the unstable region, the projection of the endpoint dimensionless angular velocity (x 1 , x 2 , x 3 ) inside (x 1 , x 2 ) plane represents a complete circle, and for initial conditions within any of the stable regions, it represents a segment of the circle.The circle's radius, denoted as F , has the flexibility to assume any length, leading to the possibility of both stable and unstable movement.
Figure 5 displays a contour map of the top area (x 1 , x 2 , x 3 > 0) for the LS and RS seen in Fig. 3.The EPs for the system are found on the two branches of the hyperbola x 2 (x 1 + b) = −1 .The equation x 1 + x 2 = 0 serves as a dividing line, creating distinct regions within the hyperbola branches known as Lyapunov stable and unstable (31) www.nature.com/scientificreports/domains, as shown in the cylindrical coordinates F = x 2 1 (0) + x 2 2 (0) and ψ 1 = tan −1 [x 1 (0) x 2 (0)] , which are presented.Within the circumference of a circle with a radius of F , the absence of EPs can be observed.
Figure 6 shows the depiction, at a moment of principal inertia D(= 2000, 1500, 1000)kg.m 2 , the initial condi- tion x(0) = (2.5, 4, 0.1) , and the value of the GM 1 = 100kg.m 2 .s−1 , for the following surfaces: 1. the circle (x 2 1 + x 2 2 ) 2 − 4 = 0 where none of the EPs are located, 2. the two branches of the hyperbola x 2 (x 1 + b) = −1 , which are where all of the system's EPs are located, 3. the line x 1 + x 2 = 0 that separates the two branches of the hyperbola x 2 (x 1 + b) = −1 into regions with Lyapunov stable and unstable regions, 4. the surface The SP, CP, and VC points connected to the LS and RS in Fig. 2 have the same positions where a circle of radius F = 1.41 intersects with the SP, CP, and VC curves in Fig. 6.
In order to find the solutions to the transcendental equations provided by Eq. (28), one must find the intersection of the two SS and the (x 1 , x 2 ) plane as follows The fact that α = α s is thought to be a solution of (29), which results in The SP, CP, and VC connected to the LS surface converge into a single focal point at (x 1 , x 2 , x 3 ) = (1, −1, 0) , whereas those connected to the RS surface converge into (x 1 , x 2 , x 3 ) = (−1, 1, 0) .The statement holds true in the (32) case where F = 1.41 .So, starting with F = 1.41 in Fig. 3, one may create the SP, CP, and VC curves by tracing the polar coordinates (F, ψ 1 = α 2) covered by the SP, CP, and VC points as F is raised.The motion that started from the separatrix area will stay trapped there and approach steadily upon the EP in the location of the intersection of (x 2 1 + x 2 2 ) 2 − 4 = 0 with the unstable branch of the hyperbola x 2 (x 1 + b) = −1 , which is encapsulated in that area.The projection of the trajectories on the (x 1 , x 2 ) plane forms a segment of a circle with a radius of F = 1.41 since it is in the stable zone bounded by the RS surface, which is plotted in Fig. 7.
By modifying system (13), one may determine the extreme values that the dimensionless angular velocity component x 1 , x 2 , and x 3 takes along an enclosed trajectory in each of the two regions of stability with constant torque along the minor axis as follows  It's crucial to mention that the extremal values for x 1 and x 2 will only take place when the dimensionless angular velocity path intersects the plane (x 1 , x 2 ) or the plane (x 1 , x 3 ) .In contrast, the extreme values for x 3 will only take place when this trajectory intersects with the two surfaces x 2 (x 1 + b) = −1 , which is a direct result of the system (34).By replacing the formulas for the constants of the motion F and G 1 , provided by Eqs. ( 18) and (30), with the system (34), and calculating the resultant transcendental equations, it is possible to determine the extremes for x 1 , x 2 , and x 3 .
(34)  Considering the following initial data to continue our estimation as the initial values of the scaled angular velocity components x(0) = (2.5, 4, 0.1) , the values of the principal inertia D(= 2000, 1500, 1000)kg m 2 , and the value of the first component of the GM 1 (= 100, 300, 500)kg m 2 s −1 .Then, the values for both constant of the motion are F = 4.72 and G 1 = 2.54.
As a result, Eq. ( 18) implies The substitution of (30) into (35) produces The conditions associated with the x 3 extremes are provided by Equation ( 38) is substituted into Eq.( 18) to get Then substituting the previous solutions for x 1 and x 2 into Eq. ( 30 to produce the required extreme values of x 3 .
Table 2 provides an analytical statement for the extreme values of several dimensionless angular velocity components x 1 , x 2 , and x 3 alongside a solution in CBFT along minor axis which is periodic.

CBFT along the major axis
In this section, an approximate solution for the angular velocities of the RB when it is subjected to CBFT along its major axis is explored.This method is employed when the analytical solution fails to be achieved and will be further discussed later on.If the constant positive torque along the major axis has the form (µ 1 , µ 2 , µ 3 ) = (1, 0, 0) , then system (5) can be rewritten as follows The equations presented in the above system illustrate their EPs, wherein they manifest as a hyperbola with the equation X 2 X 3 = −1 .Additionally, they are stable at X 3 > 1 and unstable at |X 2 | > 1 .Adding a new variable α 2 through the following equation as (38) Table 2. Presents x j extreme values corresponding to the minor axis as (µ 1 , µ 2 , µ 3 ) = (0, 0, 1).
Vol:.( 1234567890 www.nature.com/scientificreports/which also has the form dα 2 dτ = x 1 (τ ); α 2 (0) = 0 .Applying this variable will transform the system (41) into It is significant to note that the variables of system (43) can't be separated, and then it is impossible to obtain the analytical solution of this system by using the aforementioned method.Therefore, a numerical solution for the system is going to be presented and diagrammed to see the various effects of the GM on the RB motion.All numerical simulations were performed concerning the aforementioned values of the inertia's main moments D(= 2000, 1500, 1000)kg.m 2 besides the values the GM 1 (= 100, 300, 500) kg.m 2 .s−1 .
Currently, the rotational motion of the RB, with the assumption that the torque along the major axis remains constant, is going to be examined.In this particular scenario, Figs. 8, 9, 10, 11, 12, 13, 14, 15 are drawn to depict the equilibrium conditions on the specific regions mentioned in Table 1.
Figure 8 describes the 3D trajectory of the angular velocities x 1 , x 2 , and x 3 of the RB when the GM has the values 1 (= 100, 300, 500) kg m 2 s −1 .On the other hand, the projections in the planes x 1 x 3 , x 2 x 3 , and x 1 x 2 are presented in Figs. 9, 10, and 11, respectively, at the same values as Fig. 8.These figures are generated with the initial condition x(0) = (4.5, 1.8, 3.4) at the upper unstable region of the space (x 1 x 2 x 3 ) during the interval τ ∈ [0, 15] s.
Note that in all 3D graphical simulations, the motion's trajectory begins with rainbow colors and then transitions to a unique color during its path interval.
In Fig. 8, as the GM value equals 100 kg m 2 s −1 , the body begins its spin in the positive area of x 1 and remains spinning between the intervals x 2 ∈ [−4, 4] and x 3 ∈ [−4.1, 4] for its oscillation, where it eventually converges on a spin-up maneuver.As the GM value approaches to the value 300 kg m 2 s −1 , it's released that the body follows the same path as before, but it is noticed that an increase in the amplitude of the spin maneuver in the positive region of x 1 towards the negative regions of x 1 , x 2 , and x 3 , as seen in the 2D simulations of Fig. 9, 10, 11.The increasing of the frequency of the body oscillation increases the oscillation interval as x 2 ∈ [−4.2, 4] and x 3 ∈ [ −5, 4] .At 1 = 500kg.m 2 .s−1 , the body spins increasingly to maneuver around x 1 .The amplitude of the maneuver about x 2 and x 3 is also increases as x 2 ∈ [−4.5, 4] and x 3 ∈ [−6, 4] , and then the closed paths of the trajectories keep it spinning manner with a slight deformation for the one path trajectory before the increasing of the GM value, as shown in Fig. 10.The importance of these outcomes lies in their broad utilization across various domains.For instance, they can be employed to stabilize the movement of spacecraft in their orbits.This is achieved by adjusting the GM value, which impacts the component of the body's angular velocity and ensures it remains within the designated orbit. (42) Figure 8. Presents a 3D plot for x 1 , x 2 , and x 3 at a constant major axis torque when 1 (= 100, 300, 500) kg m 2 s −1 .
Figures 12, 13, 14, 15 are calculated at the stable region of the space (x 1 x 2 x 3 ) during the interval τ ∈ [0, 15]s when the initial condition x(0) = (0, −1, 5) is considered.The outcomes of this simulation are graphed at distinct values of the first component of the GM 1 .At 1 = 100 kg m 2 s −1 , the trajectory initiates a customary rotational acceleration maneuver inside positive region of angular velocity x 3 with a slight decrease in its amplitude towards the negative region and the body has a stable motion, as explored in Fig. 12a.As the GM increases to become 300 kg m 2 s −1 , the body starts its spinning motion about the positive stable region of x 3 , and then moves towards the SS with a higher oscillation amplitude and lower frequency.Therefore, it crosses into the positive unstable region of x 1 , where it converges to a pure spin-up maneuver, as shown in Fig. 12b.At 1 = 500 kg m 2 s −1 , an increase in the amplitude of the spinning about the negative region of x 1 is noticed, and the trajectory follows the same path as 1 = 300 kg m 2 s −1 with the same motion's behavior, as observed in Fig. 12c.This description can be easily noticed in the 2D representation of the body's angular velocities, as seen in Figs. 13, 14, 15.

CBFT along the middle axis
In this section, the analytical method employed to calculate the angular velocities of the RB while it experiences CBFT along its middle axis, is explored.In the scenario, we consider (µ 1 , µ 2 , µ 3 ) = (0, 1, 0) for this torque along the middle axis, and then equations of system (5) provide the following form  In this scenario, the points of equilibrium in the system form a hyperbola resembling the form (X 1 + a)X 3 = 1 .As mentioned in the preceding three sections, we consider a new variable α 3 such that (44)   Figure 15.Plots 2D graphs for the drawn angular velocities x 1 and x 2 in Fig. 12.
The enclosed pathway's 3D simulation, maintaining a constant energy G 2 is given by Eq. (71).It is calculated for the scenario that involves CBFT along central axis in addition to a starting condition |x 1 (0)| � = |x 3 (0)|.
Figure 16 provides a visual representation of the 3D simulation depicting SS of Eq. ( 71).These SSs divide the space (x 1 x 2 x 3 ) into three infinite regions: two stable regions and an unstable one.Any motion originating from one of these regions will continue to be associated with that region.If the initial conditions for motion lie within either of the stable regions, the trajectory will be closed and restricted, revolving around its corresponding center.On the other hand, if the initial conditions fall within the unstable region, a spin-up maneuver will occur.
Figure 17 shows the drawing at initial condition x(0) = (1.7,1.1, −0.2) for the following surfaces: (71)  www.nature.com/scientificreports/The system's dimensionless energy constant is denoted by G 3 , which corresponds to the system's second integral of motion.
As performed in the case of x 1 (0) = x 3 (0) , stability analysis shows that every combination of initial condi- tions concerning with x 1 (0) = x 3 (0) > 0 leads to a confined oval-shaped closed periodic solution within the plane x 1 (0) = x 3 (0) and encircles the EP ( −a+ √ a 2 +4 2 , 0, −a+ √ a 2 +4 2 ) contained within this plane, while every set of initial conditions that satisfy x 1 (0) = x 3 (0) < 0 leads to a confined oval-shaped closed periodic solution within the plane x 1 (0) = x 3 (0) and encircles the EP ) contained within this plane.Therefore, it can be inferred that, similar to the scenario of x 1 (0) = x 3 (0) , the two extremes of x 1 (0) corresponding to x 3 (0) are positioned at the point where the angular velocity trajectory intersects the plane x 2 (0) = 0 .Similarly, the two extremes for x 2 (0) are situated at the intersection points of the trajectory with the surface (x 1 + a)x 3 = 1.
Table 4 summarises the extremal values reached by x 1 = x 3 and x 2 along a closed trajectory with an initial condition at x(0) = (2.5, 1.3, 2.5).

Unstable separatrix
For the case when x 1 (0) = −x 3 (0) , equations in (47) yield Hence, equations in (50) produce Therefore, x 1 = b + x 3 , L = −M, and then N = 0 .Equation (51) results It is evident that in this scenario, there are no EPs.A closer look at Eq. (81) shows that Utilizing Eqs.(80) and (81) in order to achieve Therefore, a movement that quickly transitions into a pure spin in an upward direction around the central axis, while maintaining a constant angular acceleration is obtained.

Discussion
The investigation delves into the influence of the GM and CBFTs on the rotatory motion of an asymmetric RB, employing Euler's dynamic equations to derive the governing EOM.To reduce reliance on inertia characteristics, the equations are appropriately scaled.By identifying and expanding the controlling EOM, inertia properties are effectively removed from the equation.The EPs of the dimensionless system are determined, alongside the linearized EOM, characteristic equation, and stability properties.
2. When considering a directed CBFT along the major axis, numerical solutions are provided, along with 3D and 2D graphs depicting dimensionless angular velocities leading to a typical spin-up maneuver.Motion stability varies with increasing GM values, while maintaining the spin-up maneuver about one of the dimensionless axes.3.In the scenario of CBFT along the middle axis, an analytic solution is obtained, accompanied by complete 3D and 2D phase plane representations for various SS motions.Extreme value cases are explored, and stabilization is discussed in detail, including system solutions and tables of extreme values in the stable separatrix section.
The effects of different GM values on body paths and stabilization are thoroughly examined, yielding beneficial insights.Each instance undergoes a comprehensive assessment, analyzing obtained values at the lowest and highest points of angular velocity distinct dimensionless components along a periodic solution.

Conclusion
The influence of the GM and CBFTs on the rotatory motion of an asymmetric RB is investigated.The governing EOM has been derived using Euler's dynamic equations and scaled to reduce their dependence on inertia characteristics.To remove their reliance on inertia properties, the controlling EOM has been identified and expanded.The EPs of the dimensionless system are determined, as well as the linearized EOM, characteristic equation, and stability properties.In three distinct cases.
The remarkable significance of these results becomes apparent when considering the advancement and evaluation of systems that utilize asymmetric RBs, such as satellites and spacecraft.The study's findings hold particular promise for stabilizing spacecraft motion within their orbits.Manipulating external moment values and body parameters can further enhance stabilization efforts.Consequently, these results may drive further exploration of GM perspectives in similar scenarios, with potential impacts across diverse industries, including engineering and astrophysics applications.

Figure 1 .
Figure 1.The simulation model of the RB's motion.

Figure 6 .
Figure 6.Shows contour plots of LS and RS surfaces for a constant minor torque axis.

Figure 7 .
Figure 7. plots of typical steady trajectories for a constant minor torque axis.

Figure 9 .
Figure 9. Shows 2D plots for the represented angular velocities x 1 and x 3 in Fig. 8.

Figure 10 .
Figure 10.Shows 2D plots for the given angular velocities x 2 and x 3 in Fig. 8.

Figure 11 .
Figure 11.Shows 2D plots for the graphed angular velocities x 2 and x 3 in Fig. 8.

Figure 14 .
Figure 14.Plots 2D graphs for the drawn angular velocities x 2 and x 3 in Fig. 12.

Figure 16 .
Figure 16.Presents a 3D simulation for SS for a constant middle torque axis.

Figure 17 .
Figure 17.Shows projection of SS in the plane (x 1 x 3 ) for a constant middle torque axis.