Nonconservative extension of Keplerian integrals and a new class of integrable system

Finding first integrals is fundamental for characterizing a dynamical system. The motion is confined to submanifolds of lower dimensions on which the orbits evolve, providing an intuitive interpretation of the dynamics and reducing the complexity of the system. In addition, conserved quantities are good candidates when applying the second method of Lyapunov for stability analysis. Conservative systems related to central forces are typical examples of (Liouville) integrability, and provide useful analytic results. Hamiltonian systems have been widely analyzed in the classical and modern literature to determine adequate integrability conditions. The existence of first integrals under the action of small perturbations occupied Poincaré (1892, Chap. V) back in the 19th century. Later, Emmy Noether (1918) established in her celebrated theorem that conservation laws can be understood as the system exhibiting dynamical symmetries. In a more general framework, Yoshida (1983a, b) analyzed the conditions that yield algebraic first integrals of generic systems. He relied on the Kowalevski exponents for characterizing the singularities of the solutions and derived the necessary conditions for existence of first integrals exploiting similarity transformations.

Conservation laws are sensitive to perturbations and their generalization is not straightforward. For example, the Jacobi integral no longer holds when transforming the circular restricted three-body problem to the elliptic case (Xia, 1993). Nevertheless, Contopoulos (1967) was able to find approximate conservation laws for orbits of small eccentricities. Szebehely and Giacaglia (1964) benefited from the similarities between the elliptic and the circular problems in order to define transformations connecting them. Hénon and Heiles (1964) deepened in the nature of conservation laws and reviewed the concepts of isolating and nonisolating integrals. Their study introduced a similarity transformation that embeds one of the constants of motion and transforms the original problem into a simplified one, reducing the degrees of freedom (Arnold et al., 2007, §3.2). Carpintero (2008) proposed a numerical method for finding the dimension of the manifold in which orbits evolve, i.e. the number of isolating integrals that the system admits.

The conditions for existence of integrals of motion under nonconservative perturbations received important attention in the past due to their profound implications. Djukic and Vujanovic (1975) advanced on Noether’s theorem and included nonconservative forces in the derivation. Relying on Hamilton’s variational principle, they not only extended Noether’s theorem, but also its inverse form and the Noether-Bessel-Hagen and Killing equations. Later studies by Djukic and Sutela (1984) sought integrating factors that yield conservation laws upon integration. Examples of application of Noether’s theorem to constrained nonconservative systems can be found in the work of Bahar and Kwatny (1987). Honein et al. (1991) arrived to a compact formulation using what was later called the neutral action method. Remarkable applications exploiting Noether’s symmetries span from cosmology (Capozziello et al., 2009; Basilakos et al., 2011) to string theory (Beisert et al., 2008), field theory (Halpern et al., 1977), and fluid models (Narayan et al., 1987). In the book by Olver (2000, Chaps. 4 and 5), an exhaustive review of the connection between symmetries and conservation laws is provided within the framework of Lie algebras. We refer to Arnold et al. (2007, Chap. 3) for a formal derivation of Noether’s theorem, and a discussion on the connection between conservation laws and dynamical symmetries.

Integrals of motion are often useful for finding analytic or semi-analytic solutions to a given problem. The acclaimed solution to the satellite main problem by Brouwer (1959) is a clear example of the decisive role of conserved quantities in deriving solutions in closed form. By perturbing the Delaunay elements, Brouwer and Hori (1961) solved the dynamics of a satellite subject to atmospheric drag and the oblateness of the primary. They proved the usefulness of canonical transformations even in the context of nonconservative problems. Whittaker (1917, pp. 81–82) approached the problem of a central force depending on powers of the radial distance, $r^{n}$, and found that there are only fourteen values of $n$ for which the problem can be integrated in closed form using elementary functions or elliptic integrals. Later, he discussed the solvability conditions for equations involving square roots of polynomials (Whittaker and Watson, 1927, p. 512). Broucke (1980) advanced on Whittaker’s results and found six potentials that are a generalization of the integrable central forces discussed by the latter. These potentials include the referred fourteen values of $n$ as particular cases. Numerical techniques for shaping the potential given the orbit solution were published by Carpintero and Aguilar (1998). Classical studies on the integrability of systems governed by central forces are based strongly on Newton’s theorem of revolving orbits.¹¹1Section IX, Book I, of Newton’s Principia is devoted to the motion of bodies in moveable orbits (De Motu Corporum in Orbibus mobilibus, deq; motu Apsidum, in the original latin version). In particular, Thm. XIV states that “The difference of the forces, by which two bodies may be made to move equally, one in a quiescent, the other in the same orbit revolving, is in a triplicate ratio of their common altitudes inversely”. Newton proved this theorem relying on elegant geometric constructions. The motivation behind this result was the development of a theory for explaining the precession of the orbit of the Moon. A detailed discussion about this theorem can be found in the book by Chandrasekhar (1995, pp. 184–201) The problem of the orbital precession caused by central forces was recently recovered by Adkins and McDonnell (2007), who considered potentials involving both powers and logarithms of the radial distance, and the special case of the Yukawa potential (Yukawa, 1935). Chashchina and Silagadze (2008) relied on Hamilton’s vector to simplify the analytic solutions found by Adkins and McDonnell (2007). More elaborated potentials have been explored for modeling the perihelion precession (Schmidt, 2008). The dynamics of a particle in Schwarzschild space-time can also be regarded as orbital motion perturbed by an effective potential depending on inverse powers of the radial distance (Chandrasekhar, 1983, p. 102).

Potentials depending linearly on the radial distance appear recursively in the literature because they render constant radial accelerations, relevant for the design of spacecraft trajectories propelled by continuous-thrust systems. The pioneering work by Tsien (1953) provided the explicit solution to the problem in terms of elliptic integrals, as predicted by Whittaker (1917, p. 81). By means of a special change of variables, Izzo and Biscani (2015) arrived to an elegant solution in terms of the Weierstrass elliptic functions. These functions were also exploited by MacMillan (1908) when he solved the dynamics of a particle attracted by a central force decreasing with $r^{-5}$. Urrutxua et al. (2015a) solved the Tsien problem using the Dromo formulation, which models orbital motion with a regular set of elements (Peláez et al., 2007; Urrutxua et al., 2015b; Roa et al., 2015). Advances on Dromo can be found in the works by Baù et al. (2015) and Roa and Peláez (2015). The case of a constant radial force was approached by Akella and Broucke (2002) from an energy-driven perspective. They studied in detail the roots of the polynomial appearing in the denominator of the equation to integrate, and connected their nature with the form of the solution. General considerations on the integrability of the problem can be found in the work of San-Juan et al. (2012).

Another relevant example of an integrable system in celestial mechanics is the Stark problem, governed by a constant acceleration fixed in the inertial frame. Lantoine and Russell (2011) provided the complete solution to the motion relying extensively on elliptic integrals and Jacobi elliptic functions. A compact form of the solution involving the Weierstrass elliptic functions was later presented by Biscani and Izzo (2014), who also exploited this formalism for building a secular theory for the post-Newtonian model (Biscani and Carloni, 2012). The Stark problem provides a simplified model of radiation pressure. In the more general case, the dynamics subject to this perturbation cannot be solved in closed form. An intuitive simplification that makes the problem integrable consists in assuming that the force due to the solar radiation pressure follows the direction of the Sun vector. The dynamics are equivalent to those governed by a Keplerian potential with a modified gravitational parameter.

The present paper introduces a new class of integrable system, governed by a biparametric nonconservative perturbation. This acceleration unifies various force models, including special cases of solar radiation pressure, low-thrust propulsion, and some particular configurations in general relativity. The problem is formulated in Sec. 2, where the biparametric acceleration is defined and then reduced to a uniparametric forcing thanks to a similarity transformation. The conservation laws for the energy and angular momentum are generalized to the nonconservative case by exploiting known symmetries of Kepler’s problem. Before solving the dynamics explicitly, we will prove that there are four cases that can be solved in closed form using elementary or elliptic functions. Sections 3–6 present the properties of each family of orbits and the corresponding trajectories are derived analytically. Section 7 is a summary of the solutions, which are unified in Sec. 8 introducing the Weierstrass elliptic functions. Finally, Sec. 9 discusses the connection with known solutions to similar problems, and with Schwarzschild geodesics.

The motion of a body orbiting a central mass under the action of an arbitrary perturbation $\textbf{{a}}_{p}$ is governed by

where $\mu$ denotes the gravitational parameter of the attracting body, r is the radiusvector of the particle, and $r=||\textbf{{r}}||$. In the case $\textbf{{a}}_{p}=0$, Eq. (2.1) reduces to Kepler’s problem, which can be solved in closed form.

It is well known that under the action of a radial perturbation of the form

in which $\xi$ is an arbitrary constant, the system (2.1) remains integrable. Perturbations of this type model various physical phenomena, like the solar radiation pressure acting on a surface perpendicular to the Sun vector, or simple control laws for spacecraft with solar electric propulsion. Indeed, the perturbed problem can be written

which is equivalent to Kepler’s problem with a modified gravitational parameter, $\mu^{\ast}=\mu(1-\xi)$.

The gravitational acceleration written in the intrinsic frame $\mathfrak{I}=\{\textbf{{t}},\textbf{{n}},\textbf{{b}}\}$, with

defined in terms of the velocity and angular momentum vectors, v and h respectively, reads

The flight-direction angle $\psi$ is given by

Consequently, perturbations of the form $\textbf{{a}}_{p}=\xi\mu\,(\textbf{{r}}/r^{3})$ result in

Since the component normal to the velocity does not perform any work, Bacon (1959) neglected its contribution and found that

furnishes another integrable problem: the trajectory of the particle is a logarithmic spiral.

Motivated by this line of thought, we can treat the normal component of the acceleration in Eq. (2.7) separately by introducing a second parameter, $\eta$:

The dimensionless parameters $\xi$ and $\eta$ are assumed constant. The forcing parameter $\xi$ controls the power exerted by the perturbation,

with $\mathcal{E}_{k}$ the Keplerian energy of the system. The second parameter $\eta$ scales the terms that do not perform any work. Introducing the accelerations (2.5) and (2.9) into Eq. (2.1) yields

having replaced $\eta$ by the new parameter

Since the motion is planar we shall introduce a system of polar coordinates $(r,\theta)$ to formulate the problem. The radial velocity is defined as

whereas the projection of the velocity normal to the radiusvector takes the form

Consequently, these geometric relations unveil the dynamical equations:

The dynamics are referred to an inertial frame $\mathfrak{I}=\{\textbf{{i}}_{\mathfrak{I}},\textbf{{j}}_{\mathfrak{I}},\textbf{{k}}_{\mathfrak{I}}\}$, with $\textbf{{k}}_{\mathfrak{I}}\parallel\textbf{{h}}$ and using the $x_{\mathfrak{I}}$-axis as the origin of angles. Figure 1 depicts the configuration of the problem. We restrict the study to prograde motions ($\dot{\theta}>0$) without losing generality.

For $\gamma=1$ the integrals of motion (2.33) and (2.36) reduce to the normalized equations of the energy and angular momentum, respectively. The condition on the radius given by Eq. (2.40) becomes

For the case $\tilde{\kappa}_{1}<0$ (elliptic solution) this translates into ${\rho}\in[{\rho}_{\mathrm{min}},{\rho}_{\mathrm{max}}]$, where

These limits are none other than the periapsis and apoapsis radii, provided that $\tilde{\kappa}_{1}$ and $\tilde{\kappa}_{2}$ relate to the semimajor axis and eccentricity by means of:

For $\tilde{\kappa}_{1}=0$ (parabolic case) the semimajor axis becomes infinite, and Eq. (3.1) has only one root corresponding to ${\rho}=\tilde{\kappa}_{2}^{2}/2$. Note that it must be $\tilde{\kappa}_{2}^{2}>-1/\tilde{\kappa}_{1}$. Similarly, for hyperbolic orbits ($\tilde{\kappa}_{1}>0$) it is ${\rho}\geq{\rho}_{\mathrm{min}}$ as ${\rho}_{\mathrm{max}}$ becomes negative.

Thus, the solution is simply a conic section:

The angle $\theta_{m}=\varOmega+\omega$ defines the direction of the line of apses in frame $\mathfrak{I}$, meaning that $\theta-\theta_{m}$ is the true anomaly. If $\tilde{\kappa}_{1}<0$, then ${\rho}(\theta_{m})={\rho}_{\mathrm{min}}$, and ${\rho}(\theta_{m}+\piup)={\rho}_{\mathrm{max}}$. The velocity $\tilde{v}$ follows from the integral of the energy:

It is minimum at apoapsis and maximum at periapsis.

Applying the similarity transformation $\mathscr{S}^{-1}$ to the previous solution leads to the extended integral

The factor $\mu^{\ast}=\mu(1-\xi)$ behaves as a modified gravitational parameter. This kind of solutions arise from, for example, the effect of the solar radiation pressure directed along the Sun-line on a particle following a heliocentric orbit (McInnes, 2004, p. 121).

The case $\gamma=2$ yields the family of generalized logarithmic spirals found by Roa et al. (2016a) in the context of interplanetary mission design. An extended version and the solution to the spiral Lambert problem can be found in following sequels (Roa and Peláez, 2016; Roa et al., 2016b). Negative values of $\tilde{\kappa}_{1}$ define the so called elliptic spirals, positive values yield hyperbolic spirals, and the limit case $\tilde{\kappa}_{1}=0$ corresponds to parabolic spirals, which turn out to be pure logarithmic spirals.

The condition in Eq. (2.40), yields the polynomial inequality:

The discriminant $\Delta$ of the polynomial $P_{\mathrm{nat}}(\rho)$ predicts the nature of the roots. It is

The intermediate value theorem shows that there is at least one real root. For the elliptic case $(\tilde{\kappa}_{1}<0)$ it is $\Delta<0$, meaning that the other two roots are complex conjugates. In the hyperbolic case $(\tilde{\kappa}_{1}>0)$ the sign of the discriminant depends on the values of $\tilde{\kappa}_{2}$: if $\tilde{\kappa}_{2}^{2}>3\tilde{\kappa}_{1}$ it is $\Delta>0$, and for $\tilde{\kappa}_{2}^{2}<3\tilde{\kappa}_{1}$ the discriminant is negative. This behavior yields two types of hyperbolic solutions.