# Long term behavior of a hypothetical planet in a highly eccentric orbit

## 1 Motivation

The information on Earth’s past climate obtained from deep sea sediments and polar ice cores indicates that the global temperature oscillated with a small amplitude around a constant mean value for millions of years until about 3.2 million years ago. From then on the mean temperature decreased stepwise and the amplitude of the variations increased. For the last one million years Earth’s climate was characterized by a distinct 100 kyr periodicity in the advancing and retreating of the polar ice sheets, known as glacial–interglacial cycles. This Ice Age period ended abruptly 11.5 kyr ago. During the following time interval, the Holocene, the global mean temperature quickly recovered to a value comparable with that observed 3.2 Myr ago . Wölfli et al.  proposed that during all this time an object of planetary size, called Z, existed in a highly eccentric orbit, and showed that Z could be the cause of these changes of Earth’s climate. Because of the postulated small perihelion distance, Z was heated up by solar radiation and tidal forces so that it was surrounded by a gas cloud. Whenever the Earth crossed this cloud, molecules activating a Greenhouse effect were produced in the upper atmosphere in an amount sufficient to transiently enhance the mean surface temperature on Earth. Very close flyby events even resulted in earthquakes and volcanic activities. In extreme cases, a rotation of the entire Earth relative to its rotation axis occured in response to the transient strong gravitational interaction. These polar shifts took place with a frequency of about one in one million years on the average. The first of them was responsible for the major drop in mean temperature, whereas the last one terminated Earth’s Pleistocene Ice Age 11.5 kyr ago. At present, planet Z does not exist any more. The high eccentricity of Z’s orbit corresponds to a small orbital angular momentum. This could be transferred to one of the inner planets during a close encounter, so that Z plunged into the sun. Alternatively, and more likely, Z approached the Earth to less than the Roche limit during the last polar shift event. In this case it was split into several parts which lost material at an accelerated rate because of the reduced escape velocity, so that eventually all of these fragments evaporated during the Holocene.

Here, we describe the method used to calculate the motion of such a hypothetical object in the presence of the other planets of the solar system, and its close encounters with the Earth. The calculations neglect possible losses of mass, orbital energy and angular momentum of Z due to solar irradiation and tidal effects. We also disregard the disappearance of Z following the last polar shift event. Consequences of these effects are discussed in ref. .

## 2 The method

In order to study the motion of an additional planet Z in the gravitational field of the sun and the other planets, a set of coupled ordinary differential equations (ODE) has to be integrated numerically starting at a given time with the known or assumed set of orbital parameters of all celestial bodies of interest. For the calculation we used the Pascal program Odeint which is based on the Bulirsch-Stoer method. It includes the modified midpoint algorithm with adaptive stepsize control and Richardson’s deferred approach to the limit . The start values for the known planets were taken from ref. , where they are given in barycentric rectangular coordinates and velocity components referred to the mean equator and equinox of J2000.0, the standard epoch, corresponding to 2000 January 1.5. All values are given in astronomical units (AU) and AU/day, respectively. To obtain the heliocentric coordinates referred to the ecliptic, we subtracted the coordinates of the sun and rotated the resulting values for all planets by 2326’21.448”, the angle between the mean equatorial plane and the ecliptic of J2000.0. The Earth and the Moon were treated separately. The masses of the planets are from ref. . For reasons explained in ref.  we assume that Z was a mars–like object with $`M_Z=0.11M_E`$ and that its hypothetical orbit at the epoch J2000.0 is determined by the following heliocentric parameters:

semi-major axis: $`a`$ $`=0.978`$   
numerical eccentricity: $`e`$ $`=0.973`$   
inclination: $`i`$ $`=0^{}`$   
longitude of the perihelion: $`\mathrm{\Omega }`$ $`=0^{}`$   
argument of the perihelion: $`\omega `$ $`=0^{}`$   
mean anomaly: $`M`$ $`=270.0^{}`$   
The calculation was a classical point-mass integration without relativistic corrections. In order to save computing time we ignored the influences from the three outermost planets Uranus, Neptune and Pluto, and restricted the numerical accuracy per integration step, the tolerance level $`eps`$, to a value of $`10^{13}`$ . Several tests were made to check whether these simplifications are acceptable or not. First of all, we evaluated the osculating orbital parameters for the Earth over the past 300 kyr without Z, but for three different tolerance levels, $`eps=10^{13},\mathrm{\hspace{0.33em}10}^{15}`$ and $`10^{16}`$, and found that the positions of all planets except that of Mercury were reproducible to within less than 100 km. A comparison of the Earth’s eccentricity and inclination variations with the corresponding values published by Berger  and transferred to the invariable plane of the solar system by Quinn et al  also showed good agreement. We also determined the total angular momentum of the solar system and found a negligible small linear change at the eleventh digit of its value. All calculations were performed with a 133 MHz PC. Including Z, about 15 h were required to cover a time span of 10 kyr with an accuracy of $`10^{13}`$. The sequence of close encounters of Z with the inner planets and the Moon amplify the errors so that the calculated result represents a possible orbit only. Integrations with Z were performed both, forward (+300 kyr) and backward (-450 kyr) in time relative to J2000.0 in order to demonstrate that the behavior of the orbital parameters of Z and the Earth were not affected by the change in time direction.

## 3  Results

Fig. 1

shows the time dependent variations of the osculating orbital parameters of Z (semi-major axis, eccentricity, inclination, ascending node and argument of the perihelion) over a time period of 750 kyr. The inclination is defined here as the angle between Z’s orbital plane and the invariable plane, which is perpendicular to the orbital angular momentum of the solar system. It is close to the orbital plane of Jupiter. The movement of Z is strongly influenced by the inner planets and to a lesser extent by Jupiter. Of interest is the behavior of the inclination, which on average oscillates with a periodicity of about 7 kyr only. This period is more than an order of magnitude shorter than that of Earth’s inclination which is essentially determined by Jupiter and, without Z included in the calculation, amounts to 100 kyr. The corresponding osculating parameters for Earth’s orbit are displayed in Fig. 2;

as mentioned above, they are only marginally disturbed by the presence of Z. Sudden jumps in one or several orbital parameter values of Z indicate close encounters with one of the inner planets.

In order to find out how often Z approached the Earth to distances close enough to influence its climate or even to provoke a polar shift , we determined the distance between Z and the Earth for each integration step and fitted a parabolic function into the values close to the distance of closest approach. The minimum of this quadratic function was then identified with this distance. The upper panel of Fig. 3

shows the result of this evaluation for the Z–Earth system over the whole time range considered here. Plotted are all encounters with distances of less than 0.02 AU $`=310^6`$ km, as indicated by the endpoints of each vertical line. The lower panel of Fig. 3 shows details of the irregular structure of these encounters which are the result of the complex time dependence of the coordinates of Z and Earth. Not surprisingly, the encounter frequency is enhanced whenever the two inclinations nearly coincide. Fig. 4 (left side)

shows that Z approaches the Earth many times to less than Moon’s distance. As explained in ref.  these flyby events can excite strong earthquakes and volcanic activities. For distances smaller than 30’000 km such encounters could even result in a rotation of the Earth by as much as 20 with respect to the direction of the invariant angular momentum. Flyby events having distances larger than about the Moon-Earth distance are harmless in this respect, but they still may influence Earth’s climate. In ref.  we have shown that Z was surrounded by a gas cloud which had an estimated radius of about 2.8 Mio. km at the intersection point with Earth’s orbit. The interaction of this cloud with Earth’s atmosphere produced Greenhouse gases in sufficient amounts to significantly increase the global temperature.

Close encounters with the Earth also imply close encounters with the Moon. These are plotted on the right side of Fig. 4. Since the mass of the Moon is much smaller than that of Z, recoil effects are much larger than in the case of the Earth and, therefore, may significantly influence the Moon’s orbit. In fact, Fig. 5.

shows events in which the semi-major axis $`a`$ suddenly changes at some given time by up to 9% relative to the mean value. Since, according to Kepler’s third law, the orbital period of the Moon is proportional to $`a^{3/2}`$, the orbital angular frequency $`\omega _{orb}`$ also changes, so that the apparent rotational frequency $`\mathrm{\Omega }=\omega _{rot}\omega _{orb}`$ becomes different from zero, assuming that $`\omega _{rot}`$, the rotational angular frequency, is not influenced by such an event. In ref.  we propose that the last close encounter of Z with the Earth took place only 11.5 kyr ago, so that Moon’s orbit could also have changed at that time. Therefore, the question arises whether the tidal friction on the Moon is large enough to stop $`\mathrm{\Omega }`$ within the Holocene. In the appendix we show that this is likely to be the case.

## 4 Conclusion

The calculations presented here show that an object of planetary size in a highly eccentric orbit approaches the Earth with sufficient frequency to influence its climate and even to produce polar shifts, the last of which terminated Earth’s Ice Age period, as explained in ref. . Although the pattern of these approaches compares well with the observed pattern, we have to point out that for various reasons a one to one correspondence between “theory” and “observation” cannot be expected: First, since Z no longer exists at present, we have to calculate its long-term behavior on the basis of assumed orbital parameters and on estimations of its mass. Second, Z lost substantial amounts of mass, orbital energy and angular momentum each time it passed through the perihelion. These effects are neglected in our point–mass model calculation. They may significantly alter the orbit of Z, and have to be included in an attempt to find out whether the disturbances of the orbits of the inner planets due to Z are within the boundaries set by present day observations. In ref. 3, we also point out that Z disappeared from the solar system either during the proposed polar shift event 11.5 kyr ago or later on during the Holocene. A detailed study of the mechanisms responsible for this removal is another important task.

## Appendix: Synchronising Moon’s rotation

Since a close encounter between Z and the Moon can lead to an apparent rotational frequency $`\mathrm{\Omega }=\omega _{rot}\omega _{orb}`$ of the Moon which no longer vanishes on the average, it is important to know how fast tidal friction diminishes $`\mathrm{\Omega }`$. The tidal force field is parallel to the direction Moon-Earth and has the value

$$F=\frac{2zGM_ER_M}{R^3}.$$
(1)   
$`G`$ is the gravitational constant, $`M_E`$ the mass of the Earth, $`R_M`$ the radius of the Moon, $`R`$ the distance between the centers of Earth and Moon, and $`z`$ the cartesian coordinate in the direction Moon–Earth, measured from Moon’s center. On Moon’s surface $`z=R_M\mathrm{cos}(\gamma )`$, where $`\gamma `$ is the angle between the $`z`$-axis and the direction to the point considered. Under the influence of $`F`$ and the gravitational acceleration $`g_M`$ on Moon’s surface, its shape will deviate from that of a sphere. The deformation is in first order given by

$$H(\gamma )=H_0\left(\frac{1}{3}+\mathrm{cos}(2\gamma )\right),$$
(2)   
In equilibrium $`H_0`$ becomes

$$H_0=\frac{GR_M^2M_E}{2R^3g_M}=6.5\text{m}.$$
(3)   
In a time dependent situation elastic tensions reduce the deformation. They will fade away gradually so that equilibrium is reached, say, with a relaxation time $`\tau `$. Assuming that the presence of an apparent rotation $`\mathrm{\Omega }`$ results in a deformation which lags behind with a phase difference $`\varphi =\mathrm{\Omega }\tau `$, the force field $`F`$ acting on this deformation will exert an angular moment $`D`$ which tends to stop the rotation $`\mathrm{\Omega }`$. The integration over Moon’s surface yields

$$D=\frac{\mathrm{sin}(2\varphi )}{2}\frac{16\pi }{15}\frac{G^2M_E^2R_M^6\rho }{R^6g_M}$$
(4)   
where $`\rho `$ is the density on Moon’s surface. The apparent rotation $`\mathrm{\Omega }`$ varies in time as

$$\frac{d\mathrm{\Omega }}{dt}=\frac{D}{\mathrm{\Xi }}=K\frac{\mathrm{sin}(2\varphi )}{2}$$
(5)   
For the inertial moment of the Moon $`\mathrm{\Xi }`$ we use the value of a homogeneous sphere, $`\mathrm{\Xi }=\frac{2}{5}M_MR_M^2`$, and set $`M_M=\frac{4\pi }{3}R_M^3\rho `$. Then $`K=\frac{M_E^2GR_M^3}{M_MR^6}=1.010^{16}`$ s<sup>-2</sup>. For $`\varphi =\mathrm{\Omega }\tau 1`$ the solution of this equation is given by

$$\mathrm{\Omega }(t)=\mathrm{\Omega }(0)\text{e}^{K\tau t}$$
(6)   
The true relaxation time of this deformation is open to discussion. If, for example, this time is assumed to be $`\tau =1`$ d = 86 400 s, then the decay constant for $`\mathrm{\Omega }`$ becomes

$$\tau _\mathrm{\Omega }=\frac{1}{K\tau }=3500\text{ yr}$$
(7)   
Assuming $`\mathrm{\Omega }=0.1\omega _{rot}`$, then the phase becomes $`\varphi =0.1\omega _{rot}\tau =0.023\text{ rad}=1.3^{}`$. Thus a rather small phase lag of the tidal deformation is sufficient to synchronise the lunar rotation within the Holocene. This assumed phase lag is smaller than the value of 2.16 inferred from astronomical data regarding the actual lunar bulge . The fact that according to Eq. 6 $`\mathrm{\Omega }(t)`$ never stops, is an artefact of a model with a single decay constant. The real dynamics also involves slow relaxations, which terminate the synchronisation in a finite time.

## Acknowledgement

We thank Hans-Ude Nissen for his improvements of the manuscript.