Analytical solutions for low-thrust orbit transfers

This paper presents analytical solutions for the estimation of the ΔV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta V$$\end{document} cost of the transfer of a spacecraft subject to a low-thrust action. The equations represent an extension of solutions already available in the literature. Moreover, the paper presents novel analytical solutions for low-thrust transfers under the effect of the second-order zonal harmonics of the Earth’s gravitational potential. In particular, the paper is divided into two parts. The first part presents analytical expressions for the ΔV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta V$$\end{document} cost of transfers. All analytical equations were validated through numerical integration of the dynamics of the spacecraft. The second part of the paper introduces new analytical equations for low-thrust transfers between circular inclined orbits with different values of the right ascension of the ascending node, under the effect of the second-order zonal harmonic of the Earth’s gravitational potential. Both in the first and second parts, analytic solutions for the variation with time of the orbital elements during the transfer are presented. The proposed equations are applicable to low-thrust transfer realised through a long spiral trajectory.


Introduction
In the early stage of the definition of a space mission, it is often desirable to have a fast preliminary estimation of the V cost of a low-thrust transfer. In particular, when the transfer is realised through a long spiral trajectory, a quick estimation of the total V and transfer time can avoid potentially lengthy, albeit more accurate, calculations. For this reason a number of authors have proposed simple control laws for the variation of specific orbital elements, and, in some cases, analytical equations for the estimation of the V associated to a given type of transfer. Burt (1967) presented the secular rates of change of the orbital elements for a thrust action having radial, transversal and normal components, constant in magnitude over one revolution. Burt showed that by reversing the thrust directions at certain points in the orbit, it is possible to change the orbital elements independently one from the others.
In the work presented by Burt (1967), no closed-form expressions of the V are provided. In addition, the integral of the secular rates of changes of the orbital elements was derived only for one of the proposed control laws. Edelbaum presented a solution for the minimumtime low-thrust transfer between inclined circular orbits (Edelbaum 1961). This solution was later reformulated by Kechichian (2010), where it was further extended in order to constrain the intermediate orbits during the transfer to remain below a given altitude. In addition, Kechichian presented the solution of the problem of a transfer between circular orbits with different values of the right ascension of the ascending node. In the work presented by Kéchichian (2010), the author derived some analytical expressions for the V and for the variation of the orbital elements. Gao presented an analytic orbital averaging technique to obtain incremental changes of semi-major axis, eccentricity and inclination over an orbital revolution for specific control laws (Gao 2007). However, no closed-form expression of the V was derived. In the work of Gao and Kluever (2005), the authors use approximated integrands in order to obtain an analytical solution for the variation of the orbital elements. Also in this case, no simple equation for the V was derived. Hudson presented a method where the components of the thrust acceleration are represented by means of Fourier series in the eccentric anomaly (Hudson and Scheeres 2009). Analytic expressions for the secular rate of change of the orbital elements were given, but no simple equation for the V or for the variation of the orbital elements with time was derived. Petropolous studied specific control laws for escape trajectories; he derived equations relating the average energy and eccentricity during the transfer, and for the approximated estimation of the V (Petropolous 2002). Later, Petropoulos (2003) presented methods to determine the thrust direction to obtain specific changes in the orbital elements, leading to the definition of simple control laws. Simple control laws were also presented in Ruggiero et al. (2011). However, neither Petropoulos (2003) nor Ruggiero et al. (2011) report any analytical formula for the quick estimation of the V or for the variation of the orbital elements. Finally, Pollard (2000) presented different pitch steering cases and their effect on the variation of the orbital elements. Of particular relevance, among the profiles presented by Pollard, are a steering profile that produces the simultaneous variation of the eccentricity and inclination, and one that produces a variation of the argument of the periapsis while holding semi-major axis and eccentricity constant. The control laws were applied on thrust arcs centred at pericentre and apocentre. For this specific case, an expression for the V cost of the transfer is reported in Pollard (2000), but no equation is given for the variation of the orbital elements during the transfer.
Since not all the results available in the literature include analytical expressions for the V to realise a transfer, or for the variation of the orbital elements, the first part of this paper extends the literature by introducing a number of analytical formulas for the quick estimation of the V and for the evolution of the associated mean orbital elements. The novel equations presented in the first part of the paper have been obtained by extending the work of Burt (1967), Petropoulos (2003), Pollard (2000) and Ruggiero et al. (2011). For the control laws defined by Burt (1967), we now present analytical expression for the V and for the variation of the orbital elements during the transfer. In addition, we extend the analytical expression for the V given in Pollard (2000). While Pollard considered thrust applied on arcs centred at the apocentre and pericentre of an orbit, we now consider the same control law applied along the entire orbit. Finally, we have included the analytical expressions for the V and for the variation of the orbital elements during the transfer, for some of the control laws presented by Petropoulos (2003) and Ruggiero et al. (2011) and for the case of circular orbits. While the time variation of the orbital parameters is presented for each considered control law, it is important to stress that the main focus of the paper remains the quick estimation of the V s through simple analytical formulas. The V s obtained using the analytical equations For selected solution, a detailed analysis of the V error, and of the error associated to the equations for the evolution of the mean orbital elements, is presented. When different formulas, associated to transfers characterised by different thrust profiles, can be used to define the same variation of orbital elements, we provide a comparative analysis of the V cost of these formulas.
The solutions to the low-thrust motion presented in the first part of the paper do not consider perturbations to the orbital motion. However, gravity perturbations, such as the second zonal harmonic of the Earth's gravitational potential, J 2 , cause secular variations of the orbital elements of low Earth's orbits. Thus, in this paper, new analytical solutions are derived for the two-point boundary value problem (TPBVP) associated to low-thrust transfers between inclined circular orbits with different values of the right ascension of the ascending node, under the effect of J 2 . Analytical laws are derived for the variation of the orbital elements with time and for the cost of the transfer. The underlying assumption is that the thrust level is low, so that the transfer is represented by a long spiral.
The paper is structured as follows: Sect. 2 presents the relevant equations of motion; in the first part of the paper, Sects. 3 to 7 present the extensions to the available solutions for the variation of the semi-major axis (or combination of semi-major axis and other orbital elements), eccentricity (or combination of eccentricity and other orbital elements), inclination, right ascension of the ascending node, and argument of the periapsis. The analytical solutions presented in Sects. 3 to 7 are summarised at the end of the paper in "Appendix A", where a taxonomy of the analytical solutions is also given. The second part of the paper is constituted by Sect. 8, which introduces new analytical solutions for the transfer between circular inclined orbits with different values of the right ascension, under the effect of J 2 . To the authors' best knowledge, similar solutions are not available in the literature. Finally, Sect. 9 concludes the paper.

Equations of motion
It is assumed that the spacecraft is subject to a low-thrust acceleration f that can be expressed in a radial-transverse-normal RTN reference frame (see Section 1.4.3., Vallado 2007 and Fig. 1) as In Eq. (1) is the magnitude of the low-thrust acceleration, assumed to be constant, α is the azimuth angle, with respect to the transverse direction, and β is the elevation angle with respect to the orbit plane. Gauss' equations, expressing the time variation of the osculating classic orbital elements due to the low-thrust acceleration f, are (see Section 10.3, page 484, Battin 1987) In Eq. (2), a is the semi-major axis, e the eccentricity, p the semilatus rectum, i the inclination, Ω the right ascension of the ascending node, ω the argument of the periapsis, θ the true anomaly, h the magnitude of the angular momentum, h = √ μ p, and ν the argument of The approximation in the last of Eq.
(2) derives from the assumption that the perturbing forces are small enough to produce a negligible effect on dθ/dt, compared to the term h/r 2 , which is due to the Keplerian motion (a similar approximation is used also in Avanzini et al. 2015;Di Carlo et al. 2021;Zuiani and Vasile 2015). The variation with time of the argument of latitude is (Section 10.3, page 484, Battin 1987) The approximation in Eq.
(3) is based on the same assumptions used for the approximation in the last of Eq. (2). Another useful expression, used in the rest of this paper, is the time variation of the eccentric anomaly E (Burt 1967): The approximation in Eq. (4) holds when da/dt and de/dt are small (Burt 1967). In the rest of the paper, it is in fact assumed, unless otherwise specified, that the perturbing forces are small enough to produce negligible variations in the orbital elements over one orbital period.

Variation of the semi-major axis
In this section, we present two analytical solutions for the variation of the semi-major axis: the variation of the semi-major axis with constant tangential thrust and the variation of the semi-major axis with constant radial thrust. Both of them present no long-term variation of the eccentricity. Other solutions for the variation of the semi-major axis were presented by Edelbaum (1961) and Kéchichian (1997Kéchichian ( , 2010. In particular, Edelbaum (1961) and Kéchichian (1997Kéchichian ( , 2010 present solutions for the simultaneous variation of the semi-major axis and inclination, and of the semi-major axis and right ascension of the ascending node.

Variation of the semi-major axis with tangential thrust
The thrust angles that provide the maximum instantaneous rate of change of the orbital elements can be derived from the first derivatives of the Gauss' equations with respect to α and β, as presented by Ruggiero et al. (2011). In particular, for the maximum variation of the semi-major axis: Equation (5) provides the well-known result that the rate of change of the semi-major axis is maximum for planar thrust (β = 0) and azimuth angle equal to the flight path angle γ : In particular, α = γ to increase the semi-major axis and α = π + γ to decrease the semimajor axis. Equations (5) and (6) are given by Ruggiero et al. (2011), while, to the authors' best knowledge, the equations and derivations presented in the rest of this subsection have not been presented anywhere else in the literature. Using Eq. (6), in the general case of non-circular orbits, the time variations of semi-major axis, eccentricity and argument of periapsis can be written as: Equations (7) are derived from Eq. (2) by expressing sin α and cos α as cos α = sgn(cos α) 1 + e cos θ √ 1 + e 2 + 2e cos θ , sin α = sgn(cos α) e sin θ √ 1 + e 2 + 2e cos θ .
The variations of a, e and ω with θ are thus given by We now introduce the assumption that the variations of a, e and ω are small within one revolution, since the thrust level is assumed to be small (Petropolous 2002). Therefore, they can be kept constant when computing the average quantities: whereā,ē andω are the mean values of the semi-major axis, eccentricity and argument of periapsis. The last two integrals in Eqs. (10) give ē 2π = ω 2π = 0, meaning that there is no secular variation of the eccentricity and argument of periapsis. The equation for the semi-major axis gives instead The expression for f (ē) is In Eq. (12) F Ic and E Ic are, respectively, the complete elliptic integrals of the first and second kind: and F I and E I are the elliptic integrals of the first and second kind (Section 1.5, page 68, Battin 1987). The value of the elliptic integrals was computed with the MATLAB function ellipke which implements the method of the arithmetic-geometric mean described by Abramowitz and Stegun (1965), Section 17.6. The tolerance default value of 2.2204 × 10 −16 was used. The variation with time of the semi-major axis is and depends on the eccentricityē. In Eq. (14), T denotes the orbital period, T = 2π ā 3 /μ. Since there is no secular variation of the eccentricity during the transfer, the secular variation of the semi-major axis with time can be obtained by integrating Eq. (14) to have: Finally, integrating Eq. (14) from t 0 = 0 to ToF, withā(t 0 ) =ā 0 andā(ToF) =ā f , it is possible to state that: Proposition 1 Given Eq. (15), the cost of a transfer to change the semi-major axis fromā 0 tō a f , using the thrust profile defined in Eq. (6), and under the assumption that the variation of semi-major axis, eccentricity and argument of the periapsis is small within one revolution, is Equation (16) is not singular for circular orbits, since whenē = 0, f (ē) = 0. It is important to stress that Eqs. (16) and (15) provide only an approximation to the V and to the variation of a during the transfer, and that the approximation is incorrect when the variation of the semimajor axis over one orbital period, or the variation of the eccentricity, cannot be neglected. The validity of this assumption is illustrated hereafter for Geocentric and Heliocentric, or interplanetary, transfers. Figure 2 shows the relative difference in the V obtained using  analytical Eq. (16) and computing the cost of the transfer by numerically integrating Gauss' equations with Eq. (6), until the desired final semi-major axis is reached. More precisely, Fig. 2 shows the quantity ( V − V num )/ V num , where V is given by Eq. (16) and V num is computed as V num = ToF num . The time of flight ToF num is obtained by integrating numerically the Gauss' equations for a, e and ω using as initial conditions different values of a 0 andē 0 , as shown in Fig. 2, and usingω 0 = 0. The numerical integration is performed in MATLAB using ode45 with absolute and relative tolerance equal to 10 −12 . The considered value of is 10 −4 m/s 2 . The numerical integration is interrupted when the targeted final semi-major axisā f is reached; the corresponding final time is ToF num . Results show that the relative error is higher for higher values of the eccentricity but remains lower than 2% for Geocentric transfers with eccentricities up to 0.7. For the four points highlighted with a red mark in Fig. 2, the maximum relative errors in the expression for the variation of the orbital elements can be computed, by comparing the results of the analytical equations with those from the numerical integration. The errors are reported in Table 1, along with the values of a 0 , a and e. It is important to stress that the control law presented in this section gives higher errors in the variation of the orbital elements, compared to the control laws presented in the rest of the paper. For conciseness, therefore, no further analysis of the error will be presented in the next (a) (b)

Fig. 3
Relative difference between numeric and analytical V for Heliocentric transfers-variation of semimajor axis with tangential thrust sections. Figure 3 shows the relative difference in V for Heliocentric transfers with initial semi-major axisā 0 ∈ [1, 1.5] AU and characterised by ā ∈ [0.1, 0.5] AU; two different values of the eccentricity have been considered. The error between analytical and numerical V is bigger for Heliocentric than for Geocentric transfers. The analytical equation for the V can, therefore, be used to approximate the cost of Geocentric transfers for small values of the low-thrust acceleration. However, the analytical equation gives non-negligible errors for some interplanetary transfers. Heliocentric transfers are, in fact, characterised by longer orbital periods; over these long orbital revolutions, the value of the semi-major axis cannot be assumed to remain constant. Moreover, the ratio between the low-thrust acceleration and the local gravity field is bigger for Heliocentric than for Geocentric transfers.

Variation of semi-major axis with radial and transverse thrust
In Burt (1967), the author proposed a thrust profile with radial and transverse thrust that changes the semi-major axis of non-circular orbits, without any variation in the eccentricity. In order to obtain this thrust profile, the expressions for the variation of the orbital elements with respect to the eccentric anomaly E are required. From the Gauss' equation for da/dt (Eq. 2), one can derive a relationship for da/d E: Similarly, it is possible to obtain the equations for the variation of the eccentricity with respect to the eccentric anomaly, and, the variation of the argument of the periapsis with respect to the eccentric anomaly: The variation of a with no secular variation of e can be obtained with a radial acceleration that changes sign when E = 0 and E = π (Burt 1967). Integration of the equation for de/d E (Eq. 18) from E = 0 to E = 2π, under the assumption that the osculating elements during one orbital revolution remain equal to the elements at E = 0, results in: From Eq. (20), imposing ē 2π = 0, the following can be obtained: Expressions for f R and f T can be obtained using Eq. (21) and considering that f 2 R + f 2 T = . The expressions for the angles α and β that satisfy Eq. (21), which are not given by Burt (1967), are By using the expression in Eq. (21) to derive f R , one can obtain the equation for the secular variations of a, as reported also by Burt (1967): Analogously, it is possible to verify that: The rest of this subsection reports additional equations that are not available in the literature and that are derived here for the first time. The integration of Eq. (23) gives: from which it is possible to obtain an analytical expression for the variation ofā with time : From Eq. (26), and considering thatā f =ā(ToF), it is then possible to state the following: Proposition 2 Given Eq. (26), the cost of a transfer to change the semi-major axis fromā 0 toā f , with dē/dt = 0, using the thrust profile defined in Eqs. (21) and (22), is: Equation 27 shows no singularity forē = 0. It is important to stress that from the integration of Eq. (19), one can see that the proposed solution to the low-thrust motion does not cause any long-term variation of the argument of periapsis during the transfer: dω/dt = 0. Figure 4 compares the V required for the variation ofā using the two solutions presented in this section. The value of is set to 10 −4 m/s 2 . In particular, Fig. 4 shows the difference in V when using Eqs. (16) and (27); this difference increases with increasing eccentricity but is limited to very small values. For e = 0, the two equations coincide.

Variation of the eccentricity or combination of eccentricity and other orbital elements
In this section, analytical solutions for the variation of the eccentricity, or combinations of orbital elements including the eccentricity, are presented. Other thrust profiles, different from the ones presented in this section, have been proposed in the literature (Petropoulos 2003); however, they do not lead to a closed form analytical solution and, therefore, are not reported hereafter.

Variation of the eccentricity with no variations of semi-major axis and argument of periapsis
In this section two different strategies are presented to obtain a variation of the eccentricity without long-term variation of semi-major axis and argument of periapsis. One is based on the work presented by Pollard (2000), and the other is based on the work presented by Burt (1967).

Thrust perpendicular to the major axis of the orbit
Pollard (2000) presented a strategy to change the eccentricity keeping the semi-major axis constant, using thrusting arcs centred at the periapsis and apoapsis of the orbit. For consistency with the rest of this paper, the equations presented here have been derived for thrust continuously applied over the entire orbit; they are, therefore, to be considered a novel derivation. We consider the following thrust profile, perpendicular to the major axis of the orbit: whereē 0 is the initial eccentricity andē f is the final eccentricity.
Substituting these values for α and β in Eqs. (17), (18) and (19) results in: The variations of semi-major axis and argument of periapsis over one orbital period are ā 2π = ω 2π = 0, while the secular variation of the eccentricity is The integration of Eq. 30 allows one to obtain the expression for the variation of the eccentricity with time:ē Considering Eq. (31), it is then possible to state the following, regarding the V required to realise the transfer: Proposition 3 Given Eq. (31), the cost of the transfer to change the eccentricity fromē 0 tō e f , withā f =ā 0 andω f =ω 0 , using the thrust profile defined in Eq. (28), is: Equation (32) presents no singularity forē 0 = 0 orē f = 0.

Transverse thrust
In the work presented by Burt (1967), it was proposed to change e, without secular variations of a and ω, by applying a transverse acceleration f T reversed in sign when E = ±π/2 (crossing of the minor axis). The corresponding thrust angles are: The variations of the orbital elements over one orbital period are It is possible to obtain an analytical expression for the variation ofē with time. The secular rate of change of the eccentricity is The integration of Eq. (35) results in: Considering Eq. (36), and thatē(ToF) =ē f , it is possible to obtain the following expression for the V : Proposition 4 Given Eq. (36), the cost of a transfer to change the eccentricity fromē 0 toē f , withā f =ā 0 andω f =ω 0 , using the thrust profile defined in Eq. (33), is: Equation (37) presents no singularity forē 0 = 0 orē f = 0.

Combined variation of eccentricity and inclination without variation of other orbital elements
According to Pollard (2000), the simultaneous variation of e and i can be obtained using the following thrust profile: In Pollard (2000), the equations for this type of transfer are derived for the case of thrust applied on arcs centred at the periapsis and apoapsis of the orbit. In this paper, instead, a continuous thrust applied over the entire orbit is considered. The equations for the secular variation of a, e and ω are derived from the results presented in Sect. 4.1.1, with the addition of the term cos β in the expressions for f R and f T , and considering the addition of the term due to f N in dω/dt: The variation of the inclination with the eccentric anomaly is It is now possible to derive the secular variation of i by integrating Eq. (40) from 0 to 2π with the sign of β flipping at the points along the orbit where the spacecraft is crossing the minor axis: Combining the equations for the secular variations of eccentricity and inclination provides: It is possible to integrate Eq. (42) fromī 0 toī f and fromē 0 toē f , under the assumption that ω and β are both constant. This is true when sinω = 0 (see Eq. 39 for dω/dt), that is when ω = 0 orω = π. The integration provides an expression for the value of β: The equation for the variation of the eccentricity with time is similar to Eq. (31), but for the term in β:ē From dī/dt one can derive an expression for the variation of the inclinationī with time, which is not given in the literature: where Note that an analytic equation forī(t), corresponding to Eq. (45), was not derived in Pollard (2000). Considering Eq. (44), it is then possible to state the following: Proposition 5 Given Eq. (44), the cost of a transfer to change the eccentricity fromē 0 toē f and to change the inclination fromī 0 toī f , using the thrust profile defined in Eqs. (38) and (43), is: Equation (47) presents no singularity forē 0 = 0 orē f = 0. It is important to stress that this analytical solution, and the corresponding transfer, are valid only when there is a nonzero secular change of the eccentricity, because β = 0 whenē 0 =ē f . Therefore, this thrust profile cannot be used to obtain a pure variation of the inclination. Another important point to stress is that, despite the out-of-plane component, this thrust profile causes no variation of the right ascension. This can be verified by looking at the expression for the variation of the right ascension with the eccentric anomaly, as presented in Pollard (2000): It is easy to verify that the integral of Eq. (48) over one orbital period, with β changing sign at the crossing of the minor axis, results in Ω 2π = 0.

Comparison of laws for the variation of the eccentricity only
Figures 5 and 6 show the V required for the variation ofē without variation of the semimajor axis, for different values of the semi-major axis (7000 km, 8000 km, 9000 km and 10,000 km) and using = 10 −4 m/s 2 . In particular, Fig. 5 shows the V relative to the thrust profile defined by Pollard (Eq. 32) and Fig. 6 shows the V relative to the thrust profile defined by Burt (Eq. 37). Figures 5 and 6 show that the V cost is higher when the thrust profile defined by Burt is used.

Variation of inclination
Following the derivations presented by Ruggiero et al. (2011), the maximum instantaneous variation of inclination can be obtained using the following out of plane thrust profile: Because of the term cos ν in the Gauss' equation for di/dt, in order to obtain a nonzero variation of the inclination, β has to change sign at ν = π/2 and ν = 3π/2. Equation (49) can, therefore, be rewritten as Because of the term sin ν in the equation for dΩ/dt, this thrust profile causes no variation ofΩ. There is, however, a variation ofω due to β. No simple expression is available for the variation ofω andī with time whenē 0 = 0. Forē 0 = 0, instead, it is possible to find analytical equations for the variation of the orbital elements and for the cost of the transfer. The variation with time of the inclination is The integration of Eq. (51) gives the expression for the variation of the inclination with time: Considering Eq. (52), and ī =ī f −ī 0 =ī(ToF) −ī 0 , it is possible to state that:

Proposition 6
The cost of the transfer to change the inclination by ī following the thrust profile defined in Eq. (50) is given by: 6 Variation of the right ascension of the ascending node Considering the Gauss' equation for dΩ/dt, the instantaneous maximum variation of Ω is obtained when This value of β does not cause any variation ofī. No analytical expression is available for the time variation ofΩ, when using the values of β given in Eq. (54), in the general case when e 0 = 0. Analytical equations can be derived whenē 0 = 0. In this case: The expression forΩ(t), obtained from integration of Eq. (55), is: To the author's best knowledge, Eq. (56) Equation (57) is valid only forī = 0. Forī = 0, Ω is not defined.

Variation of argument of the periapsis
Different solutions are available for the variation of the argument of periapsis. The ones that can provide analytical solutions are presented in this section. Other thrust profiles for the variation of ω that do not result in analytical solutions are presented by Petropoulos (2003) and Ruggiero et al. (2011). The results in this sections are valid only whenē = 0. When e = 0, the argument of perigee is not defined.

Variation of the argument of the periapsis without variation of semi-major axis and eccentricity
The next subsections present thrust profiles that cause changes in the argument of the periapsis while keeping the secular variations of semi-major axis and eccentricity constant. We start from the work reported by Burt (1967) and derive analytical formulas for the secular variation of ω and for the associated V .

Transverse thrust
It is possible to changeω, without variations ofā andē, by using a thrust with nonzero transverse component f T , reversed in sign at the crossing of the major axis (that is, at E = 0 and E = π): If f T is reversed in sign at the crossing of the major axis, Eqs. (17) and (18) give while, Eq. (19) gives: The integration of Eq. (60) gives a linear variation ofω with time: Then, from Eq. (61), one can derive an analytical expression for the V cost of the transfer, as stated in the following proposition: Proposition 7 Given Eq. (61), the cost required for a transfer to change the argument of periapsis fromω 0 toω f , withā f =ā 0 andē f =ē 0 , and following the thrust profile defined in Eq. (58), is

Radial thrust
Another possible thrust profile that changesω while giving no secular variation ofā andē is the unidirectional radial acceleration (Burt 1967): In this case, the secular variations of ω are: which gives the following expression for the time variation ofω: Considering Eq. (65), it is possible to state that:

Thrust parallel to the major axis of the ellipse
The last thrusting strategy presented in this section was originally proposed by Pollard (2000); it is based on an in-plane acceleration parallel to the major axis of the ellipse: This results in the following azimuth and elevation angles: The derivation by Pollard (2000) considers a thrust applied over two symmetric arcs centred at the periapsis and apoapsis of the orbit. In the following, in order to be consistent with the rest of the paper, the thrust is assumed to be applied over the entire orbit. By using Eqs. (17) and (18) it is possible to verify that, with the thrust profile defined in Eq. (68), the secular variations of semi-major axis and eccentricity are both zero. The secular variation of the argument of the periapsis is: Equation (69) coincides with the equation given by Pollard (2000) when the two thrust arcs at periapsis and apoapsis have amplitude equal to π. From Eq. (69) one can derive the equation for the time variation of ω:ω Considering Eq. (70), it is possible to state that:

Proposition 9
The cost of a transfer to go fromω 0 toω f , withā f =ā 0 andē f =ē 0 , using the thrust profile of Eq. (68), is: Figures 7, 8 and 9 show the V cost for the variation ofω, using the analytical solution defined in Sect. 7.1, for different values of ω,ē andā. The value of is set to 10 −4 m/s 2 . Figures 7 and 8 show the results for transverse and radial low-thrust acceleration, respectively, as proposed in Burt (1967). Figure 9 presents instead the results for the case of low-thrust acceleration parallel to the major axis of the ellipse. Results show that the law with acceleration parallel to the major axis of the ellipse is the most advantageous one in terms of V , while the law with radial acceleration is the most expensive.

Comparison of laws for the variation of argument of the periapsis
The analytical solutions presented in Sects. 3 to 7 are summarised in "Appendix A", where a taxonomy of the analytical solutions is also given. The analytical solutions for low-thrust transfers presented in the previous sections are obtained assuming that the low-thrust acceleration is the only acceleration acting on the spacecraft; no other orbital perturbation is considered. In this section, new analytical solutions are derived for the Two-Point Boundary Value Problem (TPBVP) associated to the low-thrust transfer between inclined circular orbits, with different values of the initial and final right ascension of the ascending node, under the additional effect of the second order zonal harmonic of the Earth's gravitational potential, J 2 . Analytical laws are also derived for the variations of the orbital elements with time and for the cost of the transfer. As in the previous sections, it is assumed that the transfer trajectory is composed of several revolutions; during each revolution, the value of the considered orbital elements can be assumed to be constant. The proposed low-thrust transfers aim at realising the simultaneous variation of semi-major axis, inclination and right ascension of circular orbits, in a given time of flight ToF. The resulting TPBVP can be expressed as follows: The state x is: and the initial and final conditions are expressed as The right-hand side of the first equation in Eq. (72) is: The expressions for dā/dt, dī/dt and dΩ/dt can be derived from the following set of Gauss' planetary equations for e = 0 (see Section 10.2.3 at page 262 in Mengali and Quarta 2006 for the terms in J 2 ): where R ⊕ is the Earth's radius. The control vectorũ in Eq. (72) defines the thrust orientation and the duration of each thrust phase, as specified more in details in the next subsections. In the following we will devise three strategies that exploit the secular drift in Ω caused by J 2 in combination with the low-thrust action, to obtain the desired variation ofā,ī andΩ. The three strategies will be applied to the case in whichā 0 <ā f ; however, the same procedure can be used to define analogous strategies forā 0 >ā f .

Strategy 1: ¯J 2 + ( ā, ī)
For the first strategy, the transfer from the initial to the final orbit is realised in two phases. During the first phase, denoted as Ω J 2 and characterised by a time of flight ToF 1 , J 2 is exploited to obtain a given variation of the right ascension. In the second phase, denoted as ( ā, ī ) and characterised by a time of flight ToF 2 , a low-thrust action is applied to change semi-major axis and inclination from their initial to their final values. In this phase, the low-thrust action is applied along two thrust arcs of constant angular span 2ψ, with constant elevation angleβ. The variation of Ω, due to J 2 only, during the second phase is tuned in such a way that, starting from the value of Ω reached at the end of phase one, the final targeted right ascension can be achieved. The vector of controls for Strategy 1 isũ = [ψ,β, ToF 1 , ToF 2 ] T . The symbolsψ andβ are used in the rest of this section to denote constant control angles throughout the transfer. The two phases of Strategy 1 are schematically represented in Fig. 10.
The following subsections describe the two phases in more detail and derive the system of four equations required to compute the four components ofũ that solve the TPBVP (Problem Fig. 10 Representation of the two phases of Strategy 1. LT stands for low-thrust (72)). Analytical equations for the variation of the orbital elements during the transfer and for the V cost of the transfer are also provided.

First phase
During the first phase, the low-thrust engine is assumed to be off. The semi-major axis and inclination remain equal to their initial values,ā 0 andī 0 . This means that, considering the range of semi-major axes covered during the transfer [ā 0 ,ā f ] and the assumption that a 0 <ā f , the drift ofΩ due to J 2 is maximum and the resulting secular variation of Ω with time is:Ω where the time of flight associated to this phase is ToF 1 andā 1 (t) =ā 0 ,ī 1 (t) =ī 0 ∀t ∈ [0, ToF 1 ]. The right ascension of the ascending node at the end of the first phase isΩ 1 f = Ω 1 (ToF 1 ). ToF 1 and Ω 1 f satisfy the following relationship:

Second phase
During the second phase, the low-thrust engine is active during two thrust arcs per orbital revolution. During this phase, the semi-major axis and inclination change from their initial to their final values. Moreover, the drift of right ascension due to J 2 is such that, at the end of the transfer, the right ascension will have changed fromΩ 1 f toΩ f . Note that the variation of right ascension takes place only because of J 2 and not because of the low-thrust acceleration (as explained in more details later, see Eq. 85). The time of flight of the second phase is ToF 2 . The simultaneous variation ofā andī in a given time of flight ToF 2 can be obtained with two tangential thrust arcs of angular span 2ψ (constant during the transfer), constant azimuth angleᾱ = 0 and constant elevation angleβ, centred at the nodal points of the orbit, as shown in Fig. 11. The elevation angleβ has equal and opposite values along the two thrust arcs. The relevant equations for this transfer can be obtained from Eq. (76), using α = 0: a 4 3 sin i cos i sin ν) .
(79) Fig. 11 Strategy 1, phase 2: thrust arcs centred at the nodal points of the orbit The resulting secular variations of a and i are computed as follows: In Eq. (80) it is assumed that the variations of a and i in one revolution are negligible, so that they can be considered constant during the integration over ν. In order forā andī to reach their final values simultaneously, the following equation is integrated fromā 0 toā f and from i 0 toī f : From the integration of Eq. (81), one can derive the elevation angleβ required to realise the transfer fromā 0 toā f and fromī 0 toī f , as a function ofψ,ā 0 ,ā f ,ī 0 andī f : Finally, the expressions for the variation ofā andī with time, during the second phase ( ā, ī ), are (see Eqs. 80 and 81): with t ∈ [ToF 1 , ToF 1 + ToF 2 ]. From Eqs. (82) and (83), the time of flight of the second phase of the transfer can be expressed as a function ofψ and of the initial and final orbital elements: During the second phase of the transfer, the right ascension changes because of J 2 and because of the variation of the semi-major axis and inclination with time. The rate of Ω with time is given in the last equation in Eq. (76). The first term on the right-hand side of this equation cancels out when integrating in ν over one revolution, becauseβ takes opposite values on the two thrust arcs. The secular variation of the right ascension is, therefore, due only to J 2 : It is possible to write the expression of the variation ofΩ withī as Equation (86) can be integrated using the following expression for the variation of the semimajor axis with time (obtained from Eq. 83): Substituting Eqs. (87) in (86) and integrating fromΩ 1 f toΩ(ToF 1 + ToF 2 ) and fromī 0 tō i f results inΩ where In addition, it is required thatΩ(ToF 1 + ToF 2 ) =Ω f . Equation (89) becomes singular when ī f −ī 0 is small. In the caseī f =ī 0 , an alternative, non-singular formulation is available.
Whenī f =ī 0 ,β = 0 and the expression forΩ 2 f can be derived by integrating the expression of dΩ/dā which is obtained from Eqs. (80) and (85): The cost of the transfer can be computed analytically considering that the engine is active during a fraction of the orbital period equal to therefore, The final expression for the cost of the transfer is Equation (93) expresses the V required to change semi-major and inclination under the effect of J 2 . As such, it can be considered an extension to Eqs. (16), (27) and (53).

Solution method
Given the total time of flight (ToF), a complete solution for the transfer can be computed by solving the following system with respect to the vectorũ, that isβ,ψ, ToF 1 and ToF 2 . In particular, substituting the first of Eqs. (94) into the second provides a nonlinear equation inψ that can be solved using an algorithm based on a combination of bisection, secant and inverse quadratic interpolation methods, implemented in the function fzero in MATLAB. Onceψ is found, the rest of Eqs.
(94) allows one to findβ, ToF 1 and ToF 2 . Figure 12 shows an example of transfer realised using Strategy 1, with time of flight of 600 days and low-thrust acceleration equal to 1.5 10 −4 m/s 2 . The variation of the semimajor axis is from 10,000 to 24,200 km, the inclination changes from 51 to 56 • , and the right ascension changes from 0 to 150 • . The relevant parameters of the solution for the lowthrust control areψ = 67.02 • andβ = 14.1 • . The times of flight of the two phases are ToF 1 = 359.1 days and ToF 2 = 240.9 days. The cost of the transfer is V = 2.32 km/s. Figure 12 shows that, during the first phase, the semi-major axis and inclination stay constant, while the right ascension changes because of J 2 . In the second phase, the semi-major axis and inclination change from their initial to their final values and the variation of the right ascension is such as to obtain, at the end of the transfer, and by means of J 2 only, the final desired value.
The evolution of the orbital elements shown in Fig. 12 is obtained using the expressions in Eqs. (77)  the numerical integration of the Gauss' equations, using MATLAB ode45 with absolute and relative tolerance equal to 10 −10 . Results show that the absolute errors of the orbital elements were limited to the following values: semi-major axis error < 31 km; inclination error < 0.03 • ; RAAN error < 1.35 • .

Strategy 2: ( ā, ī) + ¯Ǎ
lso in the second strategy the transfer is realised in two phases. The first phase, denoted as ( ā, ī ), is similar to the second phase of Strategy 1. The low-thrust action is applied over two thrust arcs of constant angular span 2ψ 1 , with constant elevation angleβ 1 , for a time of flight equal to ToF 1 . In the second phase, of time length ToF 2 and denoted as Ω β , an out-of-plane control thrust is used to changeΩ to its final value. During the second phase, the low-thrust action is applied over two thrust arcs of constant angular span 2ψ 2 . The vectorũ associated to the TPBVP is nowũ = ψ 1 ,β 1 ,ψ 2 , ToF 1 , ToF 2 T . The two phases of Strategy 2 are schematically represented in Fig. 13 The following subsections describe the two phases in more detail. Analytical equations for the variation of the orbital elements during the transfer and for the V of the transfer are also provided.

First phase
During the first phase,ā andī are changed from their initial values,ā 0 andī 0 , to their final values,ā f andī f , as described in Sect. 8.1. The engine is operated on two arcs per orbital revolution of amplitude 2ψ 1 with elevation angleβ 1 . The time of flight of the first phase is denoted as ToF 1 . The right ascension at the end of the first phase is computed from Eq. (88) as:Ω The cost, V 1 , is computed using Eq. (93).

Second phase
During the second phase,Ω is changed by applying an out-of-plane thrust action; two thrust arcs of angular span 2ψ 2 and elevation angle β 2 = 90 • are centred at ν = 90 and ν = 270 • , as shown in Fig. 14. The sign ofβ is opposite on the two thrust arcs.
Using the Gauss' equation for Ω, Eqs. (2) and (3), it is possible to obtain the variation of Ω with the argument of latitude ν, due to the low-thrust action: four (Eqs. (82), (84), (100) and ToF 1 + ToF 2 = ToF): It is possible, however, to define different arbitrary values of ToF 1 < ToF and compute the corresponding values ofψ 1 ,β 1 ,ψ 2 and ToF 2 . In particular, it is possible to find a ToF 1 such that the V of the transfer is the minimum possible value. An example of transfer realised with this strategy is shown in Fig. 15. The initial and final orbital elements are those used for the example reported in Sect. 8.1. The minimum V transfer is realised with V = 2.31 km/s and the parameters of the low-thrust control areψ 1 = 33.75 • ,β 1 = 11.82 • andψ 2 = 0.29 • . The times of flight of the two phases are ToF 1 = 474.07 days and ToF 2 = 125.93. The variation of the orbital elements is reported in Fig. 15. During the first phase, the semi-major axis and inclination change from their initial to their final values while the variation of the right ascension is due only to J 2 . Note how the drift ofΩ is slower than in the example reported in Fig. 12. This is because the semi-major axis increases considerably during the first phase, thus reducing dΩ/dt. In the second phase the semi-major axis and inclination remain constant while the right ascension changes because of the out-of-plane thrust and because of J 2 . The minimum V transfer corresponds to a small value ofψ 2 , meaning that the out-of-plane thrust, more expensive than the in-plane thrust, is applied on very short arcs. The transfer shown in Fig. 15 is the one characterised by the lowest value of V , among all the possible solutions. The V s obtained for different values of ToF 1 are reported in Fig. 16. Figure 17 shows the variation of the right ascension during the transfer for two possible solutions: ToF 1 = 474.07 days and ToF 1 = 505 days. The higher V obtained when ToF 1 = 505 days is due to the greater variation ofΩ obtained with out-of-plane thrust during the second phase.ψ 2 is indeed equal to 61.66 • when ToF 1 = 505 days. For Strategy 2, the comparison with the numerical integration of the Gauss' equations shows that the absolute errors of the orbital elements during the transfer are limited to the following values: semi-major axis error < 8 km; inclination error < 3.4 10 −4 deg; RAAN error < 0.35 • .

Strategy 3: ( ā, ¯J
2 ) + ī Strategy 3 is composed of two phases: the first one changes semi-major axis and right ascension and is denoted as ( ā, Ω J 2 ), while the second one uses the low-thrust acceleration to change the inclination, and is denoted as ī . During the second phase,Ω changes because of J 2 . The first phase has time span ToF 1 and makes use of a low-thrust action on two thrust arcs of constant semi-amplitudeψ 1 , while the second phase has time span ToF 2 and uses arcs of constant semi-amplitudeψ 2 per orbit.  . 18 Representation of the two phases of Strategy 3. LT stands for low-thrust isũ = ψ 1 ,ψ 2 , ToF 1 , ToF 2 T . The strategy is schematically shown in Fig. 18. More details are given hereafter.

First phase
The transfer in the first phase of Strategy 3 is analogous to the transfer during the second phase of Strategy 1, withβ 1 = 0. During the first phase, a tangential thrust is used to increase the semi-major axis fromā 0 toā f . The low-thrust engine is operated during two thrust arcs per revolution, of semi-amplitudeψ 1 . Considering Eq. (84) withβ = 0, the time of flight for the first phase of the transfer is given by The cost of this phase can be computed analytically using Eq. (93) withβ = 0: Equation (104) extends Eqs. (16), (27) and (57). The right ascension at the end of the transfer is derived by integrating the expression of dΩ/dā which is obtained from the equations for dΩ/dt and dā/dt withī(t) =ī 0 (refer to Sect. 2). This results in:

Second phase
The transfer in the second phase of Strategy 3 is analogous to the transfer during the second phase of Strategy 1, withβ = 90 • . The thrust is applied during two thrust arcs per revolution, of semi-amplitudeψ 2 , with elevation angleβ 2 = 90 • . This causes the inclination to change fromī 0 toī f , while the semi-major axis stays constant atā f . The time required can be computed integrating the equation for dī/dt (Eq. 80) fromī 0 toī f , usingā =ā f : The cost of this phase is derived from Eqs. (106) and (92):

Solution method
A complete transfer using Strategy 3 can be defined by solving the following system of equations, for a given time of flight ToF, withψ 1 ,ψ 2 , ToF 1 and ToF 2 as unknowns: An example of transfer realised with this strategy is shown in Fig. 19 Figure 20 shows the variation of the orbital elements during the example of transfer considered in the previous sections, for the three proposed strategies. Figure 21 shows the V required to realise the transfer defined in Sect. 8.1.3, for different values of the times of flight and using the three proposed strategies. For Strategy 2, for which a unique solution does not exist, the solutions plotted are the ones corresponding to the values of ToF 1 and ToF 2 providing the lowest value of the V for that transfer. The fact that, for Strategy 2, ToF 1 and ToF 2 are chosen to minimise the V , explains the noisy behaviour of the curve relative to Strategy 2 in Fig. 21.

Comparison of the three strategies
Results show that the first and second strategies are those giving the lowest V s.

Conclusions
This paper has presented a number of analytical equations for the quick, approximated calculation of the V cost of low-thrust orbital transfers. The same equations provide an approximation to the variation of the mean orbital elements with time.
In particular, the first part of the paper collects a number of analytical solutions, available in the literature, for the variation of the orbital elements with a low-thrust acceleration, in the absence of orbital perturbations. The formulae in this collection have been extended by adding new analytical expressions for the total transfer V and the variations of the orbital elements. The most relevant novel equations have been presented as propositions. Two tables have been presented in Appendix to summarise the available analytical solutions. Results, on the range of case studies investigated in this paper, have shown that for most of the analytical solutions the maximum relative error in the V estimate is 0.04 %. Only in one case, the maximum error is higher than 0.04 %. The low relative error and the fact that the computation of the V simply requires the evaluation of an inexpensive analytical formula, allows one to rapidly analyse a large number of cases with good accuracy. The use of these formulae is limited by the working assumptions on the control laws and the use of standard Keplerian elements. Furthermore, the use of mean elements implies that these formulae are more appropriate to study long spirals than rendezvous where a precise calculation of the osculating elements is required.
By capitalising on the formulae and methodology presented in the first part of the paper, the second part introduced three strategies to obtain the simultaneous variation of semimajor axis, inclination and right ascension of the ascending nodes of circular orbits, by combining the effect of a low-thrust acceleration with the one of the second zonal harmonic of the gravitational potential, J 2 . Also in this case the validation of the formulae by direct numerical integration of Gauss' equations has shown that the error remains contained in all the cases analysed in this paper. Thus, the three strategies, and associated analytical expressions, provide a quick and accurate estimation of the V when a low-thrust spiral and natural dynamics can be combined to achieve the desired transfer.
In conclusion, the set of equations presented throughout this paper represent a fast and sufficiently accurate tool to obtain a preliminary estimation of the V cost of low-thrust spiral transfers.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Table 2 summarises the analytical solutions presented in this paper for the orbit transfer from (ā 0 ,ē 0 ,ī 0 ,Ω 0 ,ω 0 ) to (ā f ,ē f ,ī f ,Ω f ,ω f ), without orbital perturbations. For each solution, the number of the equations which give the thrust angles, the time evolution of the orbital elements and the V are given. When an orbital element does not change during the transfer this is directly reported in the table giving its initial value. The last column gives the number of the section where the solution is described in detail. N.A. stands for either not available or not applicable. Table 3 shows, for each analytical solution, if it was obtained computing

t)ē(t)ī(t)Ω(t)ω(t) V Section
Semi-major axis The entries correspond to the equation number for the expressions of α and β, for the time variation of the orbital elements (when applicable) and for the V equation. The last column reports the Section number of the paper Table 3 Taxonomy of analytical low-thrust solutions. For solutions 6 and 7,ē = 0 and thereforeω is not defined the maximum instantaneous rate of change of the orbital elements or using other types of guidance laws. It shows, moreover, which orbital elements are to be changed and which ones stay constant during the transfer. The last column reports if new analytical equations were derived in this paper.