Lifetime and Dynamics of Natural Orbits around Titan

: Considering the growing interest in sending probes to the natural satellite Titan, our work aims to investigate and map natural orbits around this moon. For that, we use mathematical models with forces that have symmetry/asymmetry phenomena, depending on the force, applied to orbits around Titan. We evaluated the effects due to the gravitational attraction of the Saturn, together with the perturbative effects coming from the non-sphericity of Titan (the gravitational coefﬁcient J 2 ) and the effects of the atmospheric drag present in the natural satellite. Lifetime maps were generated for different initial conﬁgurations of the orbit of the probe, which were analyzed in different scenarios of orbital perturbations. The results showed the existence of orbits surviving at least 20 years and conditions with shorter times, but sufﬁcient to carry out possible missions, including the important polar orbits. Furthermore, the investigation of the oscillation rate of the altitude of the probe, called coefﬁcient ∆ , proposed in this work, showed orbital conditions that result in more minor oscillations in the altitude of the spacecraft.


Introduction
Being one of the largest natural satellites in the solar system, Titan is an exciting body due to its similarities with Earth in many ways, making it one of the targets to be visited by space exploration missions planned for the coming years. After its discovery in the mid-17th century by Christiaan Huygens, the main findings of Titan came only centuries later during astronomical observations led by Gerard Kuiper, published in 1944. These detected methane by using sunlight reflected from Titan that passed through a spectrometer [1], and therefore it was concluded that this moon had an atmosphere. Scientists believed that this atmosphere would be dense and opaque, thanks to the advances in studies and observations made later. With the beginning of the space age, and the overflights on Titan by the Pioneer 11 and Voyager 1 and 2 missions, several predictions were made regarding its thick and extensive atmosphere, as well as its composition, in addition to other physical characteristics such as the temperature of its surface, its radius, and others [2]. However, it was with the Cassini-Huygens mission that more detailed measurements of the surface of the moon could be carried out because, due to its dense atmosphere, the measuring devices of the Voyager 1 and 2 missions were unable to obtain clear images, as was done when the Huygens probe landed on its surface. The data brought by the Huygens probe showed that the atmosphere of Titan is active [3] and that important chemical-organic processes take place in it. It was also seen that climatic cycles based on hydrocarbons occur, which are very similar to the current water cycle on Earth [4], that is, there is the formation of clouds and precipitation in addition to the existence of seas and lakes in the liquid state. Therefore gravitational perturbations coming from Jupiter, the third body, and the non-uniformity in the distribution of the mass of the satellite. The gravitational potential of the body was considered with the terms due to the second and third Zonal harmonic (J 2 and J 3 ) and to the order 2 sectoral flattening coefficient (C 22 ). It used Hamiltonian formalism, considering a Hamiltonian with explicit time dependence, as well as the method of averaging over elliptical orbits for the space probe. The results compared the useful life of the studied orbits, making calculations using one and two averages over the motion. Ref. [21] proposed to investigate the stability time of a space probe around a natural satellite in a three-body system, using the Callisto-Jupiter system as a reference but assigning different values for the eccentricity of the disturbing body, considering that the natural satellite is in the center of the reference system and Jupiter the third body. In [22], an analytical and numerical study was carried out considering the lifetime of a probe around Europa in the presence of Jupiter as a third body. The orbits of interest were orbits with low altitude but highly inclined, since this type of orbit takes advantage of the rotation of the celestial body to observe its entire surface more quickly. The simulations took into account the non-coplanarity of Jupiter and Europa. Their results showed lifetime maps for different conditions for the pericenter of probe and initial node and favorable conditions to extend the lifetime of these orbits using propulsive systems. With the same aim of studying moons of Jupiter, ref. [23] carried out studies for Enceladus in search of favorable orbit conditions for carrying out the proposed missions for the natural satellite.
Considering the growing interest in carrying out missions to Titan, our work aims to contribute to future missions through the study and mapping of orbits around this natural satellite.

Mathematical Model
Since the objective of this work is to study the stability and lifetime of the natural orbit of a probe around the natural satellite Titan, simulations were carried out assuming the presence of just three bodies: the natural satellite Titan, as the central body; a probe with negligible mass, which orbits Titan in orbit with initial three-dimensional conditions, with orbital elements (a 0 , e 0 , I 0 , ω 0 , Ω 0 , M 0 ); and the planet Saturn, acting as a third disturbing body, describing an orbit, with respect to Titan, that is circular and coplanar. Information about the masses and equatorial radii of Titan and Saturn are presented in Table 1. Table 1. Mass values, mass ratios, and the equatorial radius used in the simulations for the planet Saturn and the natural satellite Titan. µ i is the normalized mass of the bodies.

Body
Mass (kg) µ i R (km) Saturn 5.68319 × 10 26 0.99976335 5.8232 × 10 4 Titan 1.34520 × 10 23 0.00023664 2.5747 × 10 3 Considering the interest in the study of natural satellites and space probes located in regions extremely close to these satellites, the effects due to symmetry or asymmetry of the mass configurations of these natural satellites and their gravitational fields are significant. Thus, carrying out studies considering the non-uniform distribution of Titan, more precisely, we consider the second Zonal harmonic J 2 , which corresponds to the flattening of the natural satellite. Table 2 presents the non-sphericity coefficient J 2 and the orbital conditions of Titan. Previous works [13,18] showed the importance of the disruptive effects on the orbit of the probe due to the irregular shape of the natural satellites. The coefficient due to the flattening of Titan is more significant in terms of magnitude. As is known from previous work [24,25], a considerably dense atmospheric region is present in Titan. Hence, the effects due to drag are essential for an accurate determination of the dynamics of probes and space vehicles around Titan. Therefore, in our work, we also consider these effects. In general, the equations of motion for the system under study can be described according to [11,26].
being thatr andr where, in Equation (1), r and r S are the position vectors of the probe and Saturn, relative to Titan, and µ p , µ and µ S the mass ratios of the probe, Titan and Saturn, respectively. In Equation (2) to Equation (4), which correspond to the components of the acceleration due to the flattening of Titan, R is the equatorial radius of the natural satellite and J 2 is the second harmonic zonal term. In Equation (5) to Equation (7), the components due to the atmospheric drag force are present, where ρ, C d , and S are the atmospheric density, the drag coefficient, and the probe section area, respectively. The Saturn and Titan conditions were implemented in a code written in the C language, using the REBOUND package [27], which integrated the restricted three-body problem (RTBP). The simulations were performed with the IAS15 algorithm included in the REBOUND package. The IAS15 is a 15th-order integrator that can handle conservative and non-conservative forces with adaptive step-size control. It is based on Gauss-Radau quadrature [28] and is accurate down to machine precision [29]. We defined the initial orbital conditions of the probe, (a 0 , e 0 , I 0 , ω 0 , Ω 0 , M 0 ). During the temporal evolution of the orbit, the occurrence of escapes from the region of interest of the study were considered, and the simulations were carried out also assuming a collision criterion between the probe and Titan. Considering that this natural satellite has a significant atmosphere, which cannot be neglected in the dynamics, we assumed two criteria for the loss of the probe in our simulations. These collision and escape criteria are described in Sections 2.1 and 2.2.
All simulations were made using a time interval corresponding to 477 orbital periods of Titan, which corresponds approximately to 20 years. Orbits that did not escape the region of interest and were not classified as lost were considered as orbits surviving the total integration time.

The Region of Interest Ejection Criteria
In our work, we consider the region of interest for the orbits of the probes as the region of the gravitational domain of Titan, with its gravitational radius of influence being calculated based on [26], by Equation (8): where µ s is the normalized mass ratio of Titan. In our work, we established as a criterion for escape, or ejection of the probe from the region of interest, an instant when the orbit of the probe has an altitude, in relation to surface of Titan, greater than the radius of the sphere of influence of the natural satellite. Thus, when observed that the altitude of the probe, at time t, was greater than the radius of the sphere of influence of Titan, the probe was considered ejected, being removed from the simulation.

The Collision Criteria
In addition to the escape criterion, we considered two other forms of probe loss in our work. First, taking into account the presence of an atmospheric region around Titan that can be divided into a denser region, Alt ≤ 600.0 km, and a less dense region, 600.0 km < Alt ≤ 1300.0 km [24], we established a first collision criterion considering the loss of the probe upon, in the simulations, the first contact of the probe with the densest atmospheric region of Titan. Therefore, if the orbit of the probe altitude were equal to or less than 600 km, the probe would be removed from the simulation. This criterion was classified as (CollAtm) and can be illustrated by the dotted line in Figure 1. This criterion aims to analyze the lifetime of the orbit of the probe on the surface of Titan, motivated by the interest in using orbits to observe the atmosphere and the surface of the natural satellite in regions not disturbed by the action of more intense atmospheric drag. To later investigate the effects on the lifetimes of the orbits, a physical collision criterion (CollAlt) was considered in the simulations, thus taking into account the probe dynamics in the interval for Alt(t) ≤ 600 km. This criterion can be illustrated by the solid line on the surface of Titan in the scheme of Figure 1. Each of the criteria were chosen and simulated individually.

Results and Discussion
The results obtained in our investigations were gathered to present, initially, the simulations considering only the third-body perturbations due to Saturn. This is followed by the analysis considering, in addition to Saturn, the contributions due to non-sphericity of the Titan. Subsequently, the orbital regions that contribute to atmospheric drag present on Titan are also analyzed. Finally, we analyze the evolution of the orbital elements of some orbits and present some maneuver proposals.

The Third-Body Effects of Saturn on the Lifetime of Orbits
For a first analysis of the perturbative effects due to a third body, promoted by Saturn, on the orbit of the probe around Titan, simulations were initially performed considering the collision criterion in the region of the densest atmosphere (CollAlt), described in Section 2.2. In these simulations, the orbits of probe were initially considered circular, and the analysis was later expanded to eccentric initial orbits. For the case of initial circular orbits, investigations were initially carried out for an initial inclination interval between 0°≤ I 0 ≤ 180°and an initial semi-major axis interval between 1.25 R Titan ≤ I 0 ≤ 5.25 R Titan , which corresponds to initial altitudes in the range of 643.7 km ≤ Alt 0 ≤ 10,942.5 km. Then, simulations were performed with initial eccentric orbits, in the range of [0.1-0.5], considering the same initial conditions for inclination and semi-major axis, initially assuming that, for the cases of elliptical orbits, the orbital parameters pericenter argument, ω 0 , ascending node longitude, Ω 0 , and true anomaly, M 0 , were zero. The results obtained in the simulations were initially grouped into two groups: Lifetime overview and Lifetime in non-survival regions.

Lifetime Overview
For the first studies, assuming the loss criterion of the probe when entering the densest atmosphere of the Titan (CollAtm), based on Section 2.2, the simulations showed that there is a large region (see in Figure 2) over which the orbits of the probes survived the full integration time so, obtaining lifetimes of at least 20 years. In Figure 2, the limits of these regions are shown by isolines, using different colors for each initially simulated eccentricity in the range of [0.0-0.5]. These regions with orbits lasting at least 20 years are found, in the interval of inclination of 0°≤ I 0 ≤ 90°, in the lower regions of the isolines and the interval of inclination of 90°≤ I 0 ≤ 180°, in the superior regions of the isolines, corresponding to each value of initial eccentricity of the orbits. In the case of circular initial orbits (line in black) conditions extend throughout the initial semi-major axis range between 1.25 R Titan ≤ I 0 ≤ 5.25 R Titan , which corresponds to initial altitudes in the range of 643.7 km ≤ Alt 0 ≤ 10, 942.5 km. There is also a symmetry with the axis at 90°for the behavior of the surviving orbital regions. This behavior was observed for all eccentricities simulated. It is worth mentioning that, since in the simulations it is assumed that all orbits are implemented with zero true anomaly, the regions where the isolines start, for each initial simulated eccentricity, correspond to the initial orbits whose initial altitudes have a value greater than 600 km, because the regions before the beginning of the isolines on the lower inclinations, 0°, as shown on the maps, produce initials radius at the altitude of the densest atmosphere of Titan determined for this work. Based on the simulation criteria presented in Sections 2.1 and 2.2, the final conditions of the orbits were classified into three groups: orbits that collided, escaped the region of interest, or survived the total integration time. Overall, the results showed that orbits that did not survive the total integration time of 20 years were considered lost by collision. Therefore, no escapes from the region of interest were recorded. This behavior shows that these orbits entered the densest atmosphere of the Titan, Alt(t) < 600.0 km, at some point during their temporal evolution. This behavior was observed for both initial circular and eccentric orbits (0.1 ≤ e 0 ≤ 0.5).

Lifetime in Non-Survival Regions
Based on the behaviors observed in the maps of Figure 2, it was seen that for the entire initial eccentricity interval of 0.0 ≤ e 0 ≤ 0.5, the region for the initial inclination between 50°a nd 130°is a typical interval where the orbits of probes did not survive the full integration time, being considered lost for entering the atmosphere of Titan at some point in their temporal evolution. As it is interesting to investigate this region to see how the lifetimes of the orbits of probes behaved, we performed simulations in the region of 50°to 130°for the initial inclination of the interval and semi-major axis, 1.2R Titan ≤ a 0 ≤ 5.0R Titan , where R Titan corresponds to the largest equatorial radius of Titan, and eccentricity 0.0 ≤ e 0 ≤ 0.5. The initial inclinations of 65°to 120°are shown in Figure 3, where lifetime maps were built. In these maps, initial conditions whose semi-major axis already produces orbits with an initial orbital radius less than 600.0 km were omitted. The results obtained for the intervals of 50°to 65°and 120°to 130°, however, were only described in the text. The results show that, for the orbits around Titan that did not survive the total simulation time, the lifetimes show a dependence on the choice of the inclination value, the semi-major axis (orbital radius), and the initial eccentricity. There is a greater dependence on the eccentricity and the initial semi-major axis. The maps show that when the orbits are implemented with larger eccentricities, the lifetimes are smaller, and, for fixed initial eccentricities, the initial semi-major axis increases. So the increase in the initial orbital radius also contributes to shorter orbit lifetimes. These behaviors can be explained by the fact that the perturbation of the third body, Saturn, affects the eccentricity of these orbits, making them more eccentric in a shorter time, thus favoring the loss by collision in the atmosphere, or the escape of the region of interest of the probe. The influence of the choice of the initial semi-major axis can be better seen in the maps fixing the initial eccentricity because orbits with a higher initial semi-major axis have a larger initial orbital radius, thus exposing the probe to a more significant influence of Saturn, which tends to promote greater perturbations, leading to shorter lifetimes for the orbits. The choice of the initial inclination is also an important factor, although with less influence compared to the choice of a 0 and e 0 , because in the range of 65.0°≤ I 0 ≤ 120°the orbits decrease lifetime as the inclination increases from 65.0°to initial polar orbits and as they decrease from I 0 = 120°to I 0 = 90°, reinforcing the observations about the symmetry with an axis at I 0 = 90°. These results will be further detailed in Section 3.5.
In terms of the lifetimes of the probe orbits, in the range, 65.0°≤ I 0 ≤ 120°, Figure 2, shows that the initial circular orbits had more records of orbits with the longest lifetimes because orbits with durations in the range from 87 days to 3.18 years were found. For e 0 = 0.1, the lifetimes showed a sharp drop concerning the circular orbits, recording lifetimes in the range of 49 to 177 days. Orbits with durations between 17 and 104 days, for e 0 = 0.2, and between 30 and 74 days, for e 0 = 0.3, was also found. For more eccentric initial orbits, such as e 0 = 0.4, most orbits recorded time between 18 and 42 days in duration. When the simulations considered an initial eccentricity of e 0 = 0.5, the lifetimes were between 2.8 and 30 days in duration.
In the range of 50.0°≤ I 0 ≤ 60°and 120.0°≤ I 0 ≤ 130°, it is also observed that the initial circular orbits had more orbits with the longest lifetimes, because orbits with a duration in the range from 183 days to 4.3 years were found. For e 0 = 0.1, the lifetimes showed lifetimes in the range of 76 days to 1 years. Highlighting the orbit with a 0 = 5.25 R Titan and I 0 = 50°, which has a lifetime of approximately 7 years. For e 0 = 0.2 and e 0 = 0.3, orbits with lifetimes between 2 days to 1.7 years and 35 days to 1 year were registered, respectively. The lifetimes of the orbits for e 0 = 0.4 e e 0 = 0.5 have values between 18 to 55 days and 2 to 33 days, respectively.

The Effect Due to the Non-Sphericity of Titan
Considering the interest in studying orbital regions extremely close to natural satellites, the effects of the mass configurations of these natural satellites are essential. Some bibliographic references, such as [21,30,31], indicate that the flattening of the central body can offset the effects of a disturbing third body, helping to find more stable orbits. In this sense, in our work, we consider the effects coming from the non-sphericity of Titan in our simulations. More specifically, the simulations were performed considering the zonal harmonic coefficient of Titan J 2 (see Table 2). These simulations were performed, initially, considering the same conditions presented in Section 3.1.2, keeping the criterion of collision or loss of the probe by contact with the denser atmosphere of Titan. The results considering the zonal harmonic J 2 allowed us to observe that, in general, the behaviors of the lifetime maps do not change significantly in relation to the maps considering only the perturbations originated from Saturn. Based on the results obtained in the simulations considering only perturbations from Saturn and considering the perturbative actions from Saturn and the contribution due to J 2 of Titan, it was calculated the variations in these orbit lifetimes, which showed that most orbits maintained, in terms of months, the same chosen lifetimes without the action of then non-sphericity of Titan.
In terms of the initial eccentricity of the orbits, it appears that, in addition to the more significant number of long-lived orbits found in initial circular orbits, the effects due to the flattening of the natural satellite also acted more significantly in these orbits, registering slightly more significant gains or losses in the duration of the lifetimes of the orbits. However, for probes placed in an eccentric initial orbital, we observed some conditions with significant gains or losses in life were also encountered. That when the initial eccentricity increased, the variations in lifetime decreased. In terms of lifetime values, orbits implemented with an initial eccentricity e 0 = 0.5 did not gain more than 1 day in duration, and the duration of the orbit was in the range between 1 and 6 days. In the dynamics considering e 0 = 0.4, most orbits did not show gains or losses greater than 5 days. However, an orbital condition was found for the inclination of 50°and initial semi-major axis a 0 = 11126.38 km with an initial altitude of Alt 0 = 6675.8 km, where a gain of approximately 5.0 years in the life was observed. Due to flattening of Titan, the orbit had a total duration of 9.0 years with this gain. For orbits implemented with e 0 = 0.3, orbit times gains of more than 7 days were not found, as well as orbits with lifetime losses were also not larger than 5 days. Although in the simulations considering e 0 = 0.2, most orbits have registered gains or losses of less than 6 days, two orbits were found, with I 0 = 50.0°f or a 0 = 7981.6 km and I 0 = 130.0°for a 0 = 10427.5 km, with gains of 7 and 19 years. For e 0 = 0.1, most of the gains or losses were not of more than 6 days. In the case of the simulations carried out considering the initial circular orbits, the results showed that most of the orbits do not present changes in lifetimes larger than 1 day. However, for the region of lower initial altitudes, close to 600 km, the largest gains, and losses in orbit lifetimes were found. Due to the flattening of Titan, a large concentration of these changes is seen in the initial inclinations close to 50, 60, 90, 120, and 130 degrees. The orbits generated by these initial conditions had losses of up to 36 days. The observed increments were found, in general, in three large intervals: orbits with gains of up to 36 days, orbits with gains of up to 2 months, and orbits with gains of 3.5 months. Orbits were also found that, due to gravitational perturbations of Titan, became orbits surviving the total time of integration, as they showed gains of 15 and 19 years. These orbits have a 0 = 6797 km for I 0 = 50°and a 0 = 8445 km for I 0 = 130°, respectively. In general, the first results considering the zonal harmonic J 2 showed that, depending on the initial geometry of the orbit implementation possible, there are local conditions where the gains or even losses in the lifetime of the orbits can be considered for possible missions around Titan. Previous works such as [17,32] also show results showing results with similar behavior for probe studies around other natural satellites considering their effects due to J 2 . Furthermore, it was observed that the increments or losses in lifetimes in these simulations determined the duration of the orbits at altitudes above the denser atmosphere of Titan, as defined previously.

Contribution of the Altitude Region of Less Than 600 km in Lifetimes of the Probe
So far, the simulations carried out in Sections 3.1 and 3.2 had a criterion of loss of the probe as the first contact with the denser atmosphere of Titan (CollAtm). Aiming to investigate the contribution of this region with less than 600.0 km in altitude in the lifetimes of the orbits, new simulations were performed, replacing the probe loss criterion at an altitude of 600.0 km (ColLAtm) with a physical collision when the altitude of probe is zero (CollAlt), characterizing a contact between the probe and the surface of Titan. Although the simulations took into account the region below 600.0 km, the effects of atmospheric drag in the region were not considered, initially considering only the perturbations of the third body from Saturn and the effects due to the non-sphericity of Titan. Figure 4 shows the maps obtained in these new simulations.
In the maps shown in Figure 4, the regions to the right of the vertical dotted line present the results obtained with the same initial conditions as the maps shown in Figure 3. These initial conditions produced orbits with an initial radius greater than 600 km, thus being initiated outside the assumed region for the presence of atmospheric drag. Considering this initial semi-major axis range and thus the initial altitude of the orbits of probe, the results, as expected, showed an increase in the lifetime of the orbits. This increase is due precisely to considering the dynamics of the orbit of the probe in the region below the value of 600.0 km. These lifetime gains were recorded in both initial circular and eccentric orbits. In the cases of initial eccentric orbits, the gains in lifetimes are in the range of 3 days to 4 months. It was also observed that the most significant gains in lifetimes were found in regions of lower altitude for all simulated cases. Due to the slightest third-body perturbation due to Saturn, the initial circular orbits showed the most significant gains in orbit lifetimes, reaching increments of up to 11 months of orbit. It is noteworthy that, although the collision criterion has been changed, all simulated orbits, which did not survive the total integration time, were classified as collided. Thus, there is no record of escapes from the region of study. When the collision criterion was assumed to be a physical collision (CollAlt), it was possible to obtain some initial orbital conditions between surface of Titan and an altitude of 600 km, the region at the left of the red vertical dotted line. This group of orbits has the most extended lifetimes found for all simulated initial eccentricity conditions among the simulated initial conditions.
The increase in the lifetime of the orbits simulated by the new collision criterion (CollAlt) and the high times obtained in the simulations that started below the altitude of the atmosphere are factors to be observed because they show that a probe can have a specific time of its dynamics in an altitude region lower than 600.0 km. So, while many orbits have a longer orbital flight time above surface of Titan, it is also possible to observe at lower altitudes, at least for a few moments, without imminent loss of the probe. How is the presence of an atmosphere on Titan known? We propose a model for the study of orbits around Titan, Section 3.4, considering the presence of atmospheric drag over the region below 600.0 km altitude, added to third-body perturbations from Saturn and the flattening of Titan.

The Effect Due to Atmospheric Drag
The simulations considering Saturn perturbations and Titan flattening, assuming the probe's physical collision with surface of Titan, Section 3.3, show that there are orbital conditions generated with gains of up to 4.0 months. These results show that, during their orbital evolutions, the probes can stay at an altitude of less than 600.0 km. As we assume that the atmospheric drag in the region below 600.0 km is more intense, simulations were carried out considering an atmospheric drag model added to perturbations of Saturn and the non-sphericity of Titan, aiming to observe how the action of atmospheric drag in these regions can interfere in the lifetime of the orbits. The model for the atmospheric density used in our simulations was based on the studies of [24,25,33]. In [33], three models for the decay of atmospheric density of Titan as a function of the altitude of the probe are proposed, with precision for minimum, maximum, and average curves. Although [33] and other works present a model based on the Cassini-Huygens mission, more recent studies, such as [24], present empirical models that are very similar to others, albeit with certain improvements, but all within the curves established by [33]. Work carried out for this purpose shows that, although these models show the behavior of the atmospheric density until close to 1500.0 km of altitude, lower altitude regions, with 600.0 km, are important, precisely for the study and modeling of orbiter missions that can encounter the atmosphere many times or continuously [24], as shown in the results of Section 3.3. Figure 5 shows the behavior of the atmospheric density model assumed in our work. Analyzing the atmospheric density model presented in Figure 5, simulations were carried out considering the effect of the third body given by Saturn, non-sphericity of the Titan, and atmospheric drag in the region below 600 km of altitude. As well as the drag region, the criterion of physical collision (CollAlt) was considered. It was verified as possible escapes from the orbit of the probe in the region of interest or the loss by collision with the surface of Titan. The probe had a mass of 500.0 kg, a radius of 1.5 m, and a drag coefficient C d = 2.0.
In comparison to the maps obtained considering the perturbations of Saturn, the nonsphericity of Titan, and the criterion of the collision on the surface of Titan without atmospheric drag (see Figure 4), the results taking into account the effects of the atmospheric drag show that the times found in these simulations, for the eccentricities in the range from e 0 = 0.0 to e 0 = 0.2, showed significant differences. Therefore, the orbits registered contributions from the region inferior to 600.0 km smaller than those obtained without the drag effect, where gains in the order of months were registered. This shows that, in these cases, the drag acted such that it favored the loss of the probe, thus reducing the time for carrying out possible missions. However, for orbits implemented with initial eccentricity from e 0 = 0.3 to e 0 = 0.5, the results started to split into two groups: some orbits are still lost faster due to drag, but, in the more eccentric cases, the drag effects in combination with the non-sphericity of Titan and third body perturbations from Saturn, generate orbits with longer durations compared with the maps obtained considering the perturbations of Saturn, non-sphericity of Titan, and the collision criterion on the surface of Titan without atmospheric drag (see Figure 4). These results show the importance of considering the effects of drag at lower altitudes of Titan, as these effects change the scenario of the lifetimes, being an important factor to consider for possible missions.
It is interesting to note that the results obtained in these simulations were, for the range from e 0 = 0.0 to e 0 = 0.4, in most cases, considered close to the results obtained considering the probe loss model at an altitude of 600 km (CollAtm) in Sections 3.1.2 and 3.2. Figure 6 shows the cumulative distribution of duration time of orbits over the drag region in terms of days for circular and eccentric orbits. In terms of the possibility of missions around Titan, the cumulative distribution of duration time of orbits over the drag region is important, because previous work [34] shows possible designs for future missions that aim to carry out studies with space vehicles in the region atmospheric of Titan. Therefore, for these missions, the cumulative distributions presented aim to guide the duration of the probe's orbit, over the drag region, before its loss by collision with the surface of Titan. Contributing to estimating the time and planning of measurement missions and sending data in these regions. Thus, it is concluded that the approach considering the CollAtm collision criterion is, to a certain extent, a precise approach to be used in the study of natural orbits around Titan, as it is a good criterion for a first study around Titan. However, as already mentioned, the drag in the region below 600.0 km in more eccentric orbits, e 0 = 0.3 to e 0 = 0.5 and with higher altitudes, can promote lifetime increases in some cases, it is interesting to use a CollAlt collision criterion. Thus, depending on the type of mission and its interests, it is possible to use one of the criteria we used and obtain accurate results for the lifetimes of the orbits. Another essential point shown in Figure 6 is that the results show the possibility of missions in the final time of the orbits to investigate the atmosphere of Titan, since the decay times of the orbits can reach 12 days, a time that is good enough to complete more observations.

The Effect of Atmospheric Drag in the Region of Lower Atmospheric Density
Based on the results obtained in the simulations, considering the atmospheric drag in the lower region at an altitude of 600.0 km, we investigated how the atmospheric drag present in the less dense region of the atmosphere could interfere with the lifetime of the orbits. For this, we performed simulations considering the effects of the third body due to Saturn, the non-sphericity of Titan, and the atmospheric drag in the altitude region between 600.0 km and 1300.0 km, based on Figure 5. In these simulations, we kept the initial condition grids, the escape criterion based on radius of influence of Titan, and the probe loss criterion when entering the densest region of the atmosphere, CollAtm. The results obtained are shown in Figure 7.  Figure 3, the results in Figure 7 show that the general behavior of the graphs is maintained, even considering the action of atmospheric drag in the region between 600.0 km and 1300.0 km of initial altitude. There is a greater dependence on the choice of eccentricity and initial altitude of the probe and thus on the initial semi-major axis. In orbits with larger eccentricities, the lifetimes are shorter. For fixed initial eccentricities, the value of the initial semi-major axis of the orbits and, therefore, the increase in the initial orbital radius also contributes to shorter lifetimes.

Based on
One change observed was that in the case of simulations considering the Saturn effects, the zonal coefficient J 2 and an atmospheric drag between 600.0 km and 1300.0 km, for all simulated initial eccentricities, orbits with much lower initial altitudes than 1300.0 km show a significant loss in lifetime, evidencing the effects of the atmospheric drag on Titan.
In terms of the orbital lifetimes of the probe, it is observed that initial circular orbits still have the most number of orbits with longer lifetimes, as orbits lasting up to 3 years were found, in the range of 65°≤ I 0 ≤ 120°and orbits with times between 3 and 17 years in the intervals of 50°≤ I 0 ≤ 65°and 120°≤ I 0 ≤ 130°. For e 0 = 0.1 the lifetimes showed a sharp drop compared to circular orbits, with lifetimes of up to 222 days. In terms of the simulations without drag, gains of up to 10 days and losses of approximately 292 days were found.
For e 0 = 0.2, some orbits lasted up to 129 days. Orbits lasting between 17 and 80 days were also found, for e 0 = 0.3. The orbits with the highest losses were found at a 0 = 3.75 R Titan with I 0 = 50°, a 0 = 4.75 R Titan with I 0 = 55°, and a 0 = 3.65 R Titan with I 0 = 50°that recorded losses of 13, 11, and 2 years, respectively. Only one orbit, a 0 = 4.65 R Titan with I 0 = 55°, with a gain of 2 years was found.
In the case of e 0 = 0.4, the orbits, in general, had durations from 1 to 55 days. Having orbits with high lifetimes such as the cases a 0 = 4.45 R Titan with I 0 = 50°and a 0 = 4.35 R Titan with I 0 = 50°, the orbits with lifetimes of 13 and 2 years, respectively. In the case of e 0 = 0.5, the orbits recorded lifetimes in the range of 2 to 32 days. Compared to the simulations without the drag action, they had gains of up to 8 days and losses of only 3 days.
After the study and mapping of orbits as a function of initial inclination, eccentricity, and semi-major axis, the simulations began to investigate the effects, on lifetimes, of the choice of the initial pericenter argument, ω 0 , and node longitude, Ω 0 . To this end, simulations were performed around Titan, considering polar orbits and assuming a fixed value a 0 = 1.85 R Titan for the initial semi-major axis. These simulations were performed for initial eccentricities of e 0 = 0.1, e 0 = 0.2 and e 0 = 0.3. These conditions were maintained for the entire pericenter argument grid and initial node longitude in the range of [0-360] • , always keeping 0 • for the true anomaly. The maps obtained are shown in Figure 8. Figure 8a,c,e present the maps of lifetimes as a function of the pericenter argument and the ascending node longitude in terms of years. Figure 8b,d,f show the maps of the same regions in (a), (c), and (e), respectively, but highlighting other regions, showing the results in terms of days. In general, the maps show that there are regions for the initial ω 0 (pericenter argument) and Ω 0 (node longitude) arranged on the maps in several vertical bands, which give orbits with longer lifetimes. In terms of the initial pericenter argument, these regions are limited to the range of  For the case of a 0 = 1.85 R Titan , I 0 = 90°and e 0 = 0.3, Figure 8e,f, it is found orbits with lifetimes between 7 and 411 days. The orbit lifetime for ω 0 = 0.0°and Ω 0 = 0.0°is 23 days. The durations obtained, varying the value of the node longitude and the pericenter argument, show gains of up to 387 days and a loss of up to 15 days. The longest times are found in isolated regions within the vertical bands in all maps obtained. These results show that, for the simulated region, there is a greater influence on the initial argument of the pericenter of the orbit of the probe than on the longitude of the ascending node. Thus, for orbits closer to the natural satellite, there is an imbalance in relation to the choice of ω 0 and Ω 0 because even though the value of the longitude of the node does not present large differences in the lifetimes, a good choice for the pericenter argument can give considerable gains, thus allowing longer mission durations.

The Evolution of Orbital Elements and the Altitude Range of the Orbits
Considering the lifetimes observed in regions that did not survive the total integration time, a more detailed investigation was carried out to observe the temporal evolution of the orbits. The dynamics showed that, for the cases simulated considering only the perturbation due to Saturn, the eccentricity of the orbits increases as the orbit time increases. Thus, the initially circular orbits tend to become eccentric as Saturn's third-body effects act on the dynamics of the probe. This behavior can be seen in Figure 9, which presents the evolution of some orbits considering the probe losses due to collision with surface of Titan (CollAlt). The same behavior was seen in the cases considering the probe loss in the atmosphere of the Titan (CollAtm). The results also show that, although the eccentricity of the orbits increases with time, there are small oscillations during this period. These oscillations become more evident when the orbits are located at higher altitudes because, in these configurations, the third-body effects produced by Saturn are more intense due to the proximity to the planet. This behavior for the eccentricity of the orbit of the probe is more expressive for the more eccentric initial orbits. Thus, initial circular orbits generate orbits with longer lifetimes is justified. These observations reinforce the utility of circular orbits in the case of missions with orbiters. In terms of the initial inclination, it was seen that within the inclination range between 50°and 90°, the eccentricity of the orbit grows faster. It is the Kozai-Lidov effect that dominates this range of orbits. Previous works have already mentioned this effect on the orbit of probes [11,21,32], and on irregular natural satellites [35]. For the range of 90°< I 0 ≤ 130°, the opposite behavior is observed because, as the probe is implemented in more inclined orbits, the eccentricity grows at a lower rate, taking a long flight, so that the orbit becomes eccentric. Figure 9. Evolution of the inclination and eccentricity of the initial circular orbits not surviving the total integration time. Cases of semi-major axis initial equal to 1.25 R titan (a-c) and 2.65 R titan (b-d).
The evolution of these orbits also showed that, although the inclination has small oscillations (see Figure 9), the orbits tend to maintain their inclination close to the initial value. This behavior is better observed as the initial inclination of the probes approaches 90°. This is an interesting result because, although the orbits tend to decrease lifetimes as the initial inclination approaches 90°, they have durations large enough to make missions, as seen in Figure 3, in addition to keeping their inclination always close to the initial value. Thus, it is possible to carry out missions in polar and quasi-polar orbits, orbital conditions that stand out for carrying out missions with orbiters. Interestingly, this behavior was observed only in the orbits that did not survive the total integration time. The initial conditions surviving 20 years maintained the oscillations in eccentricity and inclination, as shown in figure Figure 10. However, their period of oscillation was considerably smaller compared to the cases that collided. To investigate the effects due to Saturn perturbations, the Titan flattening and the atmospheric drag on the altitude of the orbits, we define the coefficient ∆, which is determined by Equation (9). It determines the intensity of the altitude oscillation amplitude per lifetime of the orbit, based on the average altitude value.
where Alt(t) is the altitude of probe as a function of time, Alt is the average orbital altitude of the probe over the entire orbit time, and ∆t is the duration of the time interval of the orbit of the probe. Figure 11 presents the value of the coefficient ∆ for the dynamics considering only the third-body perturbations due to Saturn, (Figure 3). The results show that the orbits tend to have their altitudes less disturbed when implemented at lower altitudes. This can be seen by the darker regions on the graphs representing the smallest values of the ∆ coefficient. This behavior was observed in all cases of the initial eccentricity simulated and both probe loss regimes. This fact was also observed in the simulations considering only the perturbative effects of Saturn and in those where the non-sphericity of Titan was added. The results also show that as the orbits are implemented with larger initial eccentricities, the ∆ coefficients grow, reinforcing that the more eccentric the orbits, the larger the variations in the altitude of the probe over the lifetime of the orbits. Considering that the initial circular orbits were the ones that registered the largest number of orbits with the longest lifetimes, Figure 11a presents the coefficient ∆ for the initial circular orbits. Given the results observed, it is conclusive that the best region for circular orbits that keep the oscillations in altitude lower are those implemented with values below 3500 km. The dynamics of the orbits also allowed us to observe that even with different altitude oscillation rates, it is a behavior of the amplitude of the orbits to keep their altitudes very close to their initial values until moments close to the half-life of their orbit.  Figure 12 presents the variation of the coefficient ∆ in relation to the simulations considering the coefficient J 2 of Titan and the simulations including only the perturbations of the third body of Saturn. In these maps, it can be seen that the altitudes, in general of the probe, fall into two groups: regions where the action of J 2 dampens the third-body effects and regions where effects due to J 2 accumulate with the third body to increase the oscillations in the altitude of the probes' orbits. These regions can be located on the maps by the negative regions for an orbit in which the Titan's flattening effect decreases the average oscillations in the altitude of the probe and the positive regions on the maps, which correspond to more disturbed orbits, thus generating larger oscillations in the altitude of the orbits. In general, regions with an increase in amplitude of the oscilations are dominant, for all simulated initial eccentricities. Especially in orbits with lower initial altitudes, where the non-sphericity action of Titan is more dominant. However, it is very important to mention that the orbits that have compensation between the third-body perturbation and the J 2 are very important for real missions.
In the case of the maps of the coefficient ∆ for the simulations considering the perturbations of Saturn, the non-sphericity of Titan and the atmospheric drag are shown in Figure 13. The same behaviors can also be observed. There are many orbits with greater amplitudes of oscillation due to the action of atmospheric drag. This is an expected result since atmospheric drag contributes to the variations of the semi-major axis and, consequently, the altitude variation of the probe. Circular orbits were the least affected. These results are important and interesting, as they show the possibility of implementing low-inclination circular natural orbits or even polar orbits around Titan.

Orbital Maneuvers
The results for the lifetimes of the orbits, initially circular, around Titan, considering the perturbations of Saturn, the non-sphericity of Titan, and the atmospheric drag (Figures 7a and 8a), showed the existence of natural polar orbits within the range of [88-786] days. Previous works reinforce the importance of studying and mapping ce-lestial bodies in orbits of this type. In this sense, we propose to carry out orbital maneuvers to relocate these probes at some point before they enter the densest atmospheric region of Titan, in the initial position of their circular orbit. It means that we seek to extend the lifetime of these orbits, thus favoring the accomplishment of missions in these configurations.
We propose to perform the necessary calculations to perform the maneuver at several moments of the orbit's lifetime, trying to find the best moment to perform them. In other words, we calculate the impulses, and thus the fuel cost, for carrying out the maneuver at different times to obtain the moment with the lowest fuel consumption. The calculations were performed based on Equations (10)- (12).
To transfer the probe from its current elliptical orbit to its initial circular orbit, we try to obtain the value of the magnitude of the velocity required to move the spacecraft its current orbit from to the transfer orbit, ∆V 1 , and then the magnitude of the velocity variation required to place it in the circular target orbit, ∆V 2 . R Ap is the radius apocenter of the elliptical orbit before the maneuver, and R cir is the radius of the circular orbit after the maneuver. The chosen orbits and the results of the impulses obtained are presented in Table 3, Figures 14 and 15.    Figures 14 and 15 show the results for the total impulse required per unit of the time to perform the maneuvers to relocate the probes into their initial circular orbits. The orbits in Figure 14 are the ones with the longest lifetimes, and Figure 15 show the orbits with the lowest lifetimes obtained. The results show that as time progresses, the cost of the impulse to carry out the mission increases, showing that it is more efficient and less costly to perform the maneuver moments before the end of the life of the orbit. It is best to perform the maneuvers at times near the end of the orbits' lifetime. However, the maneuvers can also be performed in the middle of the entire life of the orbits, which would already produce an increase of at least 50% in the lifetime of the orbits.
These results are interesting because they can be extended to circular orbits with other initial inclinations, with the goal of favoring the accomplishment of possible missions, increasing the lifetime of orbits, and reducing costs related to propellants for stationrelocating maneuvers.

Final Remarks
In this work, we study the dynamics of natural orbits around the natural satellite Titan. The effects due to Saturn's gravitational attraction were considered, along with the perturbative effects of non-sphericity of Titan, represented by to the gravitational coefficient J 2 . Due to the atmosphere present on Titan, the effects due to atmospheric drag were also considered in this study.
Considering the aforementioned perturbative effects, several grids of initial conditions were simulated, (a 0 , e 0 , I 0 , ω 0 and Ω 0 ), for the orbit of the spacecraft with the objective of generating lifetime maps that would allow us to obtain an understanding of the behavior of orbits and locating the regions with the best conditions for carrying out future missions around Titan. These results reinforce the importance and feasibility of studying natural orbits for space probe missions close to natural satellites and other possible systems. In our work, we presented two criteria for the simulations for a possible loss of the probe that showed to have validity and good precision according to the objectives that can be proposed for studies and simulations. The results obtained in the simulations were presented in terms of lifetime maps that showed that, for all cases of initial eccentricities used, there are intervals for the initial inclination of the orbit of the probe, for several intervals of initial semi-major axis, which generate orbits lasting for the total integration time of 20 years. These conditions lie between the initial inclinations of [0.0-50.0]°and [120.0-130.0]°. For inclinations between [50.0-120.0]°, the orbits did not survive the simulation time used, producing orbits with varying durations between a few days and less than 20 years. The results also showed that the orbits with longer durations were those implemented with initial inclinations closer to the lower and upper limits in the range of [50.0-120.0]°. In terms of the initial altitude of the orbits, the results showed that the lower the initial altitude of the orbit, the longer the duration of the orbits. This behavior was observed for all simulated initial eccentricities.
The simulations considering the effect of the flattening of Titan, J 2 , showed that, for certain initial orbital conditions, it is possible to have a balance between the third-body effects from Saturn and the effects due to the J 2 of Titan. This balance tends to generate orbits with longer lifetimes. However, this was not the dominant behavior observed. Instead, many of the orbits maintained their lifetimes very close to the results considering only the third-body effects, which shows the dominance of the effect on the dynamics of orbits around Titan.
In the study of the effects due to the atmospheric drag present in Titan, it was observed that, in the two regions defined in our work, the drag acted, in most cases, to favor the reduction of the orbit lifetime and thus contribute to the loss of probe, reducing the durations of possible missions. These results show the importance of considering the effects of drag at the lower altitudes of Titan, as these effects change the scenario of the duration of the orbits, being an essential factor to consider for possible missions.
An imbalance was also seen in relation to the choice of ω 0 and Ω 0 because, even if the value of the longitude of the node does not present much effect in the lifetimes, a good choice for the pericenter argument can give considerable gains, thus allowing longer mission durations.
The investigation of the oscillation rate of the probe altitude, coefficient ∆, showed that the orbits are implemented with higher initial eccentricities. These coefficients grow, reinforcing that the more eccentric orbits present more significant variations in the probe altitude and the orbital lifetime of the probe. In the simulations considering the coefficient J 2 of Titan, it was observed that the altitudes of the probe, in general, are divided into two groups: regions where the action of J 2 dampens third-body effects and regions where effects due to J 2 contribute with the third body to increase the oscillations in the altitude of the probes' orbits. In the case of the maps of the coefficient ∆, for the simulations considering the perturbations of Saturn, the non-sphericity of Titan, and the atmospheric drag, the same behaviors were observed. There are many orbits with a larger amplitude of oscillation due to the action of atmospheric drag. Circular orbits were the least affected.
These results are significant and exciting, as they show the possibility of implementing natural low-inclination circular orbits or even polar orbits around Titan.
We also proposed to carry out maneuvers for some initial polar orbits to find the best moment to apply them. This aims to minimize fuel costs to return the probe to its initial configuration, thus increasing the duration of the orbits. These results are interesting because they can be extended to circular orbits at other initial inclinations, to favor the accomplishment of possible missions, increasing the lifetime of the orbits, and reducing costs related to propellants. All the results presented in this work can help plan missions to Titan, thus contributing to the achievement of the scientific objectives already mentioned in this work.