Analytical solution to logarithmic spiral trajectories with circumferential thrust and mission applications

This study made use of a shape-based method to analyze the orbital dynamics of a spacecraft subject to a continuous propulsive acceleration acting along the circumferential direction. Under the assumption of a logarithmic spiral trajectory, an exact solution to the equations of motion exists, which allows the spacecraft state variables and flight time to be expressed as a function of the angular coordinate. There is also a case characterized by specific initial conditions in which the time evolution of the state variables may be analytically determined. In this context, the presented solution is used to analyze circle-to-circle trajectories, where the combination of two impulsive maneuvers and a logarithmic spiral path are used to accomplish the transfer. The determined results are then applied to the achievement of the Earth—Mars and the Earth—Venus transfers using actual data from a recent thruster developed by NASA.


Introduction
Analytical solutions to the differential equations describing the motion of a spacecraft subject to a continuous propulsive acceleration [1,2] are known only in a few cases [3][4][5]. Nevertheless, when available, closed-form results represent a very useful tool in the preliminary phase of mission design because they simplify the parametric study of spacecraft trajectories and significantly reduce the computational cost of numerical simulations. In practice, analytical solutions may be derived using two different approaches. The first one involves solving the equations of motion for a prescribed thrust profile [6][7][8][9][10][11][12][13][14][15][16], while the second one consists in using a shape-based approach in which the thrust profile is calculated in such a way that a prescribed trajectory may be actually covered by the spacecraft [17][18][19][20][21][22][23][24]. A brief overview of the applications and a summary of the main results achievable with these two approaches were provided by Bassetto et al. [24], which the interested reader is invited to consider.
Within the context of the shape-based methods, the aim of this study was to analyze the generation of logarithmic spiral trajectories using a purely circumferential propulsive acceleration. In particular, we demonstrate the existence of an exact solution to the equations of motion that allows the spacecraft state variables to be explicitly written as a function of the angular coordinate. It is also shown that the flight time takes the form of a hypergeometric function of the polar angle and that the state variables may be explicitly obtained as a function of time under the assumption of suitable initial conditions. Finally, the propulsive acceleration required for the spacecraft to trace the desired logarithmic spiral trajectory is analytically derived.
The remainder of this paper is organized as follows. The first part of Section 2 focuses on trajectory analysis and presents the analytical solution to the polar equations of motion assuming that the spacecraft follows a logarithmic spiral path under the action of a circumferential propulsive acceleration. The second part of Section 2 derives the magnitude of the propulsive acceleration required for the spacecraft to trace the prescribed logarithmic spiral trajectory. Section 3 discusses a case study, in which the assumption is made that the spacecraft performs a circle-to-circle orbital transfer Nomenclature a Semi-major axis (km) a Propulsive acceleration vector (mm/s 2 ) a θ Circumferential component of a (mm/s 2 ) {a2, b1} Functions of α, see Eq. (20) {c1, c2} Constants of integration, see Eqs. (15) and (23)  by combining a logarithmic spiral trajectory and two impulsive maneuvers. In this context, a few cases of interplanetary transfers are discussed in detail by analyzing both the required total velocity variation and flight time. Section 4 summarizes the main results of this study and provides concluding remarks.

Mathematical model
Assume that at a given initial time t 0 ≜ 0 a spacecraft S is tracing a Keplerian parking orbit (with eccentricity e 0 ) around a primary body with a gravitational parameter µ. We introduce a polar reference frame T (O;r,θ), where the origin O coincides with the center of mass of the primary body andr (orθ) is the radial (or circumferential) unit vector (see Fig. 1). The spacecraft state is described by a set of four variables {r, θ, v r , v θ }, where r is the distance of the spacecraft from O, θ is the polar angle (measured in the counterclockwise direction from the primary-spacecraft line at time t 0 ), and v r (or v θ ) is the radial (or circumferential) component of the spacecraft velocity vector v. The spacecraft is assumed to be equipped with a primary propulsion system that provides a continuous and purely circumferential propulsive acceleration a of adjustable magnitude |a θ |, that is The two-dimensional spacecraft dynamics is described by the following nonlinear differential equations: which are completed by the initial conditions (subscript 0): Note that {v r0 , v θ0 } can be written as a function of r 0 and the initial spacecraft true anomaly

Conditions for tracing a logarithmic spiral trajectory
We now investigate the conditions under which the propelled trajectory is a logarithmic spiral in the form [17]: In this case, the constant parameter α ̸ = 0 is easily shown to coincide with the tangent of the flight path angle γ ≜ arctan(v r /v θ ) (see Fig. 1). In fact, from Eqs. (2) and (3) we have v r ≡ṙ = αr 0θ exp(αθ) v θ ≡ rθ = r 0θ exp(αθ) Clearly, the sign of α coincides with that of v r , that is, α > 0 (or α < 0) corresponds to an orbit raising (or lowering). Substituting Eq. (9) into Eq. (4), the following second-order nonlinear differential equation is obtained: where Ω ≜ µ/r 3 0 is a reference angular rate. The sign ofθ depends on the travel direction along the logarithmic spiral in the plane of the parking orbit. In particular, θ > 0 (orθ < 0) indicates that the spiral trajectory is traced counterclockwise (or clockwise). Without loss of generality, the case ofθ > 0 is assumed in the rest of this paper. Accordingly, from Eq. (12) we obtaiṅ where the auxiliary function depends on θ and contains a constant of integration c 1 . The latter is obtained by substituting the initial conditions Eq. (6) into Eqs. (13) and (14), viz.
The evolution of the state variables v r and v θ can now be obtained as a function of θ by substituting Eq. (13) into Eq. (10), that is Additionally, it should be noted that Eq. (13) gives feasible solutions provided z ⩽ 1. With the aid of Eq. (14), the latter inequality translates into the following conditions on the polar angle: where z 0 is the value of z when θ = θ(t 0 ) = 0, that is An analytical relationship exists between the time t and the polar angle θ, which may be determined by integrating Eq. (13) with the initial conditions Eq. (6). The result is c 2 is a constant of integration, and 2 F 1 (1, a 2 ; b 1 ; z) is the hypergeometric function [25], defined as where the coefficients λ n may be written in a recursive manner as In Eq. (19), the constant of integration c 2 is determined by enforcing the initial condition θ(t 0 ) = 0, viz.
where z 0 is given by Eq. (18). Figure 2 shows the variations of c 1 and c 2 as a function of ν 0 and for e 0 = {0.01, 0.05, 0.1}, assuming that the spacecraft is initially placed on a heliocentric parking orbit of period of 1 yr. It is worth noting that, in general, Eq. (19) does not allow the polar angle (nor any spacecraft state variable) to be explicitly written as a function of time. However, the functions r = r(t), θ = θ(t), v r = v r (t), and v θ = v θ (t) can be obtained in a closed form when c 1 is equal to zero. This problem is discussed in the following.
Spacecraft state variables as a function of time: case of c 1 = 0 Consider the special case in which as it follows from Eq. (15). For a given value of the initial eccentricity e 0 , the condition c 1 = 0 amounts to setting the initial position of the spacecraft along the parking orbit. In fact, recall that The last two expressions are drawn in Fig. 3 for e 0 ∈ (0, 1). Note that cos ν 0 < 0, while the sign of sin ν 0 is the same as that of α. In summary, when c 1 = 0, the value of ν 0 can be written as a function of e 0 ∈ (0, 1) as which proves that the initial position of the spacecraft along the parking orbit is univocally assigned. This result is consistent with Fig. 2(a), which shows that there are two values of ν 0 (one for the ascending phase and the other one for the descending phase) for which c 1 = 0. Finally, since α = 0 when e 0 = 0, in this case the spacecraft remains in the circular parking orbit and r(θ) ≡ r 0 for any θ ⩾ 0. Now, consider the problem of calculating how the spacecraft state variables {r, θ, v r , v θ } and the classical orbital parameters {a, e, ν, ω} vary over time, where a and ω represent the semi-major axis and argument of periapsis of the osculating orbit, respectively. According to Eq. (14), when c 1 = 0, then z = 0. In fact, substituting Eq. (9) into Eq. (13) yieldṡ which may be easily integrated to get the time variation of the polar angle as where is a dimensionless function of time. Note that Eq. (31) gives feasible values of θ for any t ∈ [0, +∞) when α > 0 and for t ∈ [0, − √ 2α 2 + 4/(3αΩ)) when α < 0. Using Eq. (31), the time variation of the spacecraft distance from the primary can be obtained from Eq. (9) as Similarly, the radial (v r ) and circumferential (v θ ) components of v are obtained from Eqs. (2) and (3) and from Eqs.
Notably, Eqs. (9), (34), and (35) can be combined to obtain useful dimensionless expressions for v r and v θ as The time variation of the semi-major axis a and eccentricity e of the osculating orbit can be easily calculated with the aid of Eqs.
from which the eccentricity of the osculating orbit turns out to be constant (i.e., e(t) ≡ e 0 ), while r(t)/a(t) ≡ 1 − e 2 0 . An expression for cos ν(t) can be derived from the polar form of the spacecraft osculating orbit as The latter expression indicates that the true anomaly of the spacecraft on the osculating orbit remains constant (i.e., ν(t) = ν 0 ). This means that the increment in θ corresponds to the rotation of the apse line of the osculating orbit. Therefore, the time variation of ω is where ω 0 ≜ ω(t 0 ) is the initial value of ω. Finally, it should be noted that the general solution (i.e., the solution obtained with c 1 ̸ = 0) may be approximated from that corresponding to c 1 = 0 when z is close to zero. The latter condition occurs when α < 0 and the exponential in Eq. (14) becomes sufficiently small as θ increases.

Propulsive acceleration analysis
This section is aimed at calculating the circumferential propulsive acceleration a θ that is necessary for a spacecraft to trace a logarithmic spiral trajectory of given characteristics. In the general case in which c 1 ̸ = 0, a θ can be easily derived from Eqs.
When substituting Eqs. (13) and (14) into Eq. (42), the required acceleration is For increasing values of θ, the exponential term in Eq. (43) diverges when α < − √ 2 or α ∈ (0, √ 2), whereas it is equal to zero when |α| = √ 2. The functions Eqs. (31) and (33)-(35) can be effectively used to evaluate a θ when c 1 = 0. In this case, Eq. (43) reduces to which points out that the necessary acceleration is inversely proportional to the square of the orbital radius and its sign agrees with that of α.

Comparison to literature results
It is noteworthy that the special case in which c 1 = 0 can also be derived from the results obtained by Bacon in 1959 [17], who followed the pioneering work by Forbes [26]. In Ref. [17], Bacon identified the logarithmic spiral as a special solution to the equations of motion of a spacecraft subject to the gravitational attraction of a primary body and to a continuous thrust of a suitably adjustable magnitude. According to Bacon [17], when a spacecraft travels a logarithmic spiral trajectory, the propulsive acceleration vector is always tangent to the velocity vector and its magnitude is inversely proportional to the square of the orbital radius. When the equations of motion in Bacon's paper [17] are adapted to the nomenclature of this work, we obtain the following pair of differential equations: As a result, the propulsive acceleration magnitude required to generate a Bacon's logarithmic spiral (subscript B) is Equations (45) and (46) can be equivalently rewritten by incorporating the radial component of the propulsive acceleration into the gravitational term through a suitable variation of µ. More precisely, denoting with the gravitational parameterμ that includes the contribution of the radial component of the propulsive acceleration, that is, definingμ then Eqs. (45) and (46) becomė which are formally equivalent to Eqs. (4) and (5) when a θ is defined as in Eq. (44). In other words, the logarithmic spiral trajectories investigated by Bacon [17] can also be described using Eqs. (31) and (33) by simply replacing µ (which is implicitly contained in Ω) withμ. The value of |a θ | in the case of a purely circumferential thrust is smaller than the propulsive acceleration ∥a B ∥ predicted by Bacon's model [17]. In fact, according to Eqs. (44) and (47), we can find that |a θ | / ∥a B ∥ = 2 √ α 2 + 1/(α 2 + 2) < 1 for any |α| ̸ = 0. Therefore, although the two trajectories have the same logarithmic spiral shape, they require two different levels of propulsive acceleration to be maintained. Furthermore, the logarithmic spirals obtained with a circumferential thrust and c 1 = 0 are traced at a higher angular rate. In fact, differentiating Eq. (31) with respect to time yieldṡ while the logarithmic spirals analyzed by Bacon [17] are traveled at an angular rateθ B given bẏ from which it is clear thatθ >θ B for any α ̸ = 0 and t ⩾ 0. According to McInnes [19], Eqs. (45) and (46) describe the two-dimensional heliocentric motion of an ideal solar sail [23] with a lightness number β and a pitch angle α n ∈ (−π/2, π/2) rad given by which are plotted in Fig. 4 for α ∈ [−5, 5]. Recall that β is the ratio of the maximum allowable thrust magnitude to the local weight and that α n describes the acute angle between the radial unit vector and the normal to the sail plane [19,22]. Therefore, the solar sail-based logarithmic spirals discussed by McInnes [19] and recently analyzed n Fig. 4 Variations of β and αn with α.
by Bassetto et al. [22] are different from those described by Eqs. (31) and (33) when µ is replaced withμ. In fact, despite having the same logarithmic spiral shape, they are traced at different angular rates except when the special conditions in Eq. (53) are met.

Mission application: Two-dimensional circle-to-circle transfers
This section discusses a possible mission application of the previously discussed model when c 1 = 0 and the parking and target orbits are both circular (with radii r 0 and r f ̸ = r 0 , respectively) and coplanar. A circle-to-circle transfer requires an entry phase into the logarithmic spiral path from the circular parking orbit and an exit phase from the logarithmic spiral path towards the circular target orbit. This is because Eq. (28) gives α = 0 when e 0 = 0, which means that the parking orbit and the target one cannot be tangent to any logarithmic spiral trajectory. In other words, a logarithmic spiral transfer between two circular coplanar orbits is unfeasible without the introduction of two "transition phases", namely a parking-to-spiral (P2S) and a spiral-to-target (S2T) phase.

Transition phases with impulsive maneuvers
For a given value of α ̸ = 0, assume that the P2S and the S2T phases are obtained with an impulsive maneuver applied at time t 0 and t f > t 0 , respectively [27]. Clearly, t f represents the overall duration of the circle-to-circle transfer, which coincides with the flight time along the logarithmic spiral trajectory. The value of t f can be obtained from Eq. (33) by enforcing the final condition r(t f ) = r f , that is A value of ρ > 1 (or ρ < 1) corresponds to an orbit raising (or lowering) [28], which is obtained through a logarithmic spiral trajectory with a value of α > 0 (or α < 0). For this reason, the right-hand side of Eq. (54) is always positive. The variation in the overall flight time with {α, ρ} is shown in Fig. 5, where T 0 ≜ 2π/Ω is the orbital period of the circular parking orbit. Since the initial spacecraft velocity vector is circumferential and its magnitude is µ/r 0 , the impulsive velocity change ∆v P2S required to enter the logarithmic spiral at time t 0 (phase P2S) is where v θ (t 0 ) is obtained from Eq. (36) by substituting the initial condition r(t 0 ) = r 0 , viz.
v θ (t 0 ) = 2 µ/r 0 √ 2α 2 + 4 (57) By substituting Eq. (57) into Eq. (56) we get ∆v P2S = µ r 0 The impulsive velocity change ∆v S2T required to accomplish the S2T phase, that is, to circularize the orbit at the end of the logarithmic spiral path (i.e., at time t f ) is Therefore, the expression of ∆v S2T is Finally, bearing in mind Eq. (44), the velocity variation ∆v S along the logarithmic spiral trajectory can be written as in accordance with the classical results of Alfano and Thorne [29] and Vallado [30]. Similar results were reported in a previous study [31] for solar sail-based missions. Therefore, the total velocity variation along the whole transfer is the sum of three contributions, that is The final result can be rewritten in a dimensionless form using Eqs. (55), (58), (61), and (62) as ∆v tot which is illustrated in Fig. 6. Equations (54) and (64) can be used to obtain the flight time and the total velocity variation in a circle-to-circle transfer when the spacecraft moves along a logarithmic spiral arc. This issue is discussed in the following section.

Cases of Earth-Mars and Earth-Venus transfers
Consider the heliocentric phases of ephemeris-free Earth- Mars and Earth-Venus transfers [32] assuming that the eccentricity and the mutual inclination of the two planetary orbits are both neglected. The heliocentric phases of the interplanetary transfers can therefore be approximated with a circle-to-circle orbit raising in the Earth-Mars case (subscript ♂, with ρ = ρ ♂ ≜ 1.5237) and a circle-to-circle orbit lowering in the Earth-Venus case (subscript ♀, with ρ = ρ ♀ ≜ 0.7233). Consider a logarithmic spiral-based transfer with a given value of α, where α > 0 in the Earth-Mars case and α < 0 in the Earth-Venus case. According to Eq. (44), in the Earth-Mars case, the maximum and minimum values of |a θ | are reached at the beginning and at the end of the transfer, respectively, that is where µ/r 2 0 ≃ 5.93 mm/s 2 . In the Earth-Venus case, instead, Eq. (44) gives The maximum change in the propulsive acceleration magnitude, that is, the difference ∆a θ ≜ max |a θ | − min |a θ | is The above equations, along with Eqs. (54) and (64) where Ω = 2π rad/yr and µ/r 0 ≃ 29.7853 km/s, can be used to evaluate the performance of the interplanetary transfer as a function of the geometrical parameter α that characterizes the logarithmic spiral shape. For the two previous mission cases, Figs. 7 and 8 show the interplanetary transfer performance in terms of dimensionless velocity variation, flight time, and maximum magnitude of the circumferential propulsive acceleration. Instead, Fig. 9 shows the variation of ∆a θ ♂ and ∆a θ ♀ as a function of α. For exemplary purposes, in both cases it is assumed that the entire mission is to be conducted within a time interval t f = 4 yr. The previously defined mathematical model provides the results summarized in Table 1 and Figs. 7-9. The corresponding transfer trajectories are plotted in Fig. 10, whereas Fig. 11 reports the time variation of a θ .
It is reasonable to ask whether the mission scenario described in Table 1 and in Figs. 10 and 11 is practically feasible for a spacecraft equipped with a real engine. In particular, we address this problem by assuming that the spacecraft uses an electric thruster. It is well known that an electric thruster has a finite number of operation points [33,34], each corresponding to a certain thrust level. For this reason, it is impossible to generate a continuously variable propulsive acceleration, as shown in Fig. 11. However, the continuous thrust variation can be replaced by a succession of discrete thrust levels that provide the required propulsive acceleration on average. For example, the NASA Evolutionary Xenon Thruster (NEXT) [34] has a total of 40 operation points I d , as illustrated in Fig. 12, with a minimum available thrust    Figure 12 also points out that the specific impulse, which ranges between 1,406 and 4,324 s, is always greater than 3,000 s within the first 33 operation points with a mean value equal toĪ sp ≜ 3,669 s. Additionally, the available thrust is always greater than 60 mN within the first 33 operation     points. To check the mission feasibility, it is necessary to calculate the maximum and minimum required thrust and verify whether these values fall within the range of permissible thrust levels.
The required thrust is equal to the product of the current spacecraft mass m and the required propulsive acceleration |a θ |, where the latter can be calculated using Eqs. (33) and (44). On the other hand, the thrust can also be thought of as the product of the mass flow rate −ṁ > 0 and the exhaust velocity relative to the spacecraft u ≜ g 0 I sp > 0, where g 0 ≜ 9.80665 m/s 2 is the standard gravity [36]. Therefore, the required thrust is obtained through the Eq. (71): where a spacecraft mass of m 0 ≜ 3,400 kg (similar to that of the NASA space probe Juno) is assumed after the P2S phase (i.e., at the beginning of the logarithmic spiral path). Equation (71) can be remarkably simplified by approximating I sp with its mean value calculated over the first 33 operation points (i.e., by setting I sp ≡Ī sp ). The integral in Eq. (71) can now be easily solved and the required thrust becomes which is shown in Fig. 13 as a function of time for the two mission cases. Table 2 shows the results in terms of maximum and minimum required thrusts and spacecraft mass just before the S2T phase for the two mission cases. One can see that the maximum and minimum thrust levels fall within the first 33 operation points of the NEXT ion engine. Finally, it should be noted that the consumed mass fraction during the logarithmic spiral path is only 14.5% in the Earth-Mars scenario and 13.5% in the Earth-Venus scenario.
For the sake of completeness, Figs. 14 and 15 show the Earth-Mars and the Earth-Venus transfer trajectories when t f = {1, 2, 3, 5} yr. Clearly, depending on the specific mission scenario, it will be necessary to select a suitable initial spacecraft mass and type of engine to ensure that the latter is able to provide (at least on average) the required propulsive acceleration along

Conclusions
This paper proposed a novel analytical solution to the two-dimensional equations of motion of a propelled spacecraft under the assumption that the spacecraft traces a logarithmic spiral trajectory by means of a continuous and circumferential thrust of adjustable magnitude. In the most general case, the flight time takes the form of a hypergeometric function of the angular coordinate, while the spacecraft state variables and the required propulsive acceleration are obtained as a function of the polar angle. When the initial orbit is closed, there is always a particular position along it starting from which all the state variables can be written as explicit functions of time. We used such functions to analyze circle-to-circle orbital transfers by assuming that the spacecraft exited from and entered the logarithmic spiral path using an impulsive maneuver. The proposed model was also applied to the study of simplified Earth-Mars and Earth-Venus transfers using realistic data from a

Funding note
Open Access funding provided by University of Pisa.