# Classical and Relativistic Orbital Motions around a Mass-Varying Body

Authors: L. Iorio
Journal: SRX Physics
 (2010)

---

## Abstract

I work out the Newtonian and general relativistic effects due to an
isotropic mass lossM˙/M of a body on the orbital motion of a test particle around it; the present analysis is also valid for a variation G˙/G of the Newtonian constant of gravitation. Concerning the Newtonian
case, I use the Gauss equations for the variation of the elements and
obtain negative secular rates for the osculating semimajor axis a, the
eccentricity e, and the mean anomaly ℳ, while the argument of pericenter ω does not experience secular precession; the longitude of the
ascending node Ω and the inclination i remain unchanged as well. The true orbit, instead, expands, as shown by a
numerical integration of the equations of motion with MATHEMATICA; in fact, this is in agreement with the seemingly counter-intuitive decreasing of a and e because they refer to the osculating Keplerian
ellipses which approximate the trajectory at each instant. A comparison with the results obtained with different approaches by other researchers is made. General relativity induces positive secular rates
of the semimajor axis and the eccentricity completely negligible in the
present and future evolution of the solar system.


---

## Body

## 1. Introduction

In this paper I investigate the orbital effects induced by an isotropic variationṀ/M of the mass of a central body on the motion of a test particle; my analysis is valid also for a change Ġ/G of the Newtonian constant of gravitation. This problem, although interesting in itself, is not only an academic one because of the relevance that it may have on the ultimate destiny of planetary companions in many stellar systems in which the host star experiences a mass loss, like our Sun [1]. With respect to this aspect, my analysis may be helpful in driving future researches towards the implementation of long-term N-body simulations including the temporal change of GM as well, especially over timescales covering paleoclimate changes, up to the Red Giant Branch (RGB) phase in which some of the inner planets should be engulfed by the expanding Sun. Another problem, linked to the one investigated here, which has recently received attention is the observationally determined secular variation of the Astronomical Unit [2–5]. Moreover, increasing accuracy in astrometry pointing towards microarcsecond level [6] and long-term stability in clocks [7] require to consider the possibility that smaller and subtler perturbations will become soon detectable in the solar system. Also future planetary ephemerides should take into account Ṁ/M. Other phenomena which may, in principle, show connections with the problem treated here are the secular decrease of the semimajor axes of the LAGEOS satellites, amounting to 1.1 mm d-1, [8] and the increase of the lunar orbit's eccentricity of 0.9×10-11 yr-1 [9].Many treatments of the mass loss-driven orbital dynamics in the framework of the Newtonian mechanics, based on different approaches and laws of variation of the central body's mass, can be found in literature; see, for example, [2, 4, 10–18] and references therein. However, they are sometimes rather confused and involved, giving unclear results concerning the behavior of the Keplerian orbital elements and the true orbit.The plan of the paper is as follows. Section2 is devoted to a theoretical description of the phenomenon in a two-body scenario. By working in the Newtonian framework, I will analytically work out the changes after one orbital revolution experienced by all the Keplerian orbital elements of a test particle moving in the gravitational field of a central mass experiencing a variation of its GM linear in time. Then, I will clarify the meaning of the results obtained by performing a numerical integration of the equations of motion in order to visualize the true trajectory followed by the planet. Concerning the method adopted, I will use the Gauss perturbation equations [19, 20], which are valid for generic disturbing accelerations depending on position, velocity, and time, the  “standard” Keplerian orbital elements (the Type I according to, e.g.,  [16]) with the eccentric anomaly E as  “fast” angular variable. Other approaches and angular variables like, for example, the Lagrange perturbation equations [19, 20], the Type II orbital elements [16], and the mean anomaly ℳ could be used, but, in my opinion, at a price of major conceptual and computational difficulties. (Think, e.g., about the cumbersome expansions in terms of the mean anomaly and the Hansen coefficients, the subtleties concerning the choice of the independent variable in the Lagrange equations for the semimajor axis and the eccentricity [19].) With respect to possible connections with realistic situations, it should be noted that, after all, the Type I orbital elements are usually determined or improved in standard data reduction analyses of the motion of planets and (natural and artificial) satellites. My approach should, hopefully, appear more transparent and easy to interpret, although, at first sight, some counter-intuitive results concerning the semimajor axis and the eccentricity will be obtained; moreover, for the chosen time variation of the mass of the primary, no approximations are used in the calculations which are quite straightforward. However, it is important to stress that such allegedly puzzling features are only seemingly paradoxical because they will turn out to be in agreement with numerical integrations of the equations of motion, as explicitly shown by the numerous pictures depicted. Anyway, the interested reader is advised to look also at  [16] for a different approach. In Section 3 I will work within the general relativistic gravitoelectromagnetic framework by calculating the gravitoelectric effects on all the Keplerian orbital elements of a freely falling test particle in a nonstationary gravitational field. Section 4 is devoted to a discussion of the findings of other researchers and contains some numerical calculations concerning the previously mentioned orbital phenomena of LAGEOS and the Moon. Section 5 summarizes my results.


## 2. Analytical Calculation of the Orbital Effects byμ̇/μ

By defining(1)μ≡GM
at a given epoch t0, the acceleration of a test particle orbiting a central body experiencing a variation of μ is, to first order in t-t0,  (2)A=-μ(t)r2r̂≈-μr2[1+(μ̇μ)(t-t0)]r̂,
with μ̇≡μ̇∣t=t0. μ̇ will be assumed constant throughout the temporal interval of interest Δt=t-t0, as it is, for example, the case for most of the remaining lifetime of the Sun as a Main Sequence (MS) star [1]. Note that μ̇ can, in principle, be due to a variation of both the Newtonian gravitational constant G and the mass M of the central body, so that  (3)μ̇μ=ĠG+ṀM.
Moreover, while the orbital angular momentum is conserved, this does not happen for the energy.By limiting ourselves to realistic astronomical scenarios like our solar system, it is quite realistic to assume that(4)(μ̇μ)(t-t0)≪1
over most of its remaining lifetime: indeed, since Ṁ/M is of the order of (about 80% of such a mass-loss is due to the core nuclear burning, while the remaining 20% is due to average solar wind)10-14 yr-1 for the Sun [1], the condition (4) is satisfied for the remaining (the age of the present-day MS Sun is 4.58 Gyr, counted from its zero-age MS star model [1]) ≈ 7.58 Gyr before the Sun will approach the RGB tip in the Hertzsprung-Russell Diagram (HRD). Thus, I can treat it perturbatively with the standard methods of celestial mechanics.The unperturbed Keplerian ellipse at epocht0, assumed coinciding with the time of the passage at perihelion tp, is characterized by (5)r=a(1-ecosE),dt=(1-ecosEn)dE,cosf=cosE-e1-ecosE,sinf=1-e2sinE1-ecosE,
where a and e are the osculating semimajor axis and the eccentricity, respectively, which fix the size and the shape of the unchanging Keplerian orbit, n=μ/a3 is its unperturbed Keplerian mean motion, f is the true anomaly, reckoned from the pericentre, and E is the eccentric anomaly. This would be the path followed by the particle for any t>tp if the mass loss would suddenly cease at tp. Instead, the true path will be different because of the perturbation induced by μ̇ and the orbital parameters of the osculating ellipses approximating the real trajectory at each instant of time will slowly change in time.


### 2.1. The Semimajor Axis and the Eccentricity

The Gauss equation for the variation of the semimajor axisa is [19, 20]  (6)dadt=2n1-e2[eArsinf+Aτ(pr)],
where Ar and Aτ are the radial and transverse, that is, orthogonal to the direction of r̂, components, respectively, of the disturbing acceleration, and p=a(1-e2) is the semilatus rectum. In the present case  (7)A≡Ar=-μ̇r2(t-tp),
that is, there is an entirely radial perturbing acceleration. For μ̇<0, that is, a decrease in the body's GM, the total gravitational attraction felt by the test particle, given by (2), is reduced with respect to the epoch tp. In order to have the rate of the semimajor axis averaged over one (Keplerian) orbital revolution, (7) must be inserted into (6), evaluated onto the unperturbed Keplerian ellipse with (5), and finally integrated over ndt/2π from 0 to 2π because n/2π=1/PKep (see below). Note that, from (5), there can be obtained  (8)t-tp=E-esinEn.
As a result, I have  (9)〈dadt〉=-eπ(μ̇μ)a∫02π(E-esinE)sinE(1-ecosE)2dE=2(e1-e)(μ̇μ)a.
Note that if μ decreases, a gets reduced as well: 〈ȧ〉<0. This may be seemingly bizarre and counter-intuitive, but, as it will be shown later, it is not in contrast with the true orbital motion.The Gauss equation for the variation of the eccentricity is [19, 20]  (10)dedt=1-e2na{Arsinf+Aτ[cosf+1e(1-ra)]}.
For A=Ar, it reduces to  (11)dedt=(1-e22ae)dadt,
so that  (12)〈dedt〉=(1+e)(μ̇μ);
also the eccentricity gets smaller for μ̇<0.As a consequence of the found variations of the osculating semimajor axis and the eccentricity, the osculating orbital angular momentum per unit mass, defined byL2=μa(1-e2), remains constant: indeed, by using  (9) and  (12), it turns out  (13)〈dL2dt〉=μ〈ȧ〉(1-e2)-2μae〈ė〉=0.The osculating total energyℰ=-μ/2a decreases according to  (14)〈dℰdt〉=μ2a2〈ȧ〉=(e1-e)μ̇a.Moreover, the osculating Keplerian period(15)PKep≡2πn=2πa3μ,
which, by definition, yields the time elapsed between two consecutive perihelion crossings in absence of perturbation, that is, it is the time required to describe a fixed osculating Keplerian ellipse, decreases according to  (16)〈dPKepdt〉=32PKep〈ȧ〉a=6πeμ̇(1-e)(aμ)3/2.
As I will show, also such a result is not in contrast with the genuine orbital evolution.


### 2.2. The Pericentre, the Node, and the Inclination

The Gauss equation for the variation of the pericentreω is [19, 20]  (17)dωdt=1-e2nae[-Arcosf+Aτ(1+rp)sinf]-cosidΩdt,
where i and Ω are the the inclination and the longitude of the ascending node, respectively, which fix the orientation of the osculating ellipse in the inertial space. Since dΩ/dt and di/dt depend on the normal component Aν of the disturbing acceleration, which is absent in the present case, and A=Ar, I have  (18)〈dωdt〉=1-e22πe(μ̇μ)∫02π(E-esinE)(cosE-e)(1-ecosE)2dE=0:
the osculating ellipse does not change its orientation in the orbital plane, which, incidentally, remains fixed in the inertial space because Aν=0 and, thus, dΩ/dt=di/dt=0.


### 2.3. The Mean Anomaly

The Gauss equation for the mean anomalyℳ, defined as ℳ=n(t-tp), [19, 20] is  (19)dℳdt=n-2naArra-1-e2(dωdt+cosidΩdt).
It turns out that, since  (20)-2naArradt=2μ̇n3a3(E-esinE)dE,〈dℳdt〉=n+2π(μ̇μ);
the mean anomaly changes uniformly in time at a slower rate with respect to the unperturbed Keplerian case for μ̇<0.


### 2.4. Numerical Integration of the Equations of Motion and Explanation of the Seeming Contradiction with the Analytical Results

At first sight, the results obtained here may be rather confusing: if the gravitational attraction of the Sun reduces in time because of its mass loss, the orbits of the planets should expand (see the trajectory plotted in Figure1, numerically integrated with MATHEMATICA), while I obtained that the semimajor axis and the eccentricity undergo secular decrements. Moreover, I found that the Keplerian period PKep decreases, while one would expect that the orbital period increases.Figure 1
Black continuous line represents the true trajectory obtained by numerically integrating with MATHEMATICA the perturbed equations of motion in Cartesian coordinates over 2 years; the disturbing acceleration (2) has been adopted. The planet starts from the perihelion on the x axis. Just for illustrative purposes, a mass loss rate of the order of 10-2 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.8 have been chosen. Red dashed line represents unperturbed Keplerian ellipse at t=t0=tp. Blue dash-dotted line represents osculating Keplerian ellipse after the first perihelion passage. As can be noted, its semimajor axis and eccentricity are clearly smaller than those of the initial unperturbed ellipse. Note also that after 2 years the planet has not yet reached the perihelion as it would have done in absence of mass loss; that is, the true orbital period is longer than the Keplerian one of the osculating red ellipse.In fact, there is no contradiction, and my analytical results do yield us realistic information on the true evolution of the planetary motion. Indeed,a, e, and PKep refer to the osculating Keplerian ellipses which, at any instant, approximate the true trajectory; it, instead, is not an ellipse, not being bounded. Let us start at tp from the osculating pericentre of the Keplerian ellipse corresponding to chosen initial conditions: let us use a heliocentric frame with the x axis oriented along the osculating pericentre. After a true revolution, that is, when the true radius vector of the planet has swept an angular interval of 2π, the planet finds itself again on the x axis, but at a larger distance from the starting point because of the orbit expansion induced by the Sun's mass loss. It is not difficult to understand that the osculating Keplerian ellipse approximating the trajectory at this perihelion passage is oriented as before because there is no variation of the (osculating) argument of pericentre, but has smaller semimajor axis and eccentricity. And so on, revolution after revolution, until the perturbation theory can be applied, that is, until μ̇/μ(t-tp)≪1. In Figure 1 the situation described so far is qualitatively illustrated. Just for illustrative purposes I enhanced the overall effect by assuming |μ̇/μ|≈10-2 yr-1 for the Sun; the initial conditions for the planet correspond to an unperturbed Keplerian ellipse with a=1 AU, e=0.8 with the present-day value of the Sun's mass in one of its foci. It is apparent that the initial osculating red dashed ellipse has larger a and e with respect to the second osculating blue dash-dotted ellipse. Note also that the true orbital period, intended as the time elapsed between two consecutive crossings of the perihelion, is larger than the unperturbed Keplerian one of the initial red dashed osculating ellipse, which would amount to 1 year for the Earth: indeed, after 2 years the planet has not yet reached the perihelion for its second passage.Now, if I compute the radial changeΔr(E) in the osculating radius vector as a function of the eccentric anomaly E, I can gain useful insights concerning how much the true path has expanded after two consecutive perihelion passages. From the Keplerian expression of the Sun-planet distance  (21)r=a(1-ecosE),
one gets the radial component of the orbital perturbation expressed in terms of the eccentric anomaly E:  (22)Δr(E)=(1-ecosE)Δa-acosEΔe+aesinEΔE;
it agrees with the results obtained by, for example, Casotto in  [21]. Since  (23)Δa=-2aen(μ̇μ)(sinE-EcosE1-ecosE),Δe=-(1-e2)n(μ̇μ)(sinE-EcosE1-ecosE),ΔE=(Δℳ+sinEΔe1-ecosE)=1n(μ̇μ)[𝒜(E)+ℬ(E)+𝒞(E)],
with  (24)𝒜(E)=E2+2e(cosE-1)1-ecosE,ℬ(E)=(1-e2e)[1+e-(1+e)cosE-EsinE(1-ecosE)2],𝒞(E)=-(1-e2)sinE(sinE-ecosE)(1-ecosE)2,
it follows  (25)Δr(E)=an(μ̇μ)[𝒟(E)+ℱ(E)],
with  (26)𝒟(E)=e[-2(sinE-EcosE)+sinE[E2+2e(cosE-1)]1-ecosE-(1-e2)sin2E(sinE-ecosE)(1-ecosE)2],ℱ(E)=(1-e21-ecosE)×{[1+e-(1+e)cosE-EsinE1-ecosE]cosE(sinE-EcosE)+sinE[1+e-(1+e)cosE-EsinE1-ecosE]}.
From (25) and  (26) it turns out that for E>0  Δr(E) never vanishes; after one orbital revolution, that is, after that an angular interval of 2π has been swept by the (osculating) radius vector, a net increase of the radial (osculating) distance occurs according to (according to (25) and  (26), Δr(0)=0)  (27)Δr(2π)-Δr(0)=Δr(2π)=-2πna(μ̇μ)(1-e).
This analytical result is qualitatively confirmed by the differenceΔr(t) between the radial distances obtained from the solutions of two numerical integrations with MATHEMATICA of the equations of motion over 3 yr with and without μ̇/μ; the initial conditions are the same. (Strictly speaking, Δr and the quantity plotted in Figure 2 are different objects, but, as the following discussion will clarify, I can assume that, in practice, they are the same.)For illustrative purposes I used a=1 AU, e=0.01, μ̇/μ=-0.1 yr-1. The result is depicted in Figure 2. Note also that  (25) and  (26) tell us that the shift at the aphelion is (28)Δr(π)=12(1+e1-e)Δr(2π),
in agreement with Figure 1 where it is 4.5 times larger than the shift at the perihelion.Figure 2
DifferenceΔr(t) between the radial distances obtained from the solutions of two numerical integrations with MATHEMATICA of the equations of motion over 3 years with and without μ̇/μ; the initial conditions are the same. Just for illustrative purposes a mass loss rate of the order of -10-1 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.01 have been chosen. The cumulative increase of the Sun-planet distance induced by the mass loss is apparent.Since Figure1 tells us that the orbital period gets larger than the Keplerian one, it means that the true orbit must somehow remain behind with respect to the Keplerian one. Thus, a negative perturbation Δτ in the transverse direction must occur as well; see Figure 3.Figure 3
Radial and transverse perturbationsΔr and Δτ of the Keplerian radius vector (in blue); the presence of the transverse perturbation Δτ makes the real orbit (in red) lagging behind the Keplerian one.Let us now analytically compute it. According to  [21], it can be used (29)Δτ=asinE1-e2Δe+a1-e2ΔE+r(Δω+ΔΩcosi).
By recalling that, in the present case, ΔΩ=0 and using  (30)Δω=-1-e2ne(μ̇μ)[1+e-(1+e)cosE-EsinE1-ecosE],
it is possible to obtain from  (23) and  (30)  (31)Δτ(E)=an(μ̇μ)1-e2(1-ecosE)×[𝒢(E)+ℋ(E)+ℐ(E)+𝒥(E)+𝒦(E)],
with  (32)𝒢(E)=sinE(EcosE-sinE),ℋ(E)=(1-ecosE)e[(1+e)(cosE-1)+EsinE],ℐ(E)=E2+2e(cosE-1),𝒥(E)=sinE[(1-e2)(ecosE-sinE)1-ecosE],𝒦(E)=(1-e2e)[(1+e)(1-ecosE)-EsinE1-ecosE].
From  (31) and  (32) it turns out that for E>0Δτ(E) never vanishes; at the (osculating) time of perihelion passage  (33)Δτ(2π)-Δτ(0)=4π2na(μ̇μ)1+e1-e<0.
This means that when the Keplerian path has reached the perihelion, the perturbed orbit is still behind it. Such features are qualitatively confirmed by Figure 1.From a vectorial point of view, the radial and transverse perturbations to the Keplerian radius vectorr yield a correction:  (34)Δ=Δrr̂+Δττ̂,
so that  (35)rpert=r+Δ.
The length of Δ is  (36)Δ(E)=Δr(E)2+Δτ(E)2.
Equations (27) and (31) tell us that at perihelion it amounts to  (37)Δ(2π)=Δr(2π)1+4π2(1+e)(1-e)3.
The angle ξ between Δ and r is given by  (38)tanξ(E)=Δτ(E)Δr(E);
at perihelion it is  (39)tanξ(2π)=-2π1+e(1-e)3/2,
that is, ξ is close to -90 deg; for the Earth it is -81.1 deg. Thus, the difference δ between the lengths of the perturbed radius vector rpert and the Keplerian one r at a given instant amounts to about  (40)δ≈Δcosξ;
if fact, this is precisely the quantity determined over 3 years by the numerical integration of Figure 2. At the perihelion I have  (41)δ=Δr(2π)1+4π2(1+e)(1-e)3cosξ;
since for the Earth  (42)1+4π2(1+e)(1-e)3cosξ=1.0037,
it holds  (43)δ≈Δr(2π).
This explains why Figure 2 gives us just Δr.Since the approximate calculations of other researchers often refer to circular orbits, and in view of the fact that when a Sun-like star evolves into a giant tidal interactions circularize the orbit of a planet [23, 24], it is interesting to consider also such limiting case in which other nonsingular osculating orbital elements must be adopted. (This fact has been quantitatively proven by the observation of convective binary stars [22].) The eccentricity and the pericentre lose their meaning: thus, it is not surprising that  (12), although formally valid for e→0, yields a meaningless result; that is, the eccentricity would become negative. Instead, the semimajor axis is still valid and  (9) predicts that 〈ȧ〉=0 for e→0. The constancy of the osculating semimajor axis is not in contrast with the true trajectory, as clearly showed by Figure 4. Again, the true orbital period is larger than the Keplerian one which, contrary to the eccentric case, remains fixed. Since 𝒟(E)=0 for e=0 and ℱ(2π)|e=0=-2π,ℱ(0)|e=0=0, the radial shift per revolution is  (44)Δr(2π)|e=0=-2πna(μ̇μ).
Also in this case the secular increase of the radial distance is present, as qualitatively shown by Figure 5. Concerning Δτ, after 2π it is  (45)Δτ(2π)=4π2na(μ̇μ);
also in this case, the orbital period is larger than the unperturbed one.Figure 4
Black continuous liner represents true trajectory obtained by numerically integrating with MATHEMATICA the perturbed equations of motion in Cartesian coordinates over 2 years; the disturbing acceleration (2) has been adopted. The planet starts from a point on the x axis. Just for illustrative purposes, a mass loss rate of the order of -10-2 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.0 have been chosen. Red dashed line represents unperturbed Keplerian circle at t=t0. Blue dash-dotted line represents osculating Keplerian circle after the first x axis crossing. As can be noted, its semimajor axis and eccentricity are equal to those of the initial unperturbed circle. Note also that after 2 years the planet has not yet reached the x axis as it would have done in absence of mass loss.Figure 5
DifferenceΔr(t) between the radial distances obtained from the solutions of two numerical integrations with MATHEMATICA of the equations of motion over 3 years with and without μ̇/μ; the initial conditions are the same. Just for illustrative purposes a mass loss rate of the order of 10-1 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.0 have been chosen. The cumulative increase of the Sun-planet distance induced by the mass loss is apparent.


## 2.1. The Semimajor Axis and the Eccentricity

The Gauss equation for the variation of the semimajor axisa is [19, 20]  (6)dadt=2n1-e2[eArsinf+Aτ(pr)],
where Ar and Aτ are the radial and transverse, that is, orthogonal to the direction of r̂, components, respectively, of the disturbing acceleration, and p=a(1-e2) is the semilatus rectum. In the present case  (7)A≡Ar=-μ̇r2(t-tp),
that is, there is an entirely radial perturbing acceleration. For μ̇<0, that is, a decrease in the body's GM, the total gravitational attraction felt by the test particle, given by (2), is reduced with respect to the epoch tp. In order to have the rate of the semimajor axis averaged over one (Keplerian) orbital revolution, (7) must be inserted into (6), evaluated onto the unperturbed Keplerian ellipse with (5), and finally integrated over ndt/2π from 0 to 2π because n/2π=1/PKep (see below). Note that, from (5), there can be obtained  (8)t-tp=E-esinEn.
As a result, I have  (9)〈dadt〉=-eπ(μ̇μ)a∫02π(E-esinE)sinE(1-ecosE)2dE=2(e1-e)(μ̇μ)a.
Note that if μ decreases, a gets reduced as well: 〈ȧ〉<0. This may be seemingly bizarre and counter-intuitive, but, as it will be shown later, it is not in contrast with the true orbital motion.The Gauss equation for the variation of the eccentricity is [19, 20]  (10)dedt=1-e2na{Arsinf+Aτ[cosf+1e(1-ra)]}.
For A=Ar, it reduces to  (11)dedt=(1-e22ae)dadt,
so that  (12)〈dedt〉=(1+e)(μ̇μ);
also the eccentricity gets smaller for μ̇<0.As a consequence of the found variations of the osculating semimajor axis and the eccentricity, the osculating orbital angular momentum per unit mass, defined byL2=μa(1-e2), remains constant: indeed, by using  (9) and  (12), it turns out  (13)〈dL2dt〉=μ〈ȧ〉(1-e2)-2μae〈ė〉=0.The osculating total energyℰ=-μ/2a decreases according to  (14)〈dℰdt〉=μ2a2〈ȧ〉=(e1-e)μ̇a.Moreover, the osculating Keplerian period(15)PKep≡2πn=2πa3μ,
which, by definition, yields the time elapsed between two consecutive perihelion crossings in absence of perturbation, that is, it is the time required to describe a fixed osculating Keplerian ellipse, decreases according to  (16)〈dPKepdt〉=32PKep〈ȧ〉a=6πeμ̇(1-e)(aμ)3/2.
As I will show, also such a result is not in contrast with the genuine orbital evolution.


## 2.2. The Pericentre, the Node, and the Inclination

The Gauss equation for the variation of the pericentreω is [19, 20]  (17)dωdt=1-e2nae[-Arcosf+Aτ(1+rp)sinf]-cosidΩdt,
where i and Ω are the the inclination and the longitude of the ascending node, respectively, which fix the orientation of the osculating ellipse in the inertial space. Since dΩ/dt and di/dt depend on the normal component Aν of the disturbing acceleration, which is absent in the present case, and A=Ar, I have  (18)〈dωdt〉=1-e22πe(μ̇μ)∫02π(E-esinE)(cosE-e)(1-ecosE)2dE=0:
the osculating ellipse does not change its orientation in the orbital plane, which, incidentally, remains fixed in the inertial space because Aν=0 and, thus, dΩ/dt=di/dt=0.


## 2.3. The Mean Anomaly

The Gauss equation for the mean anomalyℳ, defined as ℳ=n(t-tp), [19, 20] is  (19)dℳdt=n-2naArra-1-e2(dωdt+cosidΩdt).
It turns out that, since  (20)-2naArradt=2μ̇n3a3(E-esinE)dE,〈dℳdt〉=n+2π(μ̇μ);
the mean anomaly changes uniformly in time at a slower rate with respect to the unperturbed Keplerian case for μ̇<0.


## 2.4. Numerical Integration of the Equations of Motion and Explanation of the Seeming Contradiction with the Analytical Results

At first sight, the results obtained here may be rather confusing: if the gravitational attraction of the Sun reduces in time because of its mass loss, the orbits of the planets should expand (see the trajectory plotted in Figure1, numerically integrated with MATHEMATICA), while I obtained that the semimajor axis and the eccentricity undergo secular decrements. Moreover, I found that the Keplerian period PKep decreases, while one would expect that the orbital period increases.Figure 1
Black continuous line represents the true trajectory obtained by numerically integrating with MATHEMATICA the perturbed equations of motion in Cartesian coordinates over 2 years; the disturbing acceleration (2) has been adopted. The planet starts from the perihelion on the x axis. Just for illustrative purposes, a mass loss rate of the order of 10-2 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.8 have been chosen. Red dashed line represents unperturbed Keplerian ellipse at t=t0=tp. Blue dash-dotted line represents osculating Keplerian ellipse after the first perihelion passage. As can be noted, its semimajor axis and eccentricity are clearly smaller than those of the initial unperturbed ellipse. Note also that after 2 years the planet has not yet reached the perihelion as it would have done in absence of mass loss; that is, the true orbital period is longer than the Keplerian one of the osculating red ellipse.In fact, there is no contradiction, and my analytical results do yield us realistic information on the true evolution of the planetary motion. Indeed,a, e, and PKep refer to the osculating Keplerian ellipses which, at any instant, approximate the true trajectory; it, instead, is not an ellipse, not being bounded. Let us start at tp from the osculating pericentre of the Keplerian ellipse corresponding to chosen initial conditions: let us use a heliocentric frame with the x axis oriented along the osculating pericentre. After a true revolution, that is, when the true radius vector of the planet has swept an angular interval of 2π, the planet finds itself again on the x axis, but at a larger distance from the starting point because of the orbit expansion induced by the Sun's mass loss. It is not difficult to understand that the osculating Keplerian ellipse approximating the trajectory at this perihelion passage is oriented as before because there is no variation of the (osculating) argument of pericentre, but has smaller semimajor axis and eccentricity. And so on, revolution after revolution, until the perturbation theory can be applied, that is, until μ̇/μ(t-tp)≪1. In Figure 1 the situation described so far is qualitatively illustrated. Just for illustrative purposes I enhanced the overall effect by assuming |μ̇/μ|≈10-2 yr-1 for the Sun; the initial conditions for the planet correspond to an unperturbed Keplerian ellipse with a=1 AU, e=0.8 with the present-day value of the Sun's mass in one of its foci. It is apparent that the initial osculating red dashed ellipse has larger a and e with respect to the second osculating blue dash-dotted ellipse. Note also that the true orbital period, intended as the time elapsed between two consecutive crossings of the perihelion, is larger than the unperturbed Keplerian one of the initial red dashed osculating ellipse, which would amount to 1 year for the Earth: indeed, after 2 years the planet has not yet reached the perihelion for its second passage.Now, if I compute the radial changeΔr(E) in the osculating radius vector as a function of the eccentric anomaly E, I can gain useful insights concerning how much the true path has expanded after two consecutive perihelion passages. From the Keplerian expression of the Sun-planet distance  (21)r=a(1-ecosE),
one gets the radial component of the orbital perturbation expressed in terms of the eccentric anomaly E:  (22)Δr(E)=(1-ecosE)Δa-acosEΔe+aesinEΔE;
it agrees with the results obtained by, for example, Casotto in  [21]. Since  (23)Δa=-2aen(μ̇μ)(sinE-EcosE1-ecosE),Δe=-(1-e2)n(μ̇μ)(sinE-EcosE1-ecosE),ΔE=(Δℳ+sinEΔe1-ecosE)=1n(μ̇μ)[𝒜(E)+ℬ(E)+𝒞(E)],
with  (24)𝒜(E)=E2+2e(cosE-1)1-ecosE,ℬ(E)=(1-e2e)[1+e-(1+e)cosE-EsinE(1-ecosE)2],𝒞(E)=-(1-e2)sinE(sinE-ecosE)(1-ecosE)2,
it follows  (25)Δr(E)=an(μ̇μ)[𝒟(E)+ℱ(E)],
with  (26)𝒟(E)=e[-2(sinE-EcosE)+sinE[E2+2e(cosE-1)]1-ecosE-(1-e2)sin2E(sinE-ecosE)(1-ecosE)2],ℱ(E)=(1-e21-ecosE)×{[1+e-(1+e)cosE-EsinE1-ecosE]cosE(sinE-EcosE)+sinE[1+e-(1+e)cosE-EsinE1-ecosE]}.
From (25) and  (26) it turns out that for E>0  Δr(E) never vanishes; after one orbital revolution, that is, after that an angular interval of 2π has been swept by the (osculating) radius vector, a net increase of the radial (osculating) distance occurs according to (according to (25) and  (26), Δr(0)=0)  (27)Δr(2π)-Δr(0)=Δr(2π)=-2πna(μ̇μ)(1-e).
This analytical result is qualitatively confirmed by the differenceΔr(t) between the radial distances obtained from the solutions of two numerical integrations with MATHEMATICA of the equations of motion over 3 yr with and without μ̇/μ; the initial conditions are the same. (Strictly speaking, Δr and the quantity plotted in Figure 2 are different objects, but, as the following discussion will clarify, I can assume that, in practice, they are the same.)For illustrative purposes I used a=1 AU, e=0.01, μ̇/μ=-0.1 yr-1. The result is depicted in Figure 2. Note also that  (25) and  (26) tell us that the shift at the aphelion is (28)Δr(π)=12(1+e1-e)Δr(2π),
in agreement with Figure 1 where it is 4.5 times larger than the shift at the perihelion.Figure 2
DifferenceΔr(t) between the radial distances obtained from the solutions of two numerical integrations with MATHEMATICA of the equations of motion over 3 years with and without μ̇/μ; the initial conditions are the same. Just for illustrative purposes a mass loss rate of the order of -10-1 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.01 have been chosen. The cumulative increase of the Sun-planet distance induced by the mass loss is apparent.Since Figure1 tells us that the orbital period gets larger than the Keplerian one, it means that the true orbit must somehow remain behind with respect to the Keplerian one. Thus, a negative perturbation Δτ in the transverse direction must occur as well; see Figure 3.Figure 3
Radial and transverse perturbationsΔr and Δτ of the Keplerian radius vector (in blue); the presence of the transverse perturbation Δτ makes the real orbit (in red) lagging behind the Keplerian one.Let us now analytically compute it. According to  [21], it can be used (29)Δτ=asinE1-e2Δe+a1-e2ΔE+r(Δω+ΔΩcosi).
By recalling that, in the present case, ΔΩ=0 and using  (30)Δω=-1-e2ne(μ̇μ)[1+e-(1+e)cosE-EsinE1-ecosE],
it is possible to obtain from  (23) and  (30)  (31)Δτ(E)=an(μ̇μ)1-e2(1-ecosE)×[𝒢(E)+ℋ(E)+ℐ(E)+𝒥(E)+𝒦(E)],
with  (32)𝒢(E)=sinE(EcosE-sinE),ℋ(E)=(1-ecosE)e[(1+e)(cosE-1)+EsinE],ℐ(E)=E2+2e(cosE-1),𝒥(E)=sinE[(1-e2)(ecosE-sinE)1-ecosE],𝒦(E)=(1-e2e)[(1+e)(1-ecosE)-EsinE1-ecosE].
From  (31) and  (32) it turns out that for E>0Δτ(E) never vanishes; at the (osculating) time of perihelion passage  (33)Δτ(2π)-Δτ(0)=4π2na(μ̇μ)1+e1-e<0.
This means that when the Keplerian path has reached the perihelion, the perturbed orbit is still behind it. Such features are qualitatively confirmed by Figure 1.From a vectorial point of view, the radial and transverse perturbations to the Keplerian radius vectorr yield a correction:  (34)Δ=Δrr̂+Δττ̂,
so that  (35)rpert=r+Δ.
The length of Δ is  (36)Δ(E)=Δr(E)2+Δτ(E)2.
Equations (27) and (31) tell us that at perihelion it amounts to  (37)Δ(2π)=Δr(2π)1+4π2(1+e)(1-e)3.
The angle ξ between Δ and r is given by  (38)tanξ(E)=Δτ(E)Δr(E);
at perihelion it is  (39)tanξ(2π)=-2π1+e(1-e)3/2,
that is, ξ is close to -90 deg; for the Earth it is -81.1 deg. Thus, the difference δ between the lengths of the perturbed radius vector rpert and the Keplerian one r at a given instant amounts to about  (40)δ≈Δcosξ;
if fact, this is precisely the quantity determined over 3 years by the numerical integration of Figure 2. At the perihelion I have  (41)δ=Δr(2π)1+4π2(1+e)(1-e)3cosξ;
since for the Earth  (42)1+4π2(1+e)(1-e)3cosξ=1.0037,
it holds  (43)δ≈Δr(2π).
This explains why Figure 2 gives us just Δr.Since the approximate calculations of other researchers often refer to circular orbits, and in view of the fact that when a Sun-like star evolves into a giant tidal interactions circularize the orbit of a planet [23, 24], it is interesting to consider also such limiting case in which other nonsingular osculating orbital elements must be adopted. (This fact has been quantitatively proven by the observation of convective binary stars [22].) The eccentricity and the pericentre lose their meaning: thus, it is not surprising that  (12), although formally valid for e→0, yields a meaningless result; that is, the eccentricity would become negative. Instead, the semimajor axis is still valid and  (9) predicts that 〈ȧ〉=0 for e→0. The constancy of the osculating semimajor axis is not in contrast with the true trajectory, as clearly showed by Figure 4. Again, the true orbital period is larger than the Keplerian one which, contrary to the eccentric case, remains fixed. Since 𝒟(E)=0 for e=0 and ℱ(2π)|e=0=-2π,ℱ(0)|e=0=0, the radial shift per revolution is  (44)Δr(2π)|e=0=-2πna(μ̇μ).
Also in this case the secular increase of the radial distance is present, as qualitatively shown by Figure 5. Concerning Δτ, after 2π it is  (45)Δτ(2π)=4π2na(μ̇μ);
also in this case, the orbital period is larger than the unperturbed one.Figure 4
Black continuous liner represents true trajectory obtained by numerically integrating with MATHEMATICA the perturbed equations of motion in Cartesian coordinates over 2 years; the disturbing acceleration (2) has been adopted. The planet starts from a point on the x axis. Just for illustrative purposes, a mass loss rate of the order of -10-2 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.0 have been chosen. Red dashed line represents unperturbed Keplerian circle at t=t0. Blue dash-dotted line represents osculating Keplerian circle after the first x axis crossing. As can be noted, its semimajor axis and eccentricity are equal to those of the initial unperturbed circle. Note also that after 2 years the planet has not yet reached the x axis as it would have done in absence of mass loss.Figure 5
DifferenceΔr(t) between the radial distances obtained from the solutions of two numerical integrations with MATHEMATICA of the equations of motion over 3 years with and without μ̇/μ; the initial conditions are the same. Just for illustrative purposes a mass loss rate of the order of 10-1 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.0 have been chosen. The cumulative increase of the Sun-planet distance induced by the mass loss is apparent.


## 3. The General Relativistic Case

The field equations of general relativity are nonlinear, but in the slow-motion(|β|=|v|/c=1) and weak-field (U/c2=1) approximation they get linearized resembling to the linear equations of the Maxwellian electromagnetism; here v and U are the magnitudes of the typical velocities and the gravitational potential of the problem under consideration. This scenario is known as gravitoelectromagnetism [25, 26]. In this case the space-time metric is given by  (46)ds2=(1-2Φc2)c2dt2+4c(H·dr)dt-(1+2Φc2)δijdxidxj,
where, far from the source, the dominant contributions to the gravitoelectromagnetic potentials can be expressed as  (47)Φ=μr,H=GcJ×rr3.
Here J is the proper angular momentum of the central body of mass M and r is so that r≫GM/c2 and r≫J/(Mc); c is the speed of light in vacuum. For a nonstationary source the geodesic equations of motion yield [27], among other terms, -βi(3-β2)Φ,0,i=1,2,3 which, to order 𝒪(c-2), reduces to  (48)A=-3μ̇c2vr=-3GMc2(ĠG+ṀM)vr.
For μ̇<0 such a perturbing acceleration is directed along the velocity of the test particle. Although of no practical interest, being of the order of 10-24 m s-2 in the case of a typical Sun-planet system with Ṁ/M=-9×10-14 yr-1, I will explicitly work out the orbital effects of  (48); the effects of the temporal variations of J have already been worked out elsewhere [27, 28]. Also in this case I will use the Gauss perturbative case. Since the radial and transverse components of the unperturbed velocity are  (49)vr=naesinf1-e2,vτ=na(1+ecosf)1-e2,
the radial and transverse components of  (48), evaluated onto the unperturbed Keplerian orbit, are  (50)Ar=-3μ̇c2nesinE(1-ecosE)2,Aτ=-3μ̇c2n1-e2(1-ecosE)2.
After lengthy calculations they yield  (51)〈dadt〉=-6μ̇c2(21-e2-1),〈dedt〉=3μ̇μc2a4e(2-e)(1-11-e2).
Contrary to the classical case, now both the osculating semimajor axis and the eccentricity increase for μ̇<0. It turns out that the pericentre and the mean anomaly do not secularly precess. Also in this case the inclination and the node are not affected The qualitative features of the motion with the perturbation  (48) are depicted in Figure 6 in which the magnitude of the relativistic term has been greatly enhanced for illustrative purposes.Figure 6
Black continuous line represents true trajectory obtained by numerically integrating over 3 years the equations of motion perturbed by  (48). The planet starts from the perihelion on the x axis. Just for illustrative purposes, a factor -3μ̇/c2 of the order of 5×10-2 AU yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.8 have been chosen. Red dashed line represents unperturbed Keplerian ellipse at t=t0=tp. Blue dash-dotted line represents osculating Keplerian ellipse after the first perihelion passage. As can be noted, its semimajor axis and eccentricity are larger than those of the initial unperturbed ellipse.


## 4. Discussion of Other Approaches and Numerical Calculations

Here I will briefly review some of the results obtained by others by comparing with ours.Hadjidemetriou in  [14] uses a tangential perturbing acceleration proportional to the test particle's velocity v,  (52)A=-12[μ̇μ(t)]v,
and a different perturbative approach by finding that, for a generic mass loss, the semimajor axis secularly increases and the eccentricity remains constant. In fact, with the approach followed here it would be possible to show that, to first order in (μ̇/μ)(t-t0), 〈ȧ〉=-(μ̇/μ)a and 〈ė〉=0 and that the true orbit is expanding, although in a different way with respect to  (2) as depicted by Figure 7 in which the magnitude of the mass-loss has been exaggerated for better showing its orbital effects. However, it must be noted that a term like  (52) is inadmissible in any relativistic theory of gravitation because it violates the Lorentz invariance. Indeed, this fact is explicitly shown for general relativity by Bini et al. in  [27] where the full equations of motion of a test particle in a nonstationary gravitoelectromagnetic field are worked out (see [27, equation (14)]). In deriving them it is admitted that, in general, Φ=Φ(t,r), but no gravitoelectric terms like  (52) occur. Instead,  (2) is compatible with equation (14) of [27].Figure 7
Black continuous line represents true trajectory obtained by numerically integrating with MATHEMATICA the perturbed equations of motion in Cartesian coordinates over 2 years; the disturbing acceleration  (52) has been used. The planet starts from the perihelion on the x axis. Just for illustrative purposes, a mass loss rate of the order of -10-1 yr-1 has been adopted for the Sun; for the planet initial conditions corresponding to a=1 AU, e=0.8 have been chosen. Red dashed line represents unperturbed Keplerian ellipse at t=t0=tp. Blue dash-dotted line represents osculating Keplerian ellipse after the first perihelion passage. As can be noted, its semimajor axis is larger than that of the initial unperturbed ellipse, while the eccentricity remaines constant. Note also that after 2 years the planet has not yet reached the perihelion as it would have done in absence of mass loss.Schröder and Smith in [1], by assuming the conservation of the angular momentum, derive the orbital expansion by means of equations valid, instead, for orbits with constant radius only, that is, v2/r=μ(t)/r2 and L=vr. Then, they assume that not only v but also r vary and put v(t)=μ(t)/r, which is, instead, valid for circular orbits of constant radius only, into L=v(t)r(t)=vr getting μ(t)r(t)=μr, where in my notation r and μ refer to the initial epoch t0. With such an approach they obtain an expanded terrestrial orbit up to about 2 times larger than mine.Noerdlinger in [4], following Jeans [12] and Kevorkian and Cole [18], assumes for the variation of a quantity identified by him with the semimajor axis the following expression:  (53)a(t)μ(t)=aμ:
thus, his semimajor axis gets larger. Note that such an equation is the same obtained by  [1]. By assuming a variation of μ linear in time  (53) would yield an increase of a according to  (54)ȧ=-(μ̇μ)a>0;
Compare with my  (9). As a consequence of the constancy of L2=μ(t)a(t)[1-e(t)2] and of  (53) he obtains that the eccentricity remains constant, that is,  (55)ė=0;
Compare with my  (12). Moreover, another consequence of  (53) obtained by Noerdlinger in  [4] is that the Keplerian period increases as  (56)PKep(t)PKep=[μμ(t)]2;
Compare with my  (16). Should the quantities dealt with by Noerdlinger are to be identified with the usual osculating Keplerian elements, his results would be incompatible with the real dynamics of a test particle in the field of a linearly mass-losing body, as I have shown. The quantity obtained by me which exhibits the closest resemblance with (54) seems to be the secular variation of Δr(2π) for e=0. Apart from matters of interpretation, the quantitative results are different. Indeed, I obtain for the Earth a secular variation of the semimajor axis of -2×10-4 m yr-1 and a shift in the radial position along the fixed line of the apsides of +1.3×10-2 m yr-1, while Noerdlinger in [4] gets a secular rate of his semimajor axis, identified with the Astronomical Unit, of about +1×10-2 m yr-1. Note that Noerdlinger uses for the Sun Ṁ/M=-9×10-14 yr-1 as in the present work.Krasinski and Brumberg in [2] deal, among other things, with the problem of a mass-losing Sun in the framework of the observed secular increase of the Astronomical Unit for which, starting with an equation of motion like  (2), they obtain an equation like  (54). A mass-loss rate of -3×10-14 yr-1, considered somewhat underestimated by Noerdlinger [4], yields an increase of the Astronomical Unit of 3×10-3 m yr-1. With such a value for μ̇/μ I would obtain a decrease of the semimajor axis of -7×10-5 m yr-1 and an increase in r of +4×10-3 m yr-1.Concerning the observationally determined increase of the Astronomical Unit, more recent estimates from processing of huge planetary data sets by Pitjeva point towards a rate of the order of10-2 m yr-1 [29, 30]. It may be noted that my result for the secular variation of the terrestrial radial position on the line of the apsides would agree with such a figure by either assuming a mass loss by the Sun of just -9×10-14 yr-1 or a decrease of the Newtonian gravitational constant Ġ/G≈-1×10-13 yr-1. Such a value for the temporal variation of G is in agreement with recent upper limits from Lunar Laser Ranging [31] Ġ/G=(2±7)×10-13 yr-1. This possibility is envisaged in  [32] whose authors use ȧ/a=-Ġ/G by speaking about a small radial drift of -(6±13)×10-2 m yr-1 in an orbit at 1 AU. I, now, apply my analysis to the daily decrement of the semimajor axis of LAGEOS [8] whose geocentric relevant orbital parameters are a=12270 km, e=0.0045. From  (9) it turns out that a secular decrease of a of the order of ≈1 mm d-1 could only be induced by μ̇⊕/μ⊕=-3×10-6 yr-1. Clearly, it cannot be due to a variation of the Earth's mass M; if it had to be attributed to a variation of G, it would be orders of magnitude too large with respect to the bounds in  [31, 32]. Adopting  (27) does not substantially alter the situation because the required μ̇⊕/μ⊕ would be only two orders of magnitude smaller. Now I look at the anomalous increase of the eccentricity of the lunar orbit amounting to ėMoon=(0.9±0.3)×10-11 yr-1 [9]. Such a figure and  (12) yield μ̇⊕/μ⊕=8.5×10-12 yr-1; it is one order of magnitude larger than the present-day bounds on Ġ/G [31, 32].


## 5. Conclusions

I started in the framework of the two-body Newtonian dynamics by using a radial perturbing acceleration linear in time and straightforwardly treated it with the standard Gaussian scheme. I found that the osculating semimajor axisa, the eccentricity e and the mean anomaly ℳ secularly decrease while the argument of pericentre ω remains unchanged; the longitude of the ascending node Ω and the inclination i are not affected by the phenomenon considered. The radial distance from the central body, taken on the fixed line of the apsides, experiences a secular increase Δr. For the Earth, such an effect amounts to about 1.3 cm yr-1. By numerically integrating the equations of motion in Cartesian coordinates I found that the real orbital path expands after every revolution, the line of the apsides does not change, and the apsidal period is larger than the unperturbed Keplerian one. I have also clarified that such results are not in contrast with those analytically obtained for the Keplerian orbital elements which, indeed, refer to the osculating ellipses approximating the true trajectory at each instant. I also computed the orbital effects of a secular variation of the Sun's mass in the framework of the general relativistic linearized gravitoelectromagnetism which predicts a perturbing gravitoelectric tangential force proportional to v/r. I found that both the semimajor axis and the eccentricity secularly increase; the other Keplerian elements remain constant. Such effects are completely negligible in the present and future evolution of, for example, the solar system.As a suggestion to other researchers, it would be very important to complement my analytical two-body calculation by performing simultaneous long-term numerical integrations of the equations of motion of all the major bodies of the solar system by including a mass-loss term in the dynamical force models as well to see if the N-body interactions in presence of such an effect may substantially change the picture outlined here. It would be important especially in the RGB phase in which the inner regions of the solar system should dramatically change.



---

*Source: 261249-2010-03-09.xml*