Effects of a Rogue Star on Earth’s Climate
A DISSERTATION
SUBMITTED TO THE FACULTY OF THE GRADUATE SCHOOL
OF THE UNIVERSITY OF MINNESOTA
BY
Harini Chandramouli
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF
Doctor of Philosophy
Richard McGehee
August 2020
c© Harini Chandramouli 2020
ALL RIGHTS RESERVED
Acknowledgements
I would like to express my gratitude to the following people without whom I would not
be where I am today:
Dick McGehee, who has been the best advisor I could have asked for. His guidance,
kindness, strength, and unrelenting patience has taught me what it means to be a good
teacher and to never stop fighting for what I want. Bonny Flemming, who is always
willing to lend an ear and making sure that all that paperwork got submitted correctly.
Bryan Mosher, for letting me flex my teaching muscles and all the mentoring you’ve
given me over the years. Rick Moeckel, Clarence Lehman, and Paul Garrett, for serving
on my committee and making my dissertation a smooth process, even in the face of all
that 2020 had to throw at us. Arnd Scheel and Gareth Roberts, for taking the time
to write me letters of recommendation and helping me write and speak about my own
mathematics better. David Swigon, for encouraging me to do mathematical research
and being the first person to ever tell me that I should go to graduate school.
To my many friends who have been with me throughout this journey, I cannot thank
you enough. Sunita Chepuri, for always being by my side (literally). I couldn’t have done
this without you. Cora Brown and Sarah Milstein, for listening to my endless babble
and making me feel like I belonged in the world of mathematics. Alice Nadeau, for
being my graduate school inspiration and pushing me to strive for more. My academic
twin, Somyi Baek, for spending this last year helping me to stay motivated. Daniel
i
Hess, for giving me confidence when I had none and comfort when I needed it the most.
Jimmy Broomfield, Greg Michel, and Mike Loper, for helping to create best cohort that
the School of Mathematics has ever seen. My office mates past and present (Sunita
Chepuri, Cora Brown, Becky Patrias, Carolynn Johnson, Olivia Cannon, Emily Tibor,
and Dario Valdebenito), who have made Vincent Hall 552 to most enviable office in
the whole department. Therese Broomfield, Rachel Johnson, Maclane Marti, and Brian
Barthel, who have helped make Minnesota my home.
Finally, to my family whose support has seen me through numerous crises. Amma,
Appa, Anna, and Kendra, for always making time for me and helping me to keep
everything in perspective. Anand and Veera, for always providing me with the best
baby content to brighten up my many difficult days.
ii
Dedication
For giving me opportunities they never had and for their unwavering support, this work
is dedicated to my parents, Kodaganallur Shankaranarayanan Chandramouli and Latha
Chandramouli.
iii
Abstract
The details of the way in which the Earth orbits the Sun can have profound effects
on Earth’s climate. Elements such as the Earth’s tilt or how tight the orbit is can
affect temperature distribution or glacial formation. One event that could lead to such
changes is if a star passes near our solar system close enough to disturb Earth’s orbit.
Over the Earth’s 4.5 billion year history, many stars could have approached the Solar
System. Disturbances that are forced due to such an event could have effects on the
orbit of planets in the Solar System. Even small changes to the orbits of the planets
can have large effects on the climate a planet experiences as they could disrupt the
Milankovitch cycles. These effects would be visible in the Earth’s sediment core data,
which can be used to help model the Earth’s position as it orbits the Sun.
This work focuses on modeling the effects of a star passing within the Solar System
on the eccentricity, semi-major axis, and argument of the periapsis. The changes to
the Earth’s climate due to the changes on those orbital elements are also considered.
Here the focus is on what the system would look like if the star were to pass within
Kuiper Belt and the Oort Cloud. Passage within the Kuiper Belt has can change the
equilibrium temperature of the Earth up to 2◦ C without even taking into consideration
the interactions of the planets with each other on the climate. The position of a planet
as the star passes has an effect on the system, giving different results for different initial
conditions. This position is seen to lead to dramatic differences between the orbital
elements for various solutions.
iv
Contents
Acknowledgements i
Dedication iii
Abstract iv
List of Tables viii
List of Figures xi
1 Introduction 1
2 The Hyperbolic Restricted Three-Body Problem 5
2.1 Conserved Quantities, Symmetries, and Integrals of Motion . . . . . . . 8
2.2 Relative Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Regularization Methods . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Other Work on the HR3BP . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.1 Approach by Faintich . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.2 Approach by Sorokovich . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
v
3 Energy Balance Models 20
3.1 Description of the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.1 Absorbed Solar Radiation . . . . . . . . . . . . . . . . . . . . . . 22
3.1.2 Outgoing Radiation . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.1.3 Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.4 Dynamic Ice Line . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.5 Values of Parameters in the Budyko Model . . . . . . . . . . . . 26
3.2 Effects of the Semi-Major Axis and Eccentricity on Temperature Distri-
bution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4 Passage in the Kuiper Belt 30
4.1 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.1.1 Values of Parameters in Simulations . . . . . . . . . . . . . . . . 35
4.2 Effects on the Orbit of Jupiter . . . . . . . . . . . . . . . . . . . . . . . 35
4.3 Effects on the Orbit and Climate of Earth . . . . . . . . . . . . . . . . . 43
4.3.1 Global Mean Temperature Changes . . . . . . . . . . . . . . . . 49
4.3.2 Equilibrium Temperature Profile Changes . . . . . . . . . . . . . 50
4.3.3 Bifurcation Diagram in the Semi-Major Axis . . . . . . . . . . . 53
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5 Passage near the Oort Cloud 57
5.1 Varying Star’s Distance with Constant Initial Position of the Planet . . 58
5.1.1 Global Mean Temperature Changes . . . . . . . . . . . . . . . . 62
5.1.2 Equilibrium Temperature Profile Changes . . . . . . . . . . . . . 65
5.2 Varying the Initial Position of the Planet . . . . . . . . . . . . . . . . . 65
5.2.1 Global Mean Temperature Changes . . . . . . . . . . . . . . . . 66
vi
5.2.2 Equilibrium Temperature Profile Changes . . . . . . . . . . . . . 70
5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6 Discussion 73
6.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
References 76
Appendix A. Mathematica Code 84
A.1 Regularized HR3BP in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 84
A.2 Budyko-Widiasih Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
vii
List of Tables
3.1 Parameter values used in the standard Budyko-Widiasih EBM for Earth. 26
4.1 Masses of known stars that have passed or will pass within 5 light-years
of the Sun within ±3 million years measured in M, where 1 M is the
mass of the Sun. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2 Initial conditions of the HR3BP for both Jupiter and Earth as the third
body in an elliptic orbit and the star passing in the middle of the Kuiper
Belt. The system has been scaled so that Jupiter could be a distance of
1 unit away. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 Final value of eccentricity (eJ(smax)), percentage change in the semi-
major axis (% ∆aJ), final value of orbital period (TJ(smax)), and the
final value of the argument of the periapsis (ωJ(smax)) for the simulations
with initial conditions found in the first two columns, along with those
found in table 4.2. Recall that the initial value of eJ is 0.0489, the initial
value for aJ is 1, the initial value of TJ is 7.82391, and the initial value
of ωJ is 0 degrees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
viii
4.4 Final value of eccentricity (e(smax)), percentage change in the semi-major
axis (% ∆aE), final value of orbital period (TE(smax)), and the final
value of the argument of the periapsis (ωE(smax)) for the simulations
with initial conditions found in the first two columns, along with those
found in table 4.2. Recall that the initial value of eE is 0.0167, the initial
value for aE is 0.192, the initial value of TE is 0.658227, and the initial
value of ωE is 0 degrees. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.1 Initial conditions of the HR3BP for Jupiter as the third body in an elliptic
orbit starting on the horizontal axis. The system has been scaled so that
Jupiter could be a distance of 1 unit away. . . . . . . . . . . . . . . . . . 58
5.2 Final value of eccentricity (eJ(smax)), percentage change in the semi-
major axis (% ∆aJ), final value of orbital period (TJ(smax)), and final
value of the argument of the periapsis (ωJ(smax)) for the simulations
with distance of closest approach given by the first column, a, along with
those found in table 5.1. Recall that the initial value of eJ is 0.0489, the
initial value for aJ is 1, the initial value of TJ is 7.82391, and the initial
value of ωJ is 0 degrees. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3 The difference between the global mean temperature (GMT) before the
star passes to after the star passes at various ice line locations and for var-
ious closest approaches of the star. If the number is negative, the global
mean temperature increases after the star passes, and if it is positive, the
global mean temperature decreases after the star passes. . . . . . . . . 64
5.4 Orbital elements after the star passes from varying distances (closest
approach of star given in column 1). The starting locations (column 2)
that resulted in the largest changes in the eccentricity are given. . . . . 67
ix
5.5 Orbital elements after the star passes from varying distances (closest
approach of star given in column 1). The starting locations (column 2)
that resulted in the largest changes in the semi-major axis are given. . . 67
5.6 Orbital elements after the star passes from varying distances (closest
approach of star given in column 1). The starting locations (column 2)
that resulted in the largest changes in the argument of the periapsis are
given. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
x
List of Figures
2.1 Initial two-body problem that is present before the rogue star passes
through the system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Relative coordinate system transformation. . . . . . . . . . . . . . . . . 10
3.1 Sensitivity of changes in the eccentricity to changes in the semi-major
axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Effects of eccentricity on the equilibrium temperature profile for Earth. 28
4.1 The initial set-up of the problem in this chapter 4 is drawn. . . . . . . 36
4.2 Jupiter’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes for initial position (1.0489, 0) of Jupiter. . . . . . . . . . . . 38
4.3 Jupiter’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the periapsis values as the star passes for initial position
(1.0489, 0) of Jupiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.4 Jupiter’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes for initial position (−0.87498,−0.382226) of Jupiter. . . . . 41
4.5 Jupiter’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the periapsis values as the star passes for initial position
(−0.87498,−0.382226) of Jupiter. . . . . . . . . . . . . . . . . . . . . . 41
4.6 Jupiter’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes for initial position (−0.658207,−0.706261) of Jupiter. . . . 42
xi
4.7 Jupiter’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the periapsis values as the star passes for initial position
(−0.658207,−0.706261) of Jupiter. . . . . . . . . . . . . . . . . . . . . 42
4.8 Earth’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes for initial position (0.195206, 0) of Earth. . . . . . . . . . . 45
4.9 Earth’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the perihelion values as the star passes for initial position
(0.195206, 0) of Earth. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.10 Earth’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes for initial position (0.180591,−0.073465) of Earth. . . . . . 47
4.11 Earth’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the perihelion values as the star passes for initial position
(0.138971,−0.1357460) of Earth. . . . . . . . . . . . . . . . . . . . . . . 47
4.12 Earth’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes for initial position (0.195206, 0) of Earth. . . . . . . . . . . 48
4.13 Earth’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the perihelion values as the star passes for initial position
(−0.0702688,−0.17736) of Earth. . . . . . . . . . . . . . . . . . . . . . 48
4.14 Change in global mean temperature for Earth for a few values of eccen-
tricity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.15 Change in global mean temperature for Earth for a few values of the
semi-major axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.16 Equilibrium temperature profile for Earth for the current and perturbed
parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.17 Bifurcation diagram of the semi-major axis aE . All other parameters are
in table 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
xii
5.1 Jupiter’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes at a distance of 100 AUs. . . . . . . . . . . . . . . . . . . . 60
5.2 Jupiter’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the perihelion values as the star passes with its closest ap-
proach at 100 AUs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3 Jupiter’s initial orbit and final orbit are drawn in the (q1, q2) plane as the
star passes at a distance of 200 AUs. . . . . . . . . . . . . . . . . . . . 61
5.4 Jupiter’s initial orbit’s and final orbit’s eccentricity, semi-major axis, and
argument of the perihelion values as the star passes with its closest ap-
proach at 200 AUs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.5 Change in global mean temperature (GMT) based on the changes found
in table 5.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.6 Equilibrium temperature profile for Earth for the current and perturbed
parameters when star passes by at 100 AUs away. . . . . . . . . . . . . 65
5.7 Change in global mean temperature (GMT) based on the changes found
in table 5.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.8 Change in global mean temperature (GMT) based on the changes found
in table 5.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.9 Equilibrium Temperature Profile based on the changes found in table 5.4
and table 5.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
xiii
Chapter 1
Introduction
The orbits of planets in the Solar System depend upon both the Sun and the positions
of the other planets in the system. Disturbances to the orbit of any one of these planets
can have an effect on the rest of the Solar System. Gravitational perturbations can be
caused by moving objects in space, such as rogue stars. A rogue star is defined as either
stars that are free from the constraints of a galaxy (intergalactic rogue stars) or a star
that does not follow normal behavior (independent moving stars). A rogue star which
came close to the Solar System with high enough mass would be able to change the
orbits of the planets and how they interact with each other.
To study if such a star has passed by, scientists and mathematicians look to models
of the Solar System. Knowing how the planets move now and their general patterns,
the motion of the planets in the past are predicted and how these motions change could
indicate the passage of a rogue star. However, the validity of solutions of Solar System
gravitational models are constrained to about 0 - 60 million years ago. While relative
positions of the major planets remain the same over hundreds of millions of years, details
such as the relative phases along orbits can be chaotic, where initially close solutions
diverge exponentially [29, 30]. This chaos can make it difficult for scientists to test
1
2physical theories. Although progress on this front has been made in recent years, it is
still not possible to retrace the precise history of the Solar System using only current
information. On the other hand, Olsen et. al. found a way to use geological data to
constrain the astronomical solutions back in time [44].
Indeed, geological data from the last 60 million years agree with astronomical solu-
tions but beyond that there is little information on the Solar System [31, 45]. Olsen et.
al. established a method to find empirical data beyond 60 million years ago that uses a
network of highly resolved data from multiple temporally correlative and complemen-
tary records they call “the Geological Orrery.” This data comes from cores collected
from the bottoms of lakes and rivers which can be dated using uranium-lead dating
applied to the mineral zircon. The sediments are between 223 and 199 million years
old, in the Mesozoic era. The depths of the water were determined by wet and dry ages
and Olson et. al. were able to conclude that the orbital changes were responsible for
the patterns found in the sediments.
The orbital planes of planets are affected by the gravitational forces of the other
bodies in the Solar System in a quasiperiodic manner which can be broken up into a
series of secular fundamental frequencies that represent each planet’s contribution to
the changes of the orbits. These motions can be described in terms of the precession of
the periapsis, gi frequencies, in the orbital plane and the precession of the orbital plane
itself in space (represented by the precession of the ascending node), si frequencies.
Differences in the gi values yields the eccentricity cycles of 100,000 and 405,00 years
on Earth while the sums of the gi frequencies along with Earth’s axial precession con-
stant, p, yields the climate precession cycle of about 23,000 years. The sum of the si
frequencies along with p yields the obliquity periodicity of about 41,000 years. These
three cycles, obliquity, eccentricity, and precession, are the Milankovitch cycles which
inform the oscillatory nature of Earth’s climate. Olsen et. al. reported a 405,000 year
3cycle in the sediments, known as the McLaughlin cycle, which matches the eccentric-
ity cycle [44]. Unlike the McLaughlin cycle, most of the other cycles were of different
lengths in the Triassic period compared to today. The McLaughlin cycle acts as a mea-
suring stick, for example indicating the length of the precession cycle [44]. The relative
phases of these cycles can be ascribed to positions of the inner planets and Jupiter [44].
This allows for mathematicians and scientists to create conditions in the past that
the solutions of the orbits of the Solar System must satisfy. As the solutions diverge
exponentially, some of the orbits can be ruled out if they do not satisfy the conditions
set forth from this coring data [44].
If a star passed within our Solar System at a range that were to affect the Earth’s
orbit, this would appear in the work done in [44]. The changes spotted in the coring data
would then be related back to the secular fundamental frequencies and would indicate
the planetary orbits going back several million years. In this paper, the types of changes
that would be observed in the orbits of the planets that would occur from a passing star
are studied. As changes to elements of the orbit of a planet have effects on the climate
of that planet, possible climatic effects are also considered in this work.
Broadly, the point of this work is to study the changes that can be observed by
gravitational disturbances to orbits in the Solar System. The focus is placed upon how
the eccentricity, semi-major axis, and argument of the periapsis change as a rogue star
passes at varying distances.
The first part of this work sets up the equations used to model the Sun-planet-
star system that will be focused upon. This modeling is done using a special case of
the gravitational n-body problem and is done throughout the paper in two-dimensions.
Then the framework for the climate modeling is also given, with the Budyko-Widiasih
equation explained along with its relationship to orbit of planets [6, 35, 56].
In chapter 4, the effects of the rogue star passing the Kuiper Belt are studied. The
4effects of this event are calculated both on the orbit of Earth and Jupiter. Jupiter is
considered as it is the most massive of the planets in the Solar System and its orbit
sets the ecliptic. Changes to Jupiter’s orbit could change the secular fundamental
frequencies and have effects on Earth’s orbit. Additionally, Earth’s orbital changes are
also considered along with the climatic change that would result from these changes.
Here it is shown that such a close star passage almost certainly did not occur as the
changes to the eccentricity cycle on Earth would be dramatic and noticeable. As such
an occurrence has not been spotted already, it’s safe to conclude an encounter this close
did not happen.
In chapter 5, the sensitivity of the equations to the distance of the rogue star’s
passage is analyzed. The star’s orbit is moved further and further away from the Solar
System to understand where it is that subtle changes can still be observed. As the star
moves further away, its effects on the Sun-planet system are diminishing, as expected.
The star’s effects become so subtle on the orbit of the planet that the changes become
negligible. Indeed, the effects still exist but are so small that one could attribute the
variation observed to numerical error.
Chapter 2
The Hyperbolic Restricted
Three-Body Problem
An orbit is the gravitationally curved trajectory of an object, such as the trajectory of
a planet around a star or a natural satellite around a planet. Normally, orbit refers to a
regularly repeating trajectory, although it may also refer to a non-repeating trajectory.
To a close approximation, planets and satellites follow elliptic orbits, with the center
of mass being orbited at a focal point of the ellipse, as described by Kepler’s laws of
planetary motion [25]. For most situations, orbital motion is well approximated by
Newtonian mechanics, where gravity is as a force obeying an inverse-square law [26].
Before a star comes by to disturb the solar system, the planets are orbiting the
Sun. Focusing on the Sun and the Earth, there are two bodies’ motions to consider (see
figure 2.1). Note that this planet does not necessarily need to be Earth, but can be
any planet that orbits the Sun. In fact, in this paper, both Jupiter’s and Earth’s orbit
around the Sun will be studied.
The Kepler problem is the case in which two bodies interact by a central force that
varies in strength as the inverse square of the distance r between them. The problem is
5
6Figure 2.1: Initial two-body problem that is present before the rogue star passes through
the system. Here the Earth (green) is orbiting the Sun (yellow) on an orbit of eccentricity
0.0167 [60].
to find the position or speed of the two bodies over time given their masses, positions,
and velocities. Kepler motion is governed by the differential equation
x¨(t) = −µx(t)|x(t)| (2.1)
where x ∈ Rn denotes the position of the particle, the dot is the derivative with respect
to time t, and µ is the gravitational parameter. Using classical mechanics, the solution
can be expressed as a Kepler orbit, the motion of one body relative to another as a
circle, ellipse, parabola, or hyperbola, which forms a two-dimensional orbital plane in
three-dimensional space. Kepler orbits can be parametrized using six orbital elements
– eccentricity (e), semi-major axis (a), inclination (ι), longitude of the ascending node
(), argument of the periapsis (ω), and the true anomaly (ν).
The Kepler problem has a known solution, which is given in polar coordinates (r,
7ν) to be
r =
|Ω|2/µ
1 + e cos(ν)
(2.2)
where Ω represents the angular momentum [46]. Based on the value of e, the solution
is either a circle, ellipse, parabola, or hyperbola.
As the star approaches the system, three bodies are taken into consideration – the
Sun, the planet that orbits the Sun (e.g. Earth), and the passing star. Let the position
of mass mi be given by xi ∈ Rn. The equations that define the motion of the three
bodies are given below.
x¨1 = −Gm2 x1 − x2|x1 − x2|3
−Gm3 x1 − x3|x1 − x3|3
x¨2 = −Gm1 x2 − x1|x2 − x1|3
−Gm3 x2 − x3|x2 − x3|3
x¨3 = −Gm1 x3 − x1|x3 − x1|3
−Gm2 x3 − x2|x3 − x2|3
,
(2.3)
The three-body problem does have a global analytical solution in the form of a conver-
gent power series, but does not have a general analytic solution by algebraic expressions
and integrals. This power series solution as was proven by Sundman [52]. However,
the Sundman series converges slowly; therefore, it is currently necessary to approximate
solutions by numerical analysis.
In order to study the motion of the elliptic orbit of the planet when it is perturbed
by a passing star, one can reduce the problem from the complexities of a general three
body problem (equation (2.3)) with the use of the following assumptions:
1. Two finite masses, m1 and m2, have hyperbolic orbits with respect to their center
of mass, and hence, each has a hyperbolic orbit with respect to the other. These
are referred to as the primary masses.
82. A third body of infinitesimal mass is in the system which does not affect the motion
or location of the center of mass of the two primary masses. In the case studied
here, the Earth’s mass is 333,000 smaller than that of the Sun and Jupiter’s mass
is 1,047 smaller than the Sun [58, 59], allowing for such an assumption to make
sense.
The first assumption defines the Hyperbolic Three-Body Problem. The second assump-
tion defines the Restricted Three-Body Problem, where the mass of a much smaller
object is approximated to be zero. This allows the other two bodies in the system (the
primaries) to create a two-body problem, which is solvable. The third body then moves
in the vector field created by the interaction of the primaries. Together, this gives us
the Hyperbolic Restricted Three-Body Problem (HR3BP).
In this chapter, the HR3BP is set-up in a relative coordinate system that to mirror
the Sun-planet-star system considered in this paper more closely. Then, regularization
methods are applied to create more stability in the numerical integration of the system.
Lastly, other literature on the HR3BP is presented and the approaches are compared
to the one used in this paper.
2.1 Conserved Quantities, Symmetries, and Integrals of
Motion
A Hamiltonian system is a dynamical system that is completely described by a scalar
function H(x,y, t), known as the Hamiltonian, where x,y ∈ R2. The differential equa-
tions for the system are given by the Hamiltonian as follows:
x˙ =
∂H
∂y
(x,y), y˙ = −∂H
∂x
(x,y). (2.4)
9Here, (x(t),y(t)) is the solution of the IVP defined by equation (2.4), known as Hamil-
ton’s equations, and an initial condition (x(t0),y(t0)) = (x0,y0). For equation (2.3),
the Hamiltonian is given to be
H = − Gm1m2|x1 − x2| −
Gm2m3
|x3 − x2| −
Gm3m1
|x3 − x1| +
y21
2m1
+
y22
2m2
+
y23
2m3
, (2.5)
where x˙ = y and the second order equation (2.3) is converted to a system of first order
equations. Here, H does not depend explicitly on time, and as such the Hamiltonian is a
constant of motion, equaling the total energy of the system (gravitational plus kinetic).
Symmetries in the three-body problem yield global integrals of motion that simplify
the problem [39]. Translation symmetry of the problem results in the center of mass
C =
m1x1 +m2x2 +m3x3
m1 +m2 +m3
moving with constant velocity, so that C = L0t+C0, where L0 is the linear velocity and
C0 is the initial position. The constants L0 and C0 represent four integrals of motion.
Rotational symmetry results in the total angular momentum being constant
Ω = x1 × y1 + x2 × y2 + x3 × y3,
where × is the two-dimensional cross product. The two components of the total angular
momentum Ω yield two more constants of the motion.
The integrals C0, L0, and Ω provide six constants of motion, and the Hamiltonian,
H, is the seventh.
10
Figure 2.2: Coordinate system created relative to the Sun (whose mass is given by m1).
The Sun is located at the focus of the orbit of the planet. The location of the Sun is
designated as the origin of the coordinate system.
2.2 Relative Coordinate System
The problem considered in this paper is done in the plane, that is, x1,x2,x3 ∈ R2. In
this system, let m1 represent the Sun, m2 represent the passing star, and m3 represents
the planet. We are interested in solutions where the planet stays close to the Sun, and
as such we look at a coordinate system relative to the Sun (see figure 2.2).
Let
Q = x2 − x1 and q = x3 − x1,
where Q, q ∈ R2. Taking into account that m3 = 0, and choosing the units of mass
11
properly so G = 1, the equations of motion can be derived from equation (2.3).
Q¨ = −(m1 +m2)Q|Q|3 (2.6)
q¨ = −m1q|q|3 −
m2(q −Q)
|q −Q|3 −
m2Q
|Q|3 (2.7)
Equation (2.6) is the Kepler problem (see equation (2.1)), and as such, there is a closed
form solution in the form of polar coordinates (see Equation equation (2.2)). This
equation describes the relative motion of the two stars, which is set to be a hyperbola.
Such a solution can be be characterized by two parameters a > 0 and e > 1. There
is no simple formula for the solution Q(t), but defining both Q and t in terms of the
hyperbolic anomaly, E, the solution can be parametrized by the following equations
[22].
Q(t) = R
cos(ν)
sin(ν)
 (2.8)
R(E) = a (e cosh(E)− 1) (2.9)
tan
(ν
2
)
=
√
e+ 1
e− 1 tanh
(
E
2
)
(2.10)
e sinh(E)− E =
√
m1 +m2
a3
t (2.11)
Equation (2.10) can be rewritten with the substitution
τ = tan
(ν
2
)
,
then
cos(ν) =
1− τ2
1 + τ2
and sin(ν) =
2τ
1 + τ2
.
Thus, equation (2.8) can be written as a function parametrized by the hyperbolic
12
anomaly, E, and is the solution of equation (2.6).
This solution can then used to solve equation (2.7) to determine the motion of
Earth’s orbit as the rogue star passes through the system. Note that equation (2.7) is
very similar to the Kepler problem, equation (2.1), but with added term F (q, t) where
F (q, t) = −m2(q −Q)|q −Q|3 −
m2Q
|Q|3 .
2.2.1 Regularization Methods
Regularization is a transformation of both space and time variables, introduced by
Levi-Civita in the plane and generalized by Kustaanheimo and Stiefel in space [28, 33].
Historically, regularization was developed for investigating singularities of Kepler motion
and for describing collisions of two point masses as well as for improving the numerical
integration of collision and near-collision orbits.
The Sun and the rogue star are modeled in a manner that ensures that they will not
collide or get very close to each other. However, Earth’s orbit is modeled to be close to
the Sun’s, requiring regularization methods to aid the accuracy of numerical solutions.
Equation (2.7) models the Earth’s motion, relative to the Sun. Let q˙ = p, where p is
the momentum. Then equation (2.7) can be converted to a system of first order ODEs.
q˙ = p (2.12)
p˙ = −m1q|q|3 −
m2(q −Q)
|q −Q|3 −
m2Q
|Q|3 (2.13)
The time-dependent Hamiltonian of the system created by equation (2.12) and equa-
tion (2.13) is
H(q,p, t) = 1
2
|p|2 − m1|q| −
m2
|q −Q| −
m2
|Q| . (2.14)
13
Let the Kepler energy of (q,p) be given by
K(q,p) = 1
2
|p|2 − m1|q| = −h < 0
where h > 0 is a real constant. The value of K(q, p) is necessarily negative since the
planets orbit the Sun in an ellipse before the rogue star passes.
Levi-Civita Regularization
The Levi-Civita regularization introduces a new time variable and creates a conformal
squaring map. Complex notation will be used throughout, i.e. instead of the vector
x =
x1
x2
 ∈ R2, the corresponding complex coordinate x = x1 + ix2 ∈ C will be used.
The Levi-Civita transformation [33, 51] is given by the equations
q = 2z2 (2.15)
p =
w
z
. (2.16)
Additionally, in the place of the physical time t, a new independent time variable τ is
introduced by the differential equation
dt = 2|z|2dτ. (2.17)
Let derivatives with respect to τ be denoted by primes. Then the system of equa-
tion (2.12) and equation (2.13) becomes
z′(τ) =
1
2
w (2.18)
w′(τ) = −hz − 2|z|2z · F (2z2, t(τ)) (2.19)
14
with the Kepler energy of (z, w) given by
|w|2
2|z|2 −
m1
2|z|2 = −h. (2.20)
Note that this is not the Hamiltonian for the non-perturbed problem defined by
z′(τ) =
1
2
w
w′(τ) = −hz,
Instead, we find the value for the new energy surface by looking at the isoenergetic
reduction and find that the Hamiltonian for the system is actually given by [63]
K˜(z, w) = 1
2
|w|2 + h|z|2 = m1
4
. (2.21)
The system of equations defined by equation (2.18) and equation (2.19) are regular-
ized in the sense that the singularity at the origin has been eliminated from the original
system of equation (2.12) and equation (2.13).
Translating equation (2.14) to Levi-Civita coordinates yields the equation
|w|2
2|z|2 −
m1
2|z|2 −
m2
|q −Q(t)| −
m2
|Q(t)| .
The Hamiltonian for the system of equations equations (2.18) and (2.19) is given to be
H˜(z, w, τ) = 1
2
|w|2 + h|z|2 − m2|z|
2
|2z2 −Q (t(τ)) | −
m2
|Q (t(τ)) | (2.22)
In both the system of equations given by equations (2.18) and (2.19) and equation (2.22),
t is written as a function of τ . Given that the relationship between the two that is known
is given by equation (2.17), it will get added to the set of ODEs that will get solved.
15
Lastly, recall that h is no longer a constant, but a function of time as well. Putting all
this together yields the following system of equations
z′ =
w
2
w′ = −hz − 2|z|2z F (2z2, t)
h′ = 2F (2z2, t) · zw
t′ = 2|z|2
. (2.23)
Since h, t ∈ R, this is a six dimensional system.
Recall that equation (2.8) was parametrized by the hyperbolic anomaly, E. As Q
can be thought of as Q(E), The function F (q, t) can instead be thought of as F (q, E),
which becomes easier to evaluate given how Q is parametrized by E. The ODE for E′
can be substituted in the place of t′.
z′ =
w
2
w′ = −hz − 2|z|2z F (2z2, E)
h′ = 2F (2z2, E) · zw
E′ =
2|z|2
e cosh(E)− 1
(2.24)
Equations (2.24) defines the motion of the planet in the vector field created by the
interactions of the Sun and the passing star in the planar situation.
2.3 Other Work on the HR3BP
The problem of finding the general solution of the n-body problem was considered very
important and challenging and can find its roots in the time of Ptolemy. Indeed, looking
at solutions of the HR3BP is not a new problem, and various other approaches exist.
16
There’s a body of work focused on describing the motion of the third body which is
in third body moves on a straight line perpendicular to the line of motion of the two
primaries [8, 9, 10]. Work has also been done with restricting the motion of all three
bodies to the unit circle [34]. Two different studies on the HR3BP are elaborated upon
below as these both specifically consider the case where the third body is in an elliptic
orbit around one of the primaries [21, 49].
2.3.1 Approach by Faintich
Faintich studied the equations of motion for the HR3BP for a Sun-star-comet system
where the equations were studied in a set of non-rotating axes centered on the Sun with
the hyperbolic anomaly as the independent variable [21]. In this set-up, the Sun and
the passing star are in hyperbolic orbits in the plane, and the third body, the comet is
allowed to approach from outside the plane and hence the equations of motion are in
three dimensions (see equations (2.25) to (2.27)). Here, d1 represents the distance to
the star and d2 the distance to the Sun.
d2q1
dE2
− e sinh(E)
(e cosh(E)− 1)
dq1
dE
− m1a
(m1 +m2)
(
cosh(E)− e sinh
2(E)
(e cosh(E)− 1)
)
+
a3(e cosh(E)− 1)2
(m1 +m2)
(
m1(q1 + a(e− cosh(E)))
d31
+
m2q1
d32
)
= 0;
(2.25)
d2q2
dE2
− e sinh(E)
(e cosh(E)− 1))
dq2
dE
+
m1a
√
e2 − 1
(m1 +m2)
(
sinh(E)− e sinh(E) cosh(E)
(e cosh(E)− 1)
)
+
a3(e cosh(E)− 1)2
(m1 +m2)
(
m1(q2 + a
√
e2 − 1 sinh(E)
d31
+
m2q2
d32
)
= 0;
(2.26)
d2q3
dE2
− e sinh(E)
(e cosh(E)− 1)
dq3
dE
+
a3(e cosh(E)− 1)2
(m1 +m2)
(
m1q3
d31
+
m2q3
d32
)
= 0. (2.27)
Faintich’s applies these equations to a hypothetical star-sun-comet system to deter-
mine the effect of the stellar encounter on the orbit of the comet. The comet is initially
in a near-parabolic orbit around the Sun (e = 0.999980) starting around 50,000 AU
17
away fro the Sun. The passing star is considered to affect the Sun-comet system for a
finite time period.
Various numerical examples are integrated from ν = −1.5 to ν = 1.5 (see table 1
in [21]). Some initial conditions show that the stellar passage changed the path of the
comet; if unperturbed, the comet could have come into the Solar System so far as to be
detected from Earth. However, the comet’s interaction with the passing star changed
its trajectory and it would remain undetected.
2.3.2 Approach by Sorokovich
Sorokovich analyzed the single-averaged HR3BP in the plane [49]. The first body and
second body are placed in hyperbolic orbits with respect to each other, but the coordi-
nate system used places the first body at the origin (similar to the Sun in the system
studied in this paper). The third body with very small mass is assumed close to the
first body (similar to the Earth in the system studied in this paper). The perturbation
function of the problem is given by the equation
F = m2
(
1√
q2(t) +Q2(t)− 2q(t)Q(t) cos(α) −
q(t) cos(α)
Q2(t)
)
(2.28)
where α is the angle between q(t) and Q(t) and is the sum of the true anomalies of the
orbits of the second body (ν˜) and the third body (ν) and the argument of the periapsis
of the third body (ω) [18]. Assuming that q(t) Q(t), [49] expands equation (2.28) in
a series in powers of q(t)Q(t) . Discard from this expansion the term not dependent on the
coordinate of the perturbing body, the second body, the expansion is
F =
m2
Q(t)
∞∑
k=2
(
q(t)
Q(t)
)k
Pk (cos(α)) , (2.29)
18
where Pk (cos(α)) is a Legendre polynomial of k
th order. Averaging F over the mean
anomaly M yields
F =
1
2pi
2pi∫
0
F dM =
m2
a2
√
1− e2
∞∑
k=2
1
Qk+1(t)
1
2pi
2pi∫
0
qk+2(t)Pk (cos(α)) dν (2.30)
where
dM =
q2(t)dν
a2
√
1− e2 , (2.31)
where a and e are the semimajor axis and eccentricity of the third body. Assuming that
the fourth order terms are negligible due to the fact that q(t)  Q(t), only the first
order term is retained and the equation becomes
R =
m2
a2
√
1− e2Q3(t)
1
2pi
2pi∫
0
q4(t)P2 (cos(α)) dν. (2.32)
Using a series expansion on q(t), assuming it is in the form of an ellipse, and Q(t),
assuming it forms a hyperbola, the value for P2(x), and changing to ν˜ as the independent
variable instead of t, ODEs can be derived to study the motion of the third body [1, 49].
de
dν˜
=
15
4
βe
√
1− e2 (1 + e˜ cos(ν˜)) sin (2(ω − ν˜)) (2.33)
dω
dν˜
=
3
4
β
√
1− e2 (1 + e˜ cos(ν˜)) (1 + 5 cos (2(ω − ν˜))) (2.34)
where
β =
(
m2
m1
)1/2(a
˜`
)3/2
with ˜` the semi-latus rectum and e˜ is eccentricity of of the orbit of the second body.
Sorokovich’s results suggest that the stability of the unperturbed motion holds when
19
the orbital eccentricity of the third body is small and it remains close to the first body.
The exact and singly averaged equations (equation (2.33)) were integrated numerically
and it was found that while the eccentricity of the orbit of the third body changes as
the second body passes close by, it returns back to its original state eventually. More
interestingly, the argument of the periapsis changes by a small amount as the second
body passes and it does not return to its original state.
2.4 Discussion
In this chapter, the equations of motion for the Sun-planet-star system that is to be
studied were derived, see equation (2.6) and equation (2.7). These equations describe
the relative motion with respect to a non-rotating coordinate system where the Sun is
set at the origin in the plane, and at the focus of the orbit of the planet. Furthermore,
regularization techniques were applied to improve numerical integration techniques go-
ing forward (discussed in chapters 4 and 5), see equation (2.24).
Two other studies of the HR3BP were also considered in this chapter. In Faintich’s
work, the third body is in a nearly parabolic orbit. This eccentricity is much higher than
the eccentricities of orbits considered in this paper. In fact, the largest eccentricities of
the orbit of planets in the Solar System are just a little above 0.2 [61]. Faintich also
allows the comet to move in three dimensions, whereas the third body in the equations of
this paper are restricted to movement in the plane (see equations (2.12) and (2.13)) and
their equations are not in Hamiltonian form. In Sorokovich’s work, the set-up is similar
to the Sun-planet-star system of this paper but the perturbing effect of the equivalent
of the passing star is averaged over the mean anomaly of the orbit. Averaging over
the orbit allows for more analytical results to be proven about how the motion could
change, but the exact numerical values where events occur along the orbit are lost. As
this paper is is focused on these numerical values, averaging is not done.
Chapter 3
Energy Balance Models
Energy balance climate models, often referred to as energy balance models (EBMs),
gained popularity in the late 1960s with the publications of papers by Mikhail Budyko
and William Sellers [6, 47]. EBMs are low-dimensional models that attempt to model
Earth’s surface temperature distribution by the constraint of energy conservation with
a few parameters such as the solar radiation and planetary albedo (reflectivity) [6, 47,
55, 56]. As high-speed computers have become more easily accessible, larger scale Gen-
eral Circulation Models (GCMs) have become more prominent in the study of Earth’s
climate. While GCMs try to take into consideration the incredible complexity of the
real climate system, an incredible amount of data is required to constrain all the phe-
nomenological coefficients and it is not always clear how sensitive the results are to the
parameters or variables [43].
EBMs continue to be vital tools in the understanding of Earth’s climate. These
models describe global behaviors of the climate system, whereas GCMs give more local
results. For example, EBMs give results about general trends of sea-ice from the phe-
nomenon of ice-albedo feedback, but will not give results on what locations in the arctic
will have sea-ice [56]. The low-dimensional characteristic of EBMs make them easy to
20
21
modify and run to equilibrium, making them good choices for the study of paleoclima-
tology. Early research using EBMs was concerned with finding the causes of ice ages and
to investigate climate forcings over the past few hundred years [11, 13, 14, 15, 42, 48].
EBMs have also been useful in modeling the distant past [12, 16, 41].
In this chapter, the Budyko-Widiasih EBM is described along with its assumptions.
The effects of two specific orbital elements, the semi-major axis and the eccentricity, on
the equilibrium temperature profile generated from this EBM are discussed.
3.1 Description of the Model
The energy balance equation describes the evolution of the annual average temperature
T (y, t), where t denotes the time in years and y denotes the sine of the latitude (as
explained by North, the latitude θ is more easily represented by y as y is directly
proportional to surface area [40]). Here, symmetry about the equator is assumed and
so 0 ≤ θ ≤ pi2 or 0 ≤ y ≤ 1. Let β represent the obliquity. The temperature evolves
based on the equation
R
∂T
∂t
= Qs6(y, β) (1− α(y, η))− (A+BT (y, t))− C
(
T (y, t)− T (t)) (3.1)
where the change in temperature is determined by the absorbed solar radiation given
by the term Qs6(y, β) (α(y, η)), the outgoing radiation (A+BT (y, t)), and the trans-
portation of heat between the latitudes C
(
T (y, t)− T (t)). This continuous adaptation
of Budyko’s original difference model was done by Tung [6, 54]. Futher, McGehee and
Lehman combined the Budyko model version formulated by Tung with Earth’s orbital
parameters and compared it with the climate data [35]. In addition to equation (3.1),
22
Widiasih’s dynamic ice-line equation
dη
dt
= ε (T (η, t)− Tc) (3.2)
is considered [56]. The movement of the ice line is determined by the temperature at
the ice line relative to a critical temperature Tc, where ice forms
3.1.1 Absorbed Solar Radiation
The incoming energy at any point on a planet depends upon the location of the point,
the planet’s orbital parameters (e.g. obliquity), the position of the planet along its
orbit, and the solar energy output. Using Kepler’s laws and integrating over an entire
year, it can be shown that the global annual average power flux is given by
Q(a, e) =
Ka2√
1− e2 (3.3)
where K is proportional to the solar output, a is the semi-major axis of the planet’s
orbit, and e is the eccentricity [35]. Q represents the amount of energy received from
the Sun at the top of the atmosphere.
The function s(y) is the latitudinal distribution of that energy, which can be com-
puted from the Earth’s orbital elements. For planet like Earth where all longitudes see
the Sun with equal frequency, the annual average insolation distribution reduces to a
function only of obliquity β and latitude θ
s(θ, β) =
2
pi2
∫ 2pi
0
√
2− (cos(θ) sin(β) sin(γ)− sin(θ) cos(β))2 dγ (3.4)
where γ is the longitude [35]. As latitude is measured from the Equator, recall that
0 ≤ γ ≤ pi/2. Obliquity is the angle between the angular momentum vector of the
23
planetary orbit and the angular momentum vector of the planetary spin, so 0 ≤ β ≤ pi.
Although equation (3.4) gives the insolation, using the integral formula can be quite
cumbersome and time consuming. To bypass the issue, North explicitly gives a second-
degree approximation for the insolation distribution for the Earth as
s(y) = 1− 0.482P2(y)
stating that the approximation to the actual distribution is accurate to within 2% of
the integral formulation [40]. Here P2(y) represents the second Legendre polynomial.
North notes that this approximation was first given by Chylek and Coakley as a linear
interpolation of the insolation distribution, although no closed-form formula is given in
that paper [7]. It should be noted that this quadratic approximation has been used
widely. This approximation was computed only for the current obliquity of Earth, and
so it cannot be used to compute changes due to the Milankovitch cycles. Nadeau and
McGehee’s sixth degree approximation,
s6(y, β) = 1− 5
8
P2 (cos(β))P2(y)− 9
64
P4 (cos(β))P4(y)− 65
1024
P6 (cos(β))P6(y),
(3.5)
which is dependent upon obliquity β, will be used instead [38].
The fraction of the radiative energy absorbed by the planet at point y, given that
the ice line is at η, is represented by the factor
1− α(y, η).
For each η, α(y, η) is a real valued function over the y-interval. Assume that the surface
is either water (ice-free) or ice-covered and that there is only one ice-line η. Note that ice
24
reflects more sunlight than ocean does, and so the surface covered with ice has higher
albedo than that covered with ocean. The albedo function that will be used here is
piecewise constant given by
α(y, η) =

αice, y > η
αice+αH2O
2 y = η
αH2O y < η
. (3.6)
Recall that a symmetry assumption is being made, meaning that whatever happens in
the northern hemisphere also happens in the southern hemisphere. Assume the planet
has an ice cap, such as Earth does. That is, the ice is located above the ice line and
there is water below the ice line.
3.1.2 Outgoing Radiation
The incoming radiative energy described in section 3.1.1 is balanced by the reemission
of energy along with the transport term (section 3.1.3). The “energy out” is given by
the reemission term
− (A+BT (y, t)) (3.7)
which takes into account many phenomena, such as the Stefan-Boltzmann law of black
body radiation and the greenhouse gas effect on the atmosphere. The parameters A
and B used were determined empirically from satellite data [54]. Note that Kaper and
Engler mention in their book that linearizing the Stefan-Boltzmann law overestimates
the values of A and B for Earth gotten from this satellite data [24].
25
3.1.3 Transport
The energy is transported from warmer latitudes to cooler latitudes and this phe-
nomenon is approximated by the term
C
(
T (y, t)− T (t)) (3.8)
where the constant C is a parameter chosen to fit solutions to the current climate of
Earth [54]. This transport term is a simpler linear relaxation to the mean annual global
temperature, where
T (t) =
1∫
0
T (y, t) dy
is the mean annual global temperature. Many phenomena are included in this term
as well, mostly atmospheric and oceanic circulation. While it may seem a simplified
approach, many studies have shown that this approximation reproduces relevant results
when applied to Earth ([36, 56], among many others). Alternatively, one can model
heat transport using a diffusion term as Sellers did in his original work [47]. In this
paper, this possibility is not explored.
3.1.4 Dynamic Ice Line
The main idea behind the ice line equation is as follows: if the average temperature
across the ice line is above the critical temperature, some ice will melt, moving the ice
line toward the pole. If it is below the critical temperature, the ice will advance toward
the equator. Equation (3.2),
dη
dt
= ε (T (η, t)− Tc)
has two important parameters: ε and Tc. Tc, the critical temperature, is often cited as
-10 ◦C [54]. There parameter ε controls how fast the ice line changes in relation to the
26
temperature and therefore must be positive. When ε = 0, the ice lines do not move
when the temperature changes and instead the temperature moves as a result of the
location of the ice. This is not physical. In the limit as ε approaches infinity, The ice
line moves instantaneously to changes in the temperature, but this is also not realistic.
3.1.5 Values of Parameters in the Budyko Model
Parameter Name Value Units
R surface layer heat capacity 12.6 W m−2 K−1
K solar output 342.998 W m−2
a semi-major axis 1.00000011 AU
e eccentricity 0.0167 dimensionless
β obliquity 23.44 degrees
αice albedo of the ice 0.62 dimensionless
αH2O albedo of the water 0.32 dimensionless
A greenhouse gas parameter 202 W m−2
B outgoing radiation 1.9 W m−2 K−1
C heat transport 3.04 W m−2 K−1
ε ice line response to temperature change varies K−1 yr−1
Tc critical temperature -10
◦C
Table 3.1: Parameter values used in the standard Budyko-Widiasih EBM for Earth.
The values commonly used for the Budyko-Widiasih EBM parameters for Earth are
given in table 3.1, along with their units and a brief description. These values can be
found in many locations, such as [35, 54, 56, 60].
3.2 Effects of the Semi-Major Axis and Eccentricity on
Temperature Distribution
According to Laskar et. al, the semi-major axis of Earth’s orbit is essentially constant
[32]. This assumption is prevalent throughout work done using the Budyko model (for
example, in [35, 56], among many others). If a big enough star were to pass close enough
27
to our Solar System, the value of the semi-major axis would be affected. To see the
effect this orbital element would have on the climate of the planet, the relationship
between a and e are studied below in the case of constant global annual average power
flux given in equation (3.3). The climate record indicates that the eccentricity oscillates
between 0 and 0.06 for Earth’s orbit [32]. This oscillation creates sufficient changes in
Earth’s climate that are easily seen in the geological record. Examined here are what
changes in the Earth’s semi-major axis would have a corresponding impact.
Assuming that the global annual average power flux is constant in equation (3.3),
the relationship between small changes in e and small changes in a is considered as to
understand how the variations in one affect the other. To understand this, let
K(a+ ∆a)2√
1− e2 =
Ka2√
1− (e+ ∆e)2 (3.9)
where ∆e represents a small change in eccentricity and ∆a represents a small changes
in semi-major axis. Then it can be observed that
∆e =
√√√√(∆aa + 1)4 + e2 − 1(
∆a
a + 1
)4 − e, (3.10)
Figure 3.1 graphs this relationship using the orbital parameter values for Earth’s orbit.
It is noticeable that small changes to the semi-major axis correspond to stronger changes
to the eccentricity values. On Earth, changing the semi-major axis by approximately
0.00141 would be equivalent to changes in the eccentricity of the orbit by 0.06.
Additionally, eccentricity is one of the three Milankovitch cycles and has a 100,000
and 405,000 year cycle [23]. Changes in the eccentricity of the orbit would lead to
a disruption of these cycles, visible in the climate record. Eccentricity is of prime
importance to glaciation in that it alters the distance the Sun’s short wave radiation
28
0.02 0.04 0.06 0.08 0.10
Δa
0.1
0.2
0.3
0.4
0.5
Δe
Figure 3.1: Sensitivity of changes in the eccentricity to changes in the semi-major axis.
e = 0.0167086 and a = 1.00000011 AU (Earth values).
0.2 0.4 0.6 0.8 1.0
y
-20
20
40
60
80
100
Equilibrium Temperature
Figure 3.2: Equilibrium temperature profiles using the parameters of Budyko’s model
for Earth (see table 3.1) for various e values. Eccentricities drawn are e = 0 (red),
e = 0.0167 (orange), e = 0.1 (magenta), e = 0.25 (green), e = 0.5 (blue), and e = 0.75
(purple). a is held constant at 1 AU. The ice line, η is held constant here at 66.5◦ N.
29
must travel to reach Earth, subsequently reducing or increasing the amount of radiation
received at the Earth’s surface in different seasons.
The equilibrium temperature profile from Budyko’s model is given by
T ∗(η, y) =
(Q(a, e) ∗ s6(y)(1− α(η, y))−A+ C ∗ T ∗(η))
B + C
. (3.11)
A derivation of this equilibrium solution can be found in [56]. The effects on the equi-
librium temperature profile due to changes in the eccentricity values can be observed
in figure 3.2, where the vertical axis is given in degrees Celsius. In figure 3.2, equa-
tion (3.11) is plotted for various values of e and the values of Earth’s orbit were used
for all other parameters (see table 3.1). As the eccentricity increases towards 1, the
equilibrium temperature profile gets warmer.
3.3 Discussion
In this chapter, the Budyko-Widiasih EBM was presented as a way to model a planet’s
temperature distribution. The model’s terms were explained and parameter values
were given as well. The sensitivity of changes in the semi-major axis to changes in the
eccentricity of the orbit were seen (see figure 3.1). Additionally, the effect that changes
in eccentricity would have on the equilibrium temperature profile on Earth were studied.
Assuming the semi-major axis is held constant, equilibrium temperature profiles
appeared to increase as orbits became more eccentric. The relative increase in solar
irradiation at closest approach to the Sun (perihelion) compared to the irradiation
at the furthest distance (aphelion) is slightly larger than four times the eccentricity.
Additionally, when Earth’s orbit becomes more eccentric, if the semi-major axis is held
constant, this necessarily implies that the semi-minor axis shortens. This leads to
increases in the magnitude of seasonal changes [3].
Chapter 4
Passage in the Kuiper Belt
The gravitational influence of Jupiter has helped shape the solar system. The orbits of
most of the system’s planets lie closer to Jupiter’s orbital plane than the Sun’s equatorial
plane (Mercury is the only planet that is closer to the Sun’s equator in orbital tilt)
[59, 61]. In this sense, Jupiter exerts a major influence over the orbital plane of the
solar system. In fact, Jupiter contributes over 60% to the invariable, or mean, plane of
the solar system (the plane perpendicular to the total angular momentum vector of the
solar system that passes through its barycenter) [50]. Earth’s orbit is, thus, dependent
upon the orbit of Jupiter. If a star were to pass near our Solar System to exert influence
upon Earth’s climate, it would need to pass close enough to change Jupiter’s orbit.
Consider the known stars that have passed or will pass within 5 light-years of the Sun
within ±3 million years (see table 4.1) [2, 4, 27]. These stars have a higher probability of
passing near the Solar System. The average mass of these star is equal to 0.550556 M,
or approximately 55% of the mass of the Sun. Thus, this is the value that will be used
for the mass of the passing star in the problem set-up.
30
31
Star Name Mass in M
Gilese 710 0.4-0.6
Scholz’s star and companion brown dwarf A: 0.095, B: 0.063
HD 283856 ∼0.8
TYC 1662-1962-1 ∼0.8
HD 7977 ∼1.2
2MASS J2146+3813 ∼0.15
2MASS J0634-7449 ∼0.6
TYC 2730-1701-1 ∼1
2MASS J0409+0245 ∼0.4
Gliese 3649 0.49
2MASS J0605+4020 ∼0.5
2MASS J1818-4038 ∼0.5
BD-21 1529 ∼0.95
Ross 248 0.136
Proxima Centauri 0.15
Alpha Centauri AB A: 1.100, B: 0.907
Gliese 445 0.15?
HIP 117795 ∼0.5
2MASS J0625-2408 ∼0.5
Barnard’s Star 0.144
2MASS J2241-2759 ∼0.5
Gliese 3379 0.19
Zeta Leporis 2.0
Lalande 21185 0.39
2MASS J1941-4602 ∼0.15
AVERAGE 0.550556
Table 4.1: Masses of known stars that have passed or will pass within 5 light-years of
the Sun within ±3 million years measured in M, where 1 M is the mass of the Sun.
32
Let the force between the Earth and Jupiter be given by
Fd =
Gm1m2
d2
,
and the force between the passing star and Earth be given by
FD =
Gm2M
D2
,
where M is the mass of the passing star (0.550556 M), m1 is the mass of Earth (3.003×
10−6 M), m2 is the mass of Jupiter (9.548× 10−4 M), G is the universal gravitational
constant, d the shortest distance between the Earth and Jupiter (4.206 AU), and D the
shortest distance between the Earth and the passing star. For Fd  FD,
D2  M
m1
d2 =⇒ D  1800.483244 AU.
If a star were to pass within approximately 1800 AU from the Sun, it would exert more
force on the Earth than Jupiter does. Since Jupiter affects Earth’s orbit enough to
change its climate, we might expect the event of a star passing this close to the Sun
to be observable in the geological record. This distance goes far past the Kuiper Belt,
which ranges from approximately 30 AU to 50 AU away from the Sun [17]. The Kuiper
Belt is a circumstellar disc region of icy bodies beyond the orbit of Neptune made up
of objects composed largely of frozen volatiles (e.g. methane, ammonia, water) [37].
Similar to the asteroid belt, the Kuiper Belt is a region of leftover objects from the
solar system’s early history.
In this chapter, both the effects on Jupiter, along with Earth, are considered by the
passage of a star at the edge of the Kuiper Belt. First the initial conditions are laid
out, along with their specific values for simulations. Then, simulation results are given,
33
along with some discussion of these results. For Earth, the effects on the climate are
also considered. The study in this chapter is limited to looking at closer approaches
within the Kuiper Belt. Sensitivity to the distance of the star is studied in chapter 5.
4.1 Initial Conditions
The initial conditions to consider for this problem involve both the hyperbolic orbit of
the passing star and the orbit of Earth and Jupiter around the Sun. The star passes by
on a hyperbolic orbit defined by the length of the semi-major axis a and the eccentricity
e > 1. Equation (2.8) gives this relationship in terms of the hyperbolic anomaly E. The
vertex of the hyperbola is set to be on the x-axis. This will be at 50 AU from the Sun,
the edge of the Kuiper Belt.
The elliptic orbit of the planet, for example Jupiter, can be parametrized with the
focus of this orbit, where the Sun is located, at the origin. The following calculations
can also be done for Earth’s orbit, simply replacing the semi-major axis and eccentricity
values with those of the orbit of Earth. The parametrization is
q1(s) = aJ cos(s) + aJeJ (4.1)
q2(s) = aJ
√
1− e2J sin(s) (4.2)
p1(s) = −
√
m1√
aJ(1 + e cos(s))
sin(s) (4.3)
p2(s) =
√
m1
√
1− e2J√
aJ(1 + eJ cos(s))
cos(s) (4.4)
where
ds
dt
=
m
1/2
1 a
−3/2
J
1 + eJ cos(s)
. (4.5)
34
Now the initial conditions are given simply by knowing the eccentricity and semi-major
axis of the orbit, which are both known values for Jupiter and Earth [60, 58]. Equa-
tion (4.5) comes from the famous Kepler equation and its derivative,
t =
√
a3J
m1
(s+ eJ sin(s)) ,
where s is the eccentric anomaly of the orbit. Recall as well that as Jupiter is on an
elliptic orbit, the Hamiltonian will be a constant, negative number. From chapter 2, let
H˜(q,p, t) = 1
2
|p|2 − m1|q| = −h,
where h =
m1
2a
.
Equations (4.1) to (4.4) suggest that the eccentric anomaly s might be a more
suitable independent variable. In fact, this relationship is what motivates the Levi-
Civita transformation. Using the Levi-Civita transformation given by equations (2.15)
and (2.16) to find the initial conditions in the coordinate system used so far yields the
parametrization below.
z1(τ) =
1√
2
√
aJ(1 + eJ) cos
(√
m1
aJ
τ
2
)
(4.6)
z2(τ) =
1√
2
√
aJ(1− eJ) sin
(√
m1
aJ
τ
2
)
(4.7)
w1(τ) = − 1√
2
√
m1
aJ
√
aJ(1 + eJ) sin
(√
m1
aJ
τ
2
)
(4.8)
w2(τ) =
1√
2
√
m1
aJ
√
aJ(1− eJ) cos
(√
m1
aJ
τ
2
)
, (4.9)
where
s =
1
2
√
m1
aJ
τ.
35
The Hamiltonian for the Levi-Civita coordinates equation (2.21) from the unperturbed
problem gives an extra equation that the initial conditions must satisfy. Namely,
h =
m1 − |w|2
2|z|2
as the orbit starts on the energy surface that this equation satisfies. The semi-major
axis and eccentricity of Jupiter’s orbit are known values and thus the orbit of Jupiter
can be used as the initial condition [58]. An image of the set-up is given in figure 4.1.
Recall again that equations (4.6) to (4.9) are the same initial conditions used for the
problem where Earth is the third body, but replacing aJ with aE and eJ with eE .
4.1.1 Values of Parameters in Simulations
The specific values used for the initial conditions in the simulation appear in table 4.2.
In the table, a brief description of the parameters, along with their values, are given. The
table includes both initial condition values for Jupiter as the third body and Earth as
the third body. All initial conditions associated to Jupiter are noted with a J subscript,
similarly all initial conditions associated to Earth are noted with an E subscript.
4.2 Effects on the Orbit of Jupiter
First, Jupiter is considered as the third body and the effects on its orbit are studied.
Equation (2.24) were numerically solved using Mathematica, specifically the NDSolve
function, utilizing the initial conditions described in section 4.1. All initial conditions
were set at time -smax= 500. This choice was made so that the orbital parameter values
settled towards a constant (with little numerical error). The ODE was solved over the
time period [-smax, smax].
36
-4 -2 2 4 6 8 10
-10
-5
5
10
Figure 4.1: The initial set-up of the problem in this chapter is drawn with Jupiter as
the third body. The path in red is Jupiter’s initial orbit. The dashed blue line details
the path the star will travel as time evolves.
37
Parameter Name Value
a semi-major axis of orbit of passing star 9.6153846154
e eccentricity of orbit of passing star 2
E hyperbolic anomaly of orbit of passing star -3
aJ semi-major axis of orbit of Jupiter 1
eJ eccentricity of orbit of Jupiter 0.0489
TJ period of Jupiter’s orbit 7.82391
ωJ argument of the periapsis of Jupiter’s orbit 0
aE semi-major axis of orbit of Earth 0.192
eE eccentricity of orbit of Earth 0.0167
TE period of Earth’s orbit 0.658227
ωE argument of the periapsis of Earth’s orbit 0
m1 mass of Sun 0.6449299
m2 mass of passing star 0.3550701
Table 4.2: Initial conditions of the HR3BP for both Jupiter and Earth as the third body
in an elliptic orbit and the star passing in the middle of the Kuiper Belt. The system
has been scaled so that Jupiter could be a distance of 1 unit away.
In the simulation, the star take approximately a quarter of an orbit to pass, approx-
imately “three months” time. As the star passes by the solar system, there is a small
but noticeable change in Jupiter’s orbit (see figure 4.2). Here, the dashed, grey curve
indicates the original orbit while the solid, red curve indicates the final orbit. Addi-
tionally, the plot drawn is in the position plane, that is, it is drawn in the (q1, q2)-plane
and not the (z1, z2)-plane. The initial orbit is shifted slightly inward, as if the star had
pushed orbit away.
As the orbit shifts, so do the orbital elements. Specifically, the eccentricity and
semi-major axis of the orbits shift sightly. The final value for the eccentricity of the
orbit is an increase of approximately 0.01 (see figure 4.3a). Additionally, the final value
for the semi-major axis is a decrease of approximately 3.1% (see figure 4.3b).
The change in the semi-major axis results in a change to the orbital period of the
38
-1.0 -0.5 0.5 1.0 q1
-1.0
-0.5
0.5
1.0
q2
Figure 4.2: Jupiter’s initial orbit (grey, dashed) and final orbit (red) are drawn in
the (q1, q2) plane, showing the effects of the passing star. Jupiter’s initial position is
(q1(-smax), q2(-smax)) = (1.0489, 0).
-400 -200 200 400 s
0.02
0.04
0.06
0.08
0.10
eccentricity
(a)
-400 -200 200 400 s0.6
0.8
1.0
1.2
1.4
semimajoraxis
(b)
-400 -200 200 400 s
5
10
15
20
perihelion angle
(c)
Figure 4.3: (a) The orbit of Jupiter’s initial eccentricity (grey, dashed) and final ec-
centricity (red). (b) The orbit of Jupiter’s initial semi-major axis (dashed, grey) and
final semi-major axis (red). (c) The orbit of Jupiter’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (red). In (a), (b) and (c) , Jupiter’s
initial position is (q1(-smax), q2(-smax)) = (1.0489, 0).
39
planet. The orbital period is given by
TJ = 2pi
√
a3J
m1
.
For the initial value of aJ , the orbital period is 7.82391. The final value of the semi-major
axis here is 0.9690063, yielding a shorter orbital period of 7.463.
The argument of the periapsis also changes as the star passes by. The argument of
the periapsis is given by
tan(ωJ) = −eJ1
eJ2
where e = (eJ1 , eJ2). That is, eJ1 and eJ2 make up the components of the eccentricity,
or Laplace, vector. For the initial value of ωJ , the argument of the periapsis is 0 degrees
as the initial orbit begins with the major axis along the horizontal axis. The argument
of the periapsis increases to 19.1066◦ (see figure 4.3c).
In the simulations featured in figure 4.2 and figure 4.3, the initial conditions of
Jupiter’s orbit were set to
(q1(-smax), q2(-smax)) = (1.0489, 0).
This initial position informs where Jupiter is during the passing star’s closest approach
along the q1 axis. This can have an effect on how exactly Jupiter’s orbit changes as
the star passes. Starting at various points along Jupiter’s initial orbit and running the
simulation leads to different solutions. These solutions hold varying amounts of changes
to the eccentricity, semi-major axis, period, and argument of the periapsis.
Table 4.3 details the the result of running the simulation for 16 different starting
positions of Jupiter (see columns 1 and 2), equally distanced over the initial (current)
orbit of Jupiter [58]. In the table, the various final values of eccentricity, along with the
40
q1(-smax) q2(-smax) eJ(smax) % ∆aJ TJ(smax) ωJ(smax)
1.0489 0.0 0.0579576 -3.09937 7.463 19.1066
0.97278 0.382226 0.0542911 -3.40961 7.42719 23.2184
0.756007 0.706261 0.0486914 -3.16083 7.4559 23.5845
0.431583 0.922774 0.0446509 -2.41727 7.54194 18.9889
0.0489 0.998804 0.0448643 -1.30733 7.67098 12.4204
-0.333783 0.922774 0.0487823 -0.001156 7.82377 9.01296
-0.658207 0.706261 0.0534277 1.31278 7.97848 10.2666
-0.87498 0.382226 0.0562559 2.455 8.11378 14.5767
-0.9511 0.0 0.0561287 3.27651 8.21156 19.8478
-0.87498 -0.382226 0.0530516 3.66649 8.25812 24.2869
-0.658207 -0.706261 0.0481513 3.55609 8.24493 25.822
-0.333783 -0.922774 0.0437847 2.93012 8.17029 22.6209
0.0489 -0.998804 0.0430586 1.84818 8.0418 15.485
0.431583 -0.922774 0.0470903 0.463018 7.87831 9.94663
0.756007 -0.706261 0.0531985 -0.991304 7.70786 9.73182
0.97278 -0.382226 0.0575717 -2.25191 7.56112 13.6925
Table 4.3: Final value of eccentricity (eJ(smax)), percentage change in the semi-major
axis (% ∆aJ), final value of orbital period (TJ(smax)), and the final value of the argu-
ment of the periapsis (ωJ(smax)) for the simulations with initial conditions found in the
first two columns, along with those found in table 4.2. Recall that the initial value of eJ
is 0.0489, the initial value for aJ is 1, the initial value of TJ is 7.82391, and the initial
value of ωJ is 0 degrees.
percent change of the semi-major axis, the new value of the period of the orbit, and the
final argument of the periapsis, are listed. When the percent change of the semi-major
axis is positive, the axis increased by this percent. Similarly, when the percent change
of the semi-major axis is negative, this value decreased by this percent. Additionally,
the argument of the periapsis is given in degrees.
Recall, from table 4.2, that the initial value of eJ is 0.0489 and the initial value
for aJ is 1. In table 4.3, the eccentricity can be seen to both increase and decrease,
with the most dramatic change being the jump to 0.0579576, with initial conditions
41
-1.0 -0.5 0.5 1.0 q1
-1.0
-0.5
0.5
1.0
q2
Figure 4.4: Jupiter’s initial orbit (grey, dashed) and final orbit (red) are drawn in
the (q1, q2) plane, showing the effects of the passing star. Jupiter’s initial position is
(q1(-smax), q2(-smax)) = (−0.87498,−0.382226).
-400 -200 200 400 s
0.02
0.04
0.06
0.08
0.10
eccentricity
(a)
-400 -200 200 400 s
0.6
0.8
1.0
1.2
1.4
semimajoraxis
(b)
-400 -200 200 400 s
5
10
15
20
25
perihelion angle
(c)
Figure 4.5: (a) The orbit of Jupiter’s initial eccentricity (grey, dashed) and final ec-
centricity (red). (b) The orbit of Jupiter’s initial semi-major axis (dashed, grey) and
final semi-major axis (red). (c) The orbit of Jupiter’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (red). In (a), (b) and (c) , Jupiter’s
initial position is (q1(-smax), q2(-smax)) = (−0.87498− 0.382226).
42
-1.0 -0.5 0.5 1.0 q1
-1.0
-0.5
0.5
1.0
q2
Figure 4.6: Jupiter’s initial orbit (grey, dashed) and final orbit (red) are drawn in
the (q1, q2) plane, showing the effects of the passing star. Jupiter’s initial position is
(q1(-smax), q2(-smax)) = (−0.658207,−0.706261).
-400 -200 200 400 s
0.02
0.04
0.06
0.08
0.10
eccentricity
(a)
-400 -200 200 400 s
0.6
0.8
1.0
1.2
1.4
semimajoraxis
(b)
-400 -200 200 400 s
5
10
15
20
25
perihelion angle
(c)
Figure 4.7: (a) The orbit of Jupiter’s initial eccentricity (grey, dashed) and final ec-
centricity (red). (b) The orbit of Jupiter’s initial semi-major axis (dashed, grey) and
final semi-major axis (red). (c) The orbit of Jupiter’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (red). In (a), (b) and (c) , Jupiter’s
initial position is (q1(-smax), q2(-smax)) = (−0.658207,−0.706261).
43
(q1(-smax), q2(-smax)) = (1.0489, 0.0) (see figure 4.2). The biggest change in the semi-
major axis (and thus the biggest change in the orbital period) is found to be by an in-
crease of 3.66649% in the case where the initial conditions are (q1(-smax), q2(-smax)) =
(−0.87498,−0.382226) (see figures 4.4 and 4.5). Though this does not correspond to
the biggest change in the argument of the periapsis, it does correspond to the second
largest change in the argument of the periapsis. The largest change in the argument of
the periapsis occurs when (q1(-smax), q2(-smax)) = (−0.658207,−0.706261) when the
orbit rotates by almost 26◦ (see figures 4.6 and 4.7).
Not only do we observe several changes to the orbit of Jupiter here, it is important
to note that this would have an effect on Earth’s orbit as well as the secular fundamental
frequencies would be affected. The changes in the argument of the periapsis in particular
would have a strong effect on the gi frequencies, the precession of the periapsis, as this
would push these out of their regular oscillations.
4.3 Effects on the Orbit and Climate of Earth
As the Kuiper Belt is within the rage of 1800 AU, a passing star of the size considered
here will also have a direct effect on the orbit of Earth. Equation (2.24) is once again
solved but with initial conditions related to Earth given by
z1(τ) =
1√
2
√
aE(1 + eE) cos
(√
m1
aE
τ
2
)
(4.10)
z2(τ) =
1√
2
√
aE(1− eE) sin
(√
m1
aE
τ
2
)
(4.11)
w1(τ) = − 1√
2
√
m1
aE
√
aE(1 + eE) sin
(√
m1
aE
τ
2
)
(4.12)
w2(τ) =
1√
2
√
m1
aE
√
aE(1− eE) cos
(√
m1
aE
τ
2
)
, (4.13)
44
where
s =
1
2
√
m1
aE
τ.
The values of aE and eE can be found in table 4.2, along with the all other parameters
for this simulation. As the Earth is located closer to the Sun, the changes in the
orbit are far less noticeable (see figure 4.8). In fact, the orbit does not even appear to
change. This is not true of the semi-major axis, the eccentricity, and the argument of
the periapsis however. The changes are small, but do exist for those orbital elements.
The final value for the eccentricity of the orbit changes to 0.0176979 (see figure 4.9a),
an increase of approximately 0.001. Additionally, the final value for the semi-major axis
value is 0.192035, an increase of approximately 0.018% (see figure 4.9b). The change
in the semi-major axis results in a change to the orbital period of the planet. For the
initial value of aE , the orbital period is 0.658227. The final value of the semi-major axis
here yields a longer orbital period of 0.658406. Lastly, the final value of the argument of
the periapsis is 21.5522◦. Note that there is a bit of numerical error in the computation
of the argument of the perihelion here.
In the simulations featured in figures 4.8 and 4.9, the initial conditions of Earth’s
orbit were set to
(q1(-smax), q2(-smax)) = (0.195206, 0).
As in the previous section, the initial position of Earth informs how the orbit changes
as the passing star approaches. The different solutions that come for the different initial
conditions result in varying amounts of changes to the eccentricity and semi-major axis
values. Table 4.4 details the result of running the simulation for 16 different starting
positions of Earth (see columns 1 and 2), equally distanced over the initial (current)
orbit of Earth [60].
The largest change in the eccentricity and semi-major axis are observed for the initial
45
-0.2 -0.1 0.1 0.2 q1
-0.2
-0.1
0.1
0.2
q2
Figure 4.8: Earth’s initial orbit (grey, dashed) and final orbit (green) are drawn in
the (q1, q2) plane, showing the effects of the passing star. Earth’s initial position is
(q1(-smax), q2(-smax)) = (0.195206, 0).
-400 -200 200 400 s
0.015
0.020
0.025
eccentricity
(a)
-400 -200 200 400 s
0.05
0.10
0.15
0.20
0.25
0.30
semimajoraxis
(b)
-400 -200 200 400 s
5
10
15
20
perihelion angle
(c)
Figure 4.9: (a) The orbit of Earth’s initial eccentricity (grey, dashed) and final ec-
centricity (green). (b) The orbit of Earth’s initial semi-major axis (dashed, grey) and
final semi-major axis (green). (c) The orbit of Earth’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (green). In (a),(b) and (c), Earth’s
initial position is (q1(-smax), q2(-smax)) = (0.195206, 0).
46
q1(-smax) q2(-smax) eE(smax) % ∆aE TE(smax) ωE
0.195206 0.0 0.0176979 0.0181964 0.658406 21.5522
0.180591 0.073465 0.017684 0.01408 0.658366 21.5898
0.138971 0.135746 0.0176666 0.00907275 0.658316 21.5946
0.0766816 0.17736 0.0176499 0.00423917 0.658268 21.5698
0.0032064 0.191973 0.0176366 3.30223 ×10−4 0.65823 21.5211
-0.0702688 0.17736 0.0176293 -0.00227686 0.658204 21.4563
-0.132558 0.135746 0.0176283 -0.00310203 0.658196 21.3892
-0.174178 0.073465 0.0176324 -0.00198253 0.658207 21.3292
-0.188794 0.0 0.0176411 7.6627×10−4 0.658234 21.284
-0.174178 -0.073465 0.0176526 0.00489697 0.658275 21.2561
-0.132558 -0.135746 0.0176661 0.00966544 0.658322 21.2502
-0.0702688 -0.17736 0.0176801 0.0144073 0.658369 21.2666
0.0032064 -0.191973 0.0176929 0.0184611 0.658409 21.3039
0.0766816 -0.17736 0.0177027 0.0211543 0.658435 21.3592
0.138971 -0.135746 0.0177077 0.0220829 0.658445 21.4258
0.180591 -0.073465 0.0177063 0.0210528 0.658434 21.4942
Table 4.4: Final value of eccentricity (e(smax)), percentage change in the semi-major
axis (% ∆aE), final value of orbital period (TE(smax)), and the final value of the ar-
gument of the periapsis (ωE(smax)) for the simulations with initial conditions found in
the first two columns, along with those found in table 4.2. Recall that the initial value
of eE is 0.0167, the initial value for aE is 0.192, the initial value of TE is 0.658227, and
the initial value of ωE is 0 degrees.
conditions (q1(-smax), q2(-smax)) = (0.180591,−0.073465) where e becomes 0.0177063
and the semi-major axis changes by approximately 0.022% (see figures 4.10 and 4.11).
The largest change in the argument of the periapsis occurs when (−0.0702688,−0.17736)
when the orbit rotates by 21.2666◦ (see figures 4.12 and 4.13).
Laskar et. al showed that the variations in the semi-major axis of Earth are so small
that they will not create noticeable change in the average distance from the Earth to
the Sun over the last 250 million years [32]. The values seen in table 4.4 are about ten
times larger than those observed by Laskar et. al. in Earth’s natural oscillations (see
Figure 11 in [32]). This indicates that a star passing by the Kuiper Belt would cause a
noticeable change to the semi-major axis.
47
-0.2 -0.1 0.1 0.2 q1
-0.2
-0.1
0.1
0.2
q2
Figure 4.10: Earth’s initial orbit (grey, dashed) and final orbit (green) are drawn in
the (q1, q2) plane, showing the effects of the passing star. Earth’s initial position is
(q1(-smax), q2(-smax)) = (0.180591,−0.073465).
-400 -200 200 400 s
0.015
0.020
0.025
eccentricity
(a)
-400 -200 200 400 s
0.05
0.10
0.15
0.20
0.25
0.30
semimajoraxis
(b)
-400 -200 200 400 s
5
10
15
20
perihelion angle
(c)
Figure 4.11: (a) The orbit of Earth’s initial eccentricity (grey, dashed) and final ec-
centricity (green). (b) The orbit of Earth’s initial semi-major axis (dashed, grey) and
final semi-major axis (green). (c) The orbit of Earth’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (green). In (a),(b) and (c), Earth’s
initial position is (q1(-smax), q2(-smax)) = (0.138971,−0.135746).
48
-0.2 -0.1 0.1 0.2 q1
-0.2
-0.1
0.1
0.2
q2
Figure 4.12: Earth’s initial orbit (grey, dashed) and final orbit (green) are drawn in
the (q1, q2) plane, showing the effects of the passing star. Earth’s initial position is
(q1(-smax), q2(-smax)) = (−0.0702688,−0.17736).
-400 -200 200 400 s
0.015
0.020
0.025
eccentricity
(a)
-400 -200 200 400 s
0.05
0.10
0.15
0.20
0.25
0.30
semimajoraxis
(b)
-400 -200 200 400 s
5
10
15
20
perihelion angle
(c)
Figure 4.13: (a) The orbit of Earth’s initial eccentricity (grey, dashed) and final ec-
centricity (green). (b) The orbit of Earth’s initial semi-major axis (dashed, grey) and
final semi-major axis (green). (c) The orbit of Earth’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (green). In (a),(b) and (c), Earth’s
initial position is (q1(-smax), q2(-smax)) = (−0.0702688,−0.17736).
49
4.3.1 Global Mean Temperature Changes
Considering the changes that are made to the eccentricity and semi-major axis, it is
worthwhile to consider the kinds of changes that would be observed in the tempera-
ture distribution from these specific changes. Using the Budyko-Widiasih model (see
chapter 3), the global mean temperature will be plotted and analyzed.
The eccentricity of the Earth’s orbit naturally oscillates from 0.00 to 0.06 on a
100,000 year cycle [32]. Oscillations observed in table 4.4 are within this range. These
changes in the eccentricity have an effect on the global mean temperature, given by
T
∗
(η) =
Q(aE , eE) (1− α(η))−A)
B
, (4.14)
where
α(η) = αH2O
η∫
0
s6(y) dy + αice
1∫
η
s6(y) dy,
as it affects the incoming solar radiation. Let
∆GMT =
Q(0, 0.0167) (1− α(η))−A)
B
− Q(aE , eE) (1− α(η))−A)
B
. (4.15)
indicate the change in the global mean temperature from the current orbital parameter
values. As the final values for the eccentricity found in table 4.4 do not leave the range
of the regular oscillations of the Earth’s eccentricity, the range of the eccentricities that
Laskar et. al. gave are studied [32]. In figure 4.14, equation (4.15) is plotted with
eE = 0, 0.02, 0.04, and 0.06
and holding aE = 1. Changes in these values range from the order of 10
−3 to 10−1. The
changes are quite subtle around this range of differences. As such it can be assumed
50
that the changes for the values found in table 4.4 would be even more subtle, as those
changes are in the range of 0.001 for the eccentricity.
The effects on the global mean temperature from changes to the semi-major axis
are considered as well. Literature often assumes that the semi-major is constant, and
the range over which it oscillates to be negligibly small. Consider instead the changes
seen in table 4.4 as these changes are higher than those observed in [32]. In figure 4.15,
equation (4.15) is plotted holding eE = 0.0167 and
aE = 1.000220829, 1.0000966544, 0.9999801747, and 0.999969AU.
Change here range on a smaller interval, closer to the order of 10−2 or 10−3. These
images suggest that changes in the semi-major axis (aE) have a stronger effect on the
global mean temperature than changes to the eccentricity (eE). This is due to the
fact that the global mean temperature changes resulting from changing aE (figure 4.15)
and those changes resulting from changes to the eccentricity eE (figure 4.14) were on a
similar scale for vertical axis, but not the horizontal.
4.3.2 Equilibrium Temperature Profile Changes
For the most dramatic change in aE and eE seen in Earth as the star passes in ta-
ble 4.4, the equilibrium temperature profile changes can also be considered. Recall the
equilibrium temperature profile from Budyko’s model is given by equation (3.11). In
figure 4.16, both the equilibrium temperature profile for Earth’s current parameters
and the most dramatic shift seen in Earth’s parameters, aE = 1.000220829 AU and
eE = 0.0177077, are drawn. Only dramatic changes in the eccentricity and semi-major
axis are considered here as these two have an effect in the Budyko-Widiasih model.
Throughout the two temperature profiles, there is over a 1◦C difference figure 4.16.
51
0.2 0.4 0.6 0.8 1.0
η
0.005
0.010
0.015
ΔGMT
(a)
0.2 0.4 0.6 0.8 1.0
η
-0.006
-0.004
-0.002
ΔGMT
(b)
0.2 0.4 0.6 0.8 1.0
η
-0.08
-0.06
-0.04
-0.02
ΔGMT
(c)
0.2 0.4 0.6 0.8 1.0
η
-0.20
-0.15
-0.10
-0.05
ΔGMT
(d)
Figure 4.14: Change in global mean temperature as a function of the location of the ice
line, η. Here ∆GMT is the difference from the value of T
∗
with Earth’s current values
minus the value of T
∗
for various values of eE . Semi-major axis held constant at 1 AU.
(a) Eccentricity is eE = 0. (b) eE = 0.02. (c) eE = 0.04. (d) eE = 0.06.
52
0.2 0.4 0.6 0.8 1.0
η
-0.05
-0.04
-0.03
-0.02
-0.01
ΔGMT
(a)
0.2 0.4 0.6 0.8 1.0
η
-0.025
-0.020
-0.015
-0.010
-0.005
ΔGMT
(b)
0.2 0.4 0.6 0.8 1.0
η
0.001
0.002
0.003
0.004
0.005
ΔGMT
(c)
0.2 0.4 0.6 0.8 1.0
η
0.002
0.004
0.006
0.008
ΔGMT
(d)
Figure 4.15: Change in global mean temperature as a function of the location of the
ice line, η. Here ∆GMT is the difference from the value of T
∗
with Earth’s current
values minus the value of T
∗
for various values of aE . Eccentricity held constant at
0.0167. (a) Semi-major axis is aE = 1.000220829 AU. (b) aE = 1.0000966544 AU. (c)
aE = 0.9999801747 AU. (d) aE = 0.999969 AU.
53
0.2 0.4 0.6 0.8 1.0
y
-20
-10
10
20
30
40
Equilibrium Temperature
Figure 4.16: Earth’s current equilibrium temperature profile (grey, dashed) and equi-
librium temperature profile when a = 1.000220829 AU and e = 0.0177077 (green).
Temperature given in degrees Celsius. Ice line is located at 66.5◦ N.
This sort of change would be present in the sediment core data [44]. A cooling of 1
◦C is equivalent to a estimated δ18O increase of 0.22% [20]. This is an effect that would
have wide ranging repercussions on Earth. Even a 1 ◦C global change is significant
because it takes an enormous amount of heat to warm all the oceans, atmosphere, and
land by that much. In the past, a 2 to 3 ◦C drop was plunged the Earth into the Little
Ice Age [62]. A 3 to 7 ◦C drop was enough to bury a large part of North America under
ice during the Last Glacial Maximum [19]. The difference of 1.79156 ◦C observed at the
equator could be enough to dramatically change Earth’s climate.
4.3.3 Bifurcation Diagram in the Semi-Major Axis
The bifurcation diagram of the system with respect to the semi-major axis a is given in
figure 4.17. The black curve is obtained by solving for the parameter a when T ∗(η, t) is
set to the critical temperature Tc. The large blue arrows indicate the dynamics of the
ice line. If a pair of a and an initial ice line η is chosen to the right of the curve, then
54
0.96 0.98 1.00 1.02 1.04 1.06
a
0.2
0.4
0.6
0.8
1.0
η
Figure 4.17: Bifurcation diagram of the semi-major axis aE . Here eE = 0.0167, the
current value of Earth’s eccentricity. All other parameters are in table 3.1. The current
value of the Earth’s semi-major axis is 1 and marked in green. The filled dot represents
the stable equilibrium, and the open dot the unstable equilibrium.
the value of T ∗(η, t) will be greater than the critical temperature Tc. In this case, the
ice line equation dictates that the ice line moves up toward the pole; hence the upward
blue arrow to the right of the black curve. In a similar fashion, choosing a combination
left of the black curve will move the ice line toward the equator.
The value a = 1 AU is highlighted in figure 4.17. With this value, there is an
unstable big ice cap equilibrium at η ≈ 0.262 and a stable small ice cap equilibrium at
η ≈ 0.934. Furthermore, the stability of the branches above the saddle node bifurcation
point for the values near a = 0.9713 and η = 0.629 is indicated by a solid black curve
and the instability of the lower branch is represented by the dashed black curve.
There are also two bifurcation points when a ≈ 1.0517 (corresponding to η = 0)
and a ≈ 1.0149 (corresponding to η = 1), where the determination of the bifurcation
55
type is less straightforward. The physically relevant condition that the ice line is not
admissible outside of the unit interval gives a restriction on the dynamics of the ice line
at the boundaries. That is [56],
at the equator: η = 0,
dη
dt
= max {0, ε (T (η, t)− Tc)}
at the pole: η = 1,
dη
dt
= min {0, ε (T (η, t)− Tc)} .
This restriction determines the stability of the branches in the bifurcation diagram at
the physical boundary. The purple solid lines in figure 4.17 at the equator suggests a
ice covered planet (often referred to as “Snowball Earth”) whereas the solid line at the
pole suggests an ice-free Earth. The purple dashed lines in the same figure indicate an
unstable branch of the bifurcation diagram for a.
The stable equilibrium for Earth’s current semi-major axis value, aE = 1 AU, is when
the ice line is located at 68.9938◦ N (filled green dot in figure 4.17). If the semi-major
axis changes to aE = 1.000220829 AU, the ice line moves to 69.1728
◦ N. Approximating
Earth as a sphere, the surface area of ice lost here would be 745790 km2.
4.4 Discussion
In this chapter, the effects of a star passing at the edge of the Kuiper Belt on the orbits
of Jupiter and Earth were studied. The initial position of the planet had an effect on
the changes to the orbit that would be seen. If a star about half the size of the Sun were
to pass by the edge of the Kuiper Belt, there would be noticeable changes to the orbital
elements of Jupiter’s orbit, as can be observed from table 4.3, and subtler changes on
the orbital elements of Earth’s orbit, as can be observed from table 4.4. The effects that
these changes would have on the temperature distribution on Earth were also evaluated.
The equilibrium temperature profile appears to move over 1◦ C, close to 2◦ C at many
56
points, which is a noticeable difference. Lastly, a bifurcation diagram was drawn to
observe the way the semi-major axis length would affect the location of Earth’s ice line.
In the case of Jupiter’s orbit, there were some significant changes that occurred. The
eccentricity, semi-major axis, and argument of the periapsis saw changes that would be
visible in the sedimentary core data. While the range of the changes to eccentricity that
were observed were within the natural oscillations of the eccentricity parameter, the
semi-major axis changed at a rate higher than its natural oscillations [32]. This would
surely disrupt the 400,000 year eccentricity cycle for planetary orbits. Additionally the
changes to the argument of the periapsis would have an effect on the precession of the
periapsis frequency, altering Jupiter’s relationship to Earth.
In the case of Earth’s orbit, without the compounding effects from the changes
observed to Jupiter’s orbit, the changes are quite small. There is still an effect on the
equilibrium temperature profile, but the line does not significantly change, only moving
0.2◦ N even in the most dramatic of circumstances. This reveals that perhaps such an
occurrence, without being able to account for the effects due to the changes of Jupiter,
would not significantly contribute to glaciation or deglaciation. The biggest effects
would be from the changes to the secular fundamental frequencies of Jupiter’s orbit.
Give the changes observed to Jupiter’s orbit, it can be safely assumed that a star
of mass approximately half the size of the Sun did not pass within the Kuiper Belt
during the last 65 million years as this would be observable in the climate record. The
climate record prior to 65 million years ago has large gaps, but the work by Olsen et.
al. extending back to the Mesozoic era (250 to 65 million years ago) shows that the
basic Milankovitch cycles were much of what we see today, indicating that a Kuiper
Belt flyby probably has not occurred during the last 250 million years.
Chapter 5
Passage near the Oort Cloud
The distance that the passing star is in relation to the Sun and the planet will determine
the kinds of effects that are observed. While marked differences in the orbits of Earth
and Jupiter were observed in chapter 4 when the star passed by the edge of the Kuiper
Belt, it begs the question: how far away the star can pass and still cause observable
changes to orbits of planets in the Solar System?
The Oort cloud is a theoretical cloud of predominantly icy objects proposed to
surround the Sun at distances ranging from approximately 3,000 to 50,000 AU 1 [57].
The Oort cloud lies beyond the heliosphere, in interstellar space. The Kuiper belt on
the other hand is less than one thousandth as far from the Sun as the Oort cloud. The
outer limit of the Oort cloud defines the boundary of the Solar System and the extent
of the Sun’s region of attraction, or Hill sphere. Unlike the orbits of the planets and
the Kuiper Belt, which lie mostly in the same flat disk around the Sun, the Oort Cloud
is believed to be a giant spherical shell surrounding the rest of the Solar System.
In chapter 4, it was shown that the Earth would be affected by a star passing within
1The bounds of the Oort cloud are not distinct and are often disputed. Some say it starts at 10,000
AU and some say it ends out at 100,000 AU. In this paper, 3,000 to 50,000 AU will be considered the
Oort cloud.
57
58
approximately 1800 AU from the Sun. In this chapter, distances ranging from 100 AU
all the way out to 50,000 AU will be considered for the closest approach of the passing
star. The study will focus on the passing star’s effects on Jupiter’s orbit in this chapter
as Earth cannot be affected alone at some of these more distant ranges. In section 1,
the effect of the distance of the star on Jupiter’s orbital elements is considered. Then
in section 2, potential possible effects on the climate due to such changes are analyzed.
5.1 Varying Star’s Distance with Constant Initial Position
of the Planet
The distance of the star is controlled by the parameter a, the semi-major axis of the
hyperbolic orbit of the star. By changing the position of the closest approach, the
effect on the orbital elements as a result of the distance of the can be studied. Other
parameters must be held constant, so all orbits will start along the horizontal axis at
-smax. Initial conditions remain similar as a result and can be found in table 5.1.
Parameter Name Value
e eccentricity of orbit of passing star 2
E hyperbolic anomaly of orbit of passing star -3
aJ semi-major axis of orbit of Jupiter 1
eJ eccentricity of orbit of Jupiter 0.0489
TJ period of Jupiter’s orbit 7.82391
ωJ argument of the periapsis of Jupiter’s orbit 0
m1 mass of Sun 0.6449299
m2 mass of passing star 0.3550701
(z1, z2) initial position of Jupiter (1.0489, 0.0)
Table 5.1: Initial conditions of the HR3BP for Jupiter as the third body in an elliptic
orbit starting on the horizontal axis. The system has been scaled so that Jupiter could
be a distance of 1 unit away.
59
As the star is moved further and further away, it makes sense that the effects on the
orbit of Jupiter are diminishing. In figures 5.1 to 5.4, the orbit, eccentricity,semi-major
axis, and argument of the periapsis of Jupiter are drawn from the values a = 100 AU
(figures 5.1 and 5.2a to 5.2c) and for a = 200 AU (figures 5.3 and 5.4a to 5.4c). In
figure 5.1, the initial and final orbit of Jupiter do not appear to change at all though the
eccentricity and semi-major axis plots in figures 5.2a and 5.2b do show small changes
that are still perceivable to the eye. Even in figure 5.4a, there is a small bump that can
be seen near the beginning of the graph as the star passes. Notice additionally that the
argument of the periapsis graph is accumulating more error as the value of a increases,
making it a less reliable metric.
These effects are summarized in table 5.2. Around 1000 AU, the star’s effect partic-
ularly appears to drop off. Though the changes have all been small, there is a noticeable
change in the order of the percent change in the semi-major axis of the orbit of Jupiter
between the star’s closest approach being at 750 AU versus at 1000 AU.
After a = 30000 AU, the value of the eccentricity of Jupiter’s orbit appears constant,
as do the period of Jupiter’s orbit. While the value aJ appears to continue changing,
the value is so small that it has no effect on TJ . After a = 30000 AU, the value of
ωJ stays constant for four decimal places. Table 5.2 implies that after 30000 AU, the
effects of the passing star are imperceptible. Though it appears that changes are still
happening, these minuscule changes can be the result of numerical error. This means
that subtle changes on the orbit of Jupiter due to a passing star may still be present
so long as the star passes within approximately 30000 AU. While there are still subtle
changes happening around a = 20000 AU and a = 10000 AU, those changes are also
quite small and could be very difficult to notice in the geological record.
60
-1.0 -0.5 0.5 1.0 q1
-1.0
-0.5
0.5
1.0
q2
Figure 5.1: Jupiter’s initial orbit (grey, dashed) and final orbit (red) are drawn in the
(q1, q2) plane, showing the effects of the passing star. The star’s closest approach is at
100 AUs.
-400 -200 200 400 s
0.02
0.04
0.06
0.08
0.10
eccentricity
(a)
-400 -200 200 400 s
0.6
0.8
1.0
1.2
1.4
semimajoraxis
(b)
-400 -200 200 400 s
1
2
3
4
5
perihelion angle
(c)
Figure 5.2: (a) The orbit of Jupiter’s initial eccentricity (grey, dashed) and final ec-
centricity (red). (b) The orbit of Jupiter’s initial semi-major axis (dashed, grey) and
final semi-major axis (red). (c) The orbit of Jupiter’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (red). In (a),(b) and (c), the star’s
closest approach occurs at 100 AUs.
61
-1.0 -0.5 0.5 1.0 q1
-1.0
-0.5
0.5
1.0
q2
Figure 5.3: Jupiter’s initial orbit (grey, dashed) and final orbit (red) are drawn in the
(q1, q2) plane, showing the effects of the passing star. The star’s closest approach is at
200 AUs.
-400 -200 200 400 s
0.02
0.04
0.06
0.08
0.10
eccentricity
(a)
-400 -200 200 400 s
0.6
0.8
1.0
1.2
1.4
semimajoraxis
(b)
-400 -200 200 400 s
0.5
1.0
1.5
perihelion angle
(c)
Figure 5.4: (a) The orbit of Jupiter’s initial eccentricity (grey, dashed) and final ec-
centricity (red). (b) The orbit of Jupiter’s initial semi-major axis (dashed, grey) and
final semi-major axis (red). (c) The orbit of Jupiter’s initial argument of the periapsis
(dashed, grey) and final argument of the periapsis (red). In (a),(b) and (c), the star’s
closest approach occurs at 200 AUs.
62
a eJ(smax) % ∆aJ TJ(smax) ωJ(smax)
100 0.0505012 -0.775703 7.73305 5.28593
200 0.0492389 -0.191121 7.80149 1.34754
300 0.0490427 -0.0833261 7.81413 0.606828
400 0.0489771 -0.0457074 7.81854 0.347067
500 0.0489473 -0.0283277 7.82058 0.226833
750 0.0489182 -0.0111906 7.82259 0.108161
1000 0.0489081 -0.00520196 7.8233 0.0666694
2000 0.0488984 0.000564855 7.8239 0.0267012
3000 0.0488966 0.0016314 7.8241 0.0193084
4000 0.048896 0.00200461 7.82414 0.0167222
5000 0.0488957 0.00217728 7.82416 0.0155255
7500 0.0488954 0.00234748 7.82418 0.0143421
10000 0.0488953 0.00240714 7.82419 0.0139286
20000 0.0488952 0.00246467 7.82419 0.0135299
30000 0.0488952 0.00247562 7.8242 0.0134578
40000 0.0488952 0.00247935 7.8242 0.0134319
50000 0.0488952 0.0024807 7.8242 0.0134183
100000 0.0488952 0.00248307 7.8242 0.0134024
Table 5.2: Final value of eccentricity (eJ(smax)), percentage change in the semi-major
axis (% ∆aJ), final value of orbital period (TJ(smax)), and final value of the argument
of the periapsis (ωJ(smax)) for the simulations with distance of closest approach given
by the first column, a, along with those found in table 5.1. Recall that the initial value
of eJ is 0.0489, the initial value for aJ is 1, the initial value of TJ is 7.82391, and the
initial value of ωJ is 0 degrees.
5.1.1 Global Mean Temperature Changes
To study more closely if effects of a distant passing star could be noticeable in the climate
record, the global mean temperature is be plotted. As the parameters for the Budyko
model are known for Earth, the assumption will be made that the changes observed to
Jupiter will occur at the same percentage to Earth. For example, consider the case of
a = 100 AU for the closest approach of the star. For eJ = 0.0505012, meaning that
the eccentricity of Jupiter’s orbit has increased by approximately 3.27%. Additionally,
it is known that a decreased by 0.775703%. Thus, eE will be increased by the same
63
percentage and eE ≈ 0.017247 and aE will be decrease by the same amount yielding
aE = 0.99224297 AU. These new values are plugged into equation (4.14).
0.2 0.4 0.6 0.8 1.0
η
0.5
1.0
1.5
2.0
ΔGMT
(a)
0.2 0.4 0.6 0.8 1.0
η
0.1
0.2
0.3
0.4
ΔGMT
(b)
0.2 0.4 0.6 0.8 1.0
η
0.05
0.10
0.15
0.20
ΔGMT
(c)
0.2 0.4 0.6 0.8 1.0
η
0.02
0.04
0.06
0.08
0.10
ΔGMT
(d)
Figure 5.5: Change in global mean temperature as a function of the location of the ice
line, η based off the values found in table 5.2. Here ∆GMT is the difference from the
value of T
∗
with Earth’s current values minus the value of T
∗
for various values of aE .
(a) a = 100 AU; (b) a = 200 AU; (c) a = 300 AU; (d) a = 400 AU.
In figure 5.5, the change in global mean temperature is plotted (see equation (4.15))
for the closest approach of the star being a = 100 AU (figure 5.5a), a = 200 AU
(figure 5.5b), a = 300 AU (figure 5.5c), and a = 400 AU (figure 5.5d).
For the case of a = 100 in figure 5.5a, there is a marked difference between the global
mean temperature for Earth’s current orbital elements versus the orbital elements after
the star passes. If the ice line were at the equator (η = 0), the global mean temperature
would decrease by 0.896963◦C and if the ice line were at the pole (η = 1), the global
64
mean temperature would decrease by 1.60509◦C. For the case of a = 200 in figure 5.5b,
the biggest difference between the two graphs is 0.468609◦C at the pole. In figure 5.5d
for a = 400 AU has a difference of 0.112154◦C for the global mean temperature if the
ice line were at the pole.
a = 100 a = 500 a = 1000 a = 5000 a = 10000
∆ GMT(0.0) 1.05965 0.0388468 0.00713469 -0.00298598 -0.00330124
∆ GMT(0.2) 1.26244 0.0462811 0.0085001 -0.00355742 -0.00393302
∆ GMT(0.4) 1.45808 0.0534533 0.00981736 -0.00410872 -0.00454252
∆ GMT(0.6) 1.63743 0.0600283 0.0110249 -0.0046141 -0.00510126
∆ GMT(0.8) 1.7883 0.0655591 0.0120407 -0.00503923 -0.00557128
∆ GMT(1.0) 1.89622 0.0695153 0.0127673 -0.00534333 -0.00590748
Table 5.3: The difference between the global mean temperature (GMT) before the
star passes to after the star passes at various ice line locations and for various closest
approaches of the star. If the number is negative, the global mean temperature increases
after the star passes, and if it is positive, the global mean temperature decreases after
the star passes.
Table 5.3 lists values of equation (4.15) at various ice line locations and for various
closest approaches of the star where the columns designate the distance of the closest
approach of the passing star. Even when the star is 10000 AU away, there are still
minute changes occurring in the global mean temperature profile. For a = 5000 AU
and a = 10000 AU, the global mean temperature appears to increase slightly whereas
for the smaller values of a the global mean temperature decreased. This is due to the
fact that starting with a = 2000 AU, the value of the semi-major axis increases slightly
when the star passes (as can be seen in table 5.2). As was previously noted, in particular
the changes are quite minor when the star approaches further away than 1000 AU as
all changes occur on an order of 10−3.
65
5.1.2 Equilibrium Temperature Profile Changes
As there are notable changes to the global mean temperature when the star passes at
100 AUs as its closest approach, the equilibrium temperature profile will be drawn for
this case. Recall that the equilibrium temperature profile is given by equation (3.11).
In this case, the semi-major axis shrinks and the eccentricity increases.
0.2 0.4 0.6 0.8 1.0
sin(lat)
-30-20
-10
10
20
30
Equilibrium Temperature
Figure 5.6: Earth’s current equilibrium temperature profile (grey, dashed) and equilib-
rium temperature profile when aE = 0.99224297 AU and eE ≈ 0.0172468 (green) which
occurs when the closest approach is at a = 100 AUs. Temperature given in degrees
Celsius. Ice line is located at 66.5◦ N here.
The equilibrium temperature profile given in figure 5.6 differs by over 2 ◦C at the
equator (η = 0) and by around 1.3 ◦C at the pole (η = 1). This is a dramatic enough
change to the the temperature profile that it would be spotted in the climate data.
5.2 Varying the Initial Position of the Planet
As seen from chapter 4, the initial position of Jupiter can have a significant effect on the
way that the way the orbital elements change as the rogue star passes. In this section,
the initial position of the planet will be varied along with the closest approach of the
passing star. The closest approach of the star was varied from a = 100 AU to a = 1000
66
AU and the orbit of Jupiter was split into 100 equally spaced starting points. The
starting position of Jupiter which resulted in the largest variation in either eccentricity,
semi-major axis, and the argument of the periapsis is listed in tables 5.4 to 5.6.
Table 5.4 holds the initial conditions at which the eccentricity showed the greatest
variation for the case of a found in the first column. Similarly, table 5.5 does this for the
greatest variation of the semi-major axis (in absolute deviation) and table 5.6 does this
for the greatest variation of the argument of the perihelion. While the largest variation
in the eccentricity is approximately 0.0001 larger than the eccentricity variation found
in table 5.2, the largest variation in the argument of the periapsis is over 1◦ larger and
the largest variation in the semi-major axis changes by over 0.1%. The change in the
argument of the periapsis in particular is the most interesting as this would change
Jupiter’s resonance with Earth, having an effect on the orbit and climate.
5.2.1 Global Mean Temperature Changes
The global mean temperature is considered for the changes in tables 5.4 and 5.5 as the
eccentricity and the semi-major axis length contribute to the incoming solar radiation.
Once again, as the parameters for the Budyko-Widiasih model are known for Earth, the
changes observed are assumed to happen in the same percentage. In figure 5.7, changes
in the global mean temperature are drawn for the first three rows of table 5.4 where
a = 100, 200, and 300 AU. In figure 5.8, changes in the global mean temperature are
drawn for the first three rows of table 5.5 where a = 100, 200, and 300 AU.
In figures 5.7 and 5.8, the global mean temperatures move with the value of the semi-
major axis. When the semi-major axis increases, the global mean temperature profile
increases whereas when the semi-major axis decreases, the global mean temperature
profile decreases. The biggest change occurs in figure 5.8a, where the star passes at
a = 100 AU and the largest change is in the semi-major axis and the argument of the
67
a (q1(−smax), q2(−smax)) eJ(smax) % ∆aJ TJ(smax) ωJ(smax)
100 (1.01748, -0.248392) 0.0507187 -0.6462 7.7481910 4.30786
200 (0.999957, -0.308647) 0.0493079 -0.148068 7.80653488 1.03733
300 (0.999957, -0.308647) 0.0490767 -0.0640961 7.81638448 0.46893
400 (0.999957, -0.308647) 0.048998 -0.0348672 7.81981389 0.269818
500 (0.999957, -0.308647) 0.048962 -0.0213799 7.82139652 0.177741
750 (-0.902157, 0.308647) 0.0489279 0.0145121 7.825608682 0.068049
1000 (-0.933387, 0.187157) 0.0489166 0.00994914 7.825073147 0.0367391
Table 5.4: Orbital elements after the star passes from varying distances (closest approach
of star given in column 1). The starting locations (column 2) that resulted in the largest
changes in the eccentricity are given.
a (q1(−smax), q2(−smax)) eJ(smax) % ∆aJ TJ(smax) ωJ(smax)
100 (-0.855927, -0.42527) 0.0489522 0.895905 7.9292828 6.38302
200 (-0.855927, -0.42527) 0.0488604 0.224227 7.8502352 1.57297
300 (-0.855927, -0.42527) 0.048884 0.100982 7.83575960 0.691083
400 (-0.855927, -0.42527) 0.0488946 0.0579254 7.83070453 0.383763
500 (-0.855927, -0.42527) 0.0488999 0.0380126 7.82836703 0.241752
750 (-0.855927, -0.42527) 0.0489053 0.0183561 7.826059846 0.101644
1000 (-0.855927, -0.42527) 0.0489073 0.0114788 7.825252675 0.0526429
Table 5.5: Orbital elements after the star passes from varying distances (closest approach
of star given in column 1). The starting locations (column 2) that resulted in the largest
changes in the semi-major axis are given.
a (q1(−smax), q2(−smax)) eJ(smax) % ∆aJ TJ(smax) ωJ(smax)
100 (-0.827407, -0.481177) 0.0487435 0.894821 7.9291550 6.402
200 (0.925207, 0.481177) 0.0488555 -0.218216 7.7983100 1.58209
300 (0.925207, 0.481177) 0.0488688 -0.0957713 7.81266861 0.708449
400 (0.953727, 0.42527) 0.0488917 -0.0527492 7.81771575 0.402492
500 (0.953727, 0.42527) 0.0488911 -0.0328775 7.82004736 0.261065
750 (0.953727, 0.42527) 0.0488909 -0.0132791 7.822347136 0.121253
1000 (0.978676, 0.367684) 0.0488937 -0.00635297 7.823159937 0.072523
Table 5.6: Orbital elements after the star passes from varying distances (closest approach
of star given in column 1). The starting locations (column 2) that resulted in the largest
changes in the argument of the periapsis are given.
68
0.2 0.4 0.6 0.8 1.0
η
0.5
1.0
1.5
ΔGMT
(a)
0.2 0.4 0.6 0.8 1.0
η
0.1
0.2
0.3
ΔGMT
(b)
0.2 0.4 0.6 0.8 1.0
η
0.05
0.10
0.15
ΔGMT
(c)
Figure 5.7: Change in global mean temperature as a function of the location of the ice
line, η based off the values found in table 5.4. Here ∆GMT is the difference from the
value of T
∗
with Earth’s current values minus the value of T
∗
. (a) (1.01748,−0.248392),
(b) (0.999957,−0.308647), (c) (0.999957,−0.308647).
69
0.2 0.4 0.6 0.8 1.0
η
-2.0
-1.5
-1.0
-0.5
ΔGMT
(a)
0.2 0.4 0.6 0.8 1.0
η
-0.5
-0.4
-0.3
-0.2
-0.1
ΔGMT
(b)
0.2 0.4 0.6 0.8 1.0
η
-0.25
-0.20
-0.15
-0.10
-0.05
ΔGMT
(c)
Figure 5.8: Change in global mean temperature as a function of the location of
the ice line, η based off the values found in table 5.4. Here ∆GMT is the differ-
ence from the value of T
∗
with Earth’s current values minus the value of T
∗
. (a)
(−0.855927,−0.42527), (b) (−0.855927,−0.42527), (c) (−0.855927,−0.42527)
70
periapsis. The semi-major axis changes continue to have a strong effect on the global
mean temperature. While the argument of the periapsis does not come into play in
equation (4.14), it affects the resonances in the Solar System.
5.2.2 Equilibrium Temperature Profile Changes
For the global mean temperature profiles studied in the previous section where a = 100
AU for the closest approach of the passing star, the equilibrium temperature profiles
are plotted. Figure 5.9a is for the case of the biggest change in the eccentricity and
figure 5.9b is for the case of the biggest change in the semi-major axis.
In all the curves of figure 5.9, the dashed, grey curve denotes the equilibrium temper-
ature profile for Earth’s current orbital elements with aE = 1 AU and eE = 0.0167. In
figure 5.9a, a = 100 AU for the initial condition (1.01748,−0.248392). The two profiles
differ by over 1.5◦ C, with a maximum difference of approximately 2.4◦ C when η = 0.
Figure 5.9b has a difference of around 1.1◦ C - 1.7◦ C in the equilibrium temperature
profile with Earth’s current orbit.
5.3 Discussion
In this chapter, the sensitivity of the Sun-planet-star system to the distance of the star’s
closest approach was analyzed. First, the effects on the orbital elements were noted,
followed by the possible effects those could have on Earth’s climate. As the passing
star moves further away from the Solar System, the perturbing effects on Jupiter’s orbit
are weakened, as expected. When the star’s closest approach is around 1000 AU, the
effect particularly appears to drop off as the order of the percent change in the semi-
major axis changes at this point. The semi-major axis is of particular note as it has
the strongest influence on the global mean temperature and equilibrium temperature
71
0.2 0.4 0.6 0.8 1.0
sin(lat)
-30
-20
-10
10
20
30
Equilibrium Temperature
(a)
0.2 0.4 0.6 0.8 1.0
sin(lat)
-30
-20
-10
10
20
30
Equilibrium Temperature
(b)
Figure 5.9: Earth’s current equilibrium temperature profile (grey, dashed) and equilib-
rium temperature profile when (a) (1.01748,−0.248392), a = 100 AU, biggest change
in eccentricity; (b) (−0.855927,−0.42527), a = 100 AU, biggest change in semi-major
axis. Temperature given in degrees Celsius. Ice line is located at 66.5◦ N here.
72
profiles of the planet. The changes in the argument of the periapsis are very minimal
at these distances, shifting at most approximately 5◦ - 6◦.
Different starting points were also considered. A smaller range of the closest ap-
proach was considered here as opposed to the earlier section as it was previously noted
that 1000 AU is where there appears to be a drop off in the effect of the passing star.
Larger effects by the passing star were observed at various starting locations. The
changes in the semi-major axis continued to have the largest effect on the planet’s tem-
perature profile. Some of the changes observed in this chapter could have a 1◦ C to 2◦
C change in the temperature profile, a significant change as noted in chapter 4.
As the parameters for the Budyko-Widiasih model that are known are for Earth,
the study of the climate in this chapter was based off assumptions of how the Earth’s
orbits might change as a result of Jupiter’s orbit. It did not take into direct account the
effects of the secular fundamental frequencies but instead considered that case of the
same percentage changes happening to Earth’s orbit that happened to Jupiter’s. The
global mean temperature profile had changes ranging from 2 ◦C for the case of a = 100
AU to around 0.1 ◦C when the star passes at a = 300 AU from the Sun. Equilibrium
temperature profiles shifted by around 1 ◦C - 2 ◦C in the most dramatic cases when
a = 100 AU.
Chapter 6
Discussion
In this work, the effects of the passage of a star on the orbits of a planet were studied
using a special case of the n-body problem called the hyperbolic restricted 3-body
problem (HR3BP). An application of these changes were then studied, looking at how
it could affect the climate on Earth.
One consideration to be made is the likelihood of a rogue star passing by the Solar
System. Based on results from the Gaia telescope’s second data release in April 2018,
an estimated 694 stars will possibly approach our Solar System to less than 1.543×1014
km (1.0315× 106 AU) in the next 15 million years. Of these, 26 have a greater than 50
% chance of coming within 3.12×1013 km (209000 AU) and another 7 within 1.51×1013
km (101000 AU) [2]. This number is likely much higher, due to the sheer number of
starts needed to be surveyed. The present rate of encounters within 206564 AU (1
parsec) is estimated to be 19.7 ± 2.2 per million years, a quantity expected to scale
quadratically with the encounter distance out to at least several parsecs [2].
The first part of this work focuses on the star passing within 50 AU of the Sun in
the Kuiper Belt. This distance is close enough to have an effect on Earth stronger than
that of Jupiter. This means that not only will the changes to Jupiter’s orbit exist, these
73
74
will compound with the changes to Earth’s orbit from the star. This effect was observed
to be dramatic enough that such a passage this close could be ruled out. An event such
as this would been easily spotted in the geological record and as no such effect has been
spotted, the conclusion stands that it probably did not occur. Indeed the number of
stars passing within the outer regions of the Kuiper Belt (50 AU) during the life of the
solar system is 0.0003, a highly unlikely occurrence.
In the second part of this work, the focus was on the sensitivity of the planetary
orbits to the star’s distance. The goal was to find where it is that minute changes are
still perceptible as the star moves out towards the Oort Cloud. The effects of the passing
star diminish as the star moves farther and farther away and the effects are unlikely to
be spotted. The number of stars passing within the outer reaches of the Oort Cloud
(50,000 AU) during the last 65 million years is about 75 and during the last 250 million
years is about 289. The number of stars passing within the inner reaches of the Oort
Cloud (3000 AU) during the last 250 million years is about 1. Over the life of the Solar
System, which is approximately 4.568 billion years old, this would be about 19 stars [5].
While there is a much higher probability that a star passed within this region of the
Solar System, the effects of the star would be negligible and difficult to perceive.
The results of this work, in addition to the cited survey of nearby stars, indicate
that it is unlikely that one of these stars passed close enough to the Solar System that
its effects would appear in the Earth’s geological record. There does remain the distinct
possibility that a rogue star did pass near enough to do so. Other interstellar objects
have passed along hyperbolic orbits within the Solar System, making it more likely that
a larger object could follow a similar path [53].
75
6.1 Future Work
From the results found here, more work is indicated on this question. There are several
directions that one can explore based on these results. One immediate extension that
arises is the expansion to the 3-dimensional case. The full 3-dimensional problem should
be considered to see what happens when the star approaches from out of the plane
could have effects on the obliquity, a Milankovitch cycle whose changes would have
strong climatic effects. Another immediate improvement in the study could be made
by studying the effects of some parameters that were held constant here. In this study,
both the eccentricity of the passing star’s orbit and its mass were held constant. The
eccentricity of the passing star’s orbit will change how close the rogue star would be to
the planetary system. Additionally, rogue stars are stars that are roaming across the
galaxy and therefore could be of any mass. One future study should focus on the effects
that the mass of the star have on the distance of perceivable changes to the orbits of
the planets in the Solar System.
References
[1] E. P. Aksenov, Averaged, restricted, circular three-body problem (averaged re-
stricted circular three body problem, analyzing small mass motion about two bodies
revolving about each other in circular orbits), Universitet Druzhby Narodov, Trudy,
21 (1967), pp. 184–202.
[2] C. A. L. Bailer-Jones, J. Rybizki, R. Andrae, and M. Fouesneau, New
stellar encounters discovered in the second gaia data release, Astronomy & Astro-
physics, 616 (2018), p. A37.
[3] A. Berger, M. Loutre, and J. L. Me´lice, Equatorial insolation: from pre-
cession harmonics to eccentricity frequencies, Climate of the Past, European Geo-
sciences Union (EGU), 2 (2006), pp. 131–136.
[4] V. V. Bobylev, Searching for stars closely encountering with the solar system,
Astronomy Letters, 36 (2010), pp. 220 – 226.
[5] A. Bouvier and M. Wadhwa, The age of the Solar System redefined by the oldest
Pb–Pb age of a meteoritic inclusion, Nature Geoscience, 3 (2010), pp. 637–641.
[6] M. I. Budyko, The effect of solar radiation variations on the climate of the earth,
Tellus, 21 (1969), pp. 611–619.
76
77
[7] P. Chy`lek and J. A. Coakley, Analytical analysis of a budyko-type climate
model, Journal of the Atmospheric Sciences, 32 (1975), pp. 675–679.
[8] J. M. Cors and J. Llibre, Chaos in the hyperbolic restricted 3-body problem, in
From Newton to Chaos, Springer, 1995, pp. 315–318.
[9] , The global flow of the hyperbolic restricted three-body problem, Archive for
rational mechanics and analysis, 131 (1995), pp. 335–358.
[10] , Qualitative study of the hyperbolic collision restricted three-body problem,
Nonlinearity, 9 (1996), p. 1299.
[11] T. J. Crowley and K. Kim, Milankovitch forcing of the last interglacial sea level,
Science, 265 (1994), pp. 1566–1568.
[12] , Comparison of longterm greenhouse projections with the geologic record, Geo-
physical Research Letters, 22 (1995), pp. 933–936.
[13] , Comparison of proxy records of climate change and solar forcing, Geophysical
Research Letters, 23 (1996), pp. 359–362.
[14] , Modeling the temperature response to forced climate change over the last six
centuries, Geophysical Research Letters, 26 (1999), pp. 1901–1904.
[15] T. J. Crowley, K. Kim, J. G. Mengel, and D. A. Short, Modeling
100,000-year climate fluctuations in pre-pleistocene time series, Science, 255 (1992),
pp. 705–707.
[16] T. J. Crowley, D. A. Short, J. G. Mengel, and G. R. North, Role of
seasonality in the evolution of climate during the last 100 million years, Science,
231 (1986), pp. 579–584.
78
[17] M. De Sanctis, M. Capria, and A. Coradini, Thermal evolution and differ-
entiation of edgeworth-kuiper belt objects, The Astronomical Journal, 121 (2001),
p. 2792.
[18] G. N. Duboshin, Celestial mechanics: Basic problems and methods, Izdatel’stvo
Nauka, 1975.
[19] A. S. Dyke, J. T. Andrews, P. U. Clark, J. H. England, G. H. Miller,
J. Shaw, and J. J. Veillette, The laurentide and innuitian ice sheets during
the last glacial maximum, Quaternary Science Reviews, 21 (2002), pp. 9–31.
[20] S. Epstein, R. Buchsbaum, H. A. Lowenstam, and H. C. Urey, Revised
carbonate-water isotopic temperature scale, Geological Society of America Bulletin,
64 (1953), pp. 1315–1326.
[21] M. B. Faintich, Application of the restricted hyperbolic three-body problem to a
star-sun-comet system, Celestial mechanics, 6 (1972), pp. 22–26.
[22] R. Fitzpatrick, An Introduction to Celestial Mechanics, Cambridge University
Press, 2012.
[23] J. D. Hays, J. Imbrie, and N. J. Shackleton, Variations in the earth’s orbit:
pacemaker of the ice ages, science, 194 (1976), pp. 1121–1132.
[24] H. Kaper and H. Engler, Mathematics and climate, vol. 131, Siam, 2013.
[25] J. Kepler, Astronomia Nova, Green Lion Press, 1609. Translated by W.H. Don-
ahue (2015).
[26] , Harmonice Mundi, vol. 6 of Gesammelte Werke, 1619. reprinted (1837).
[27] P. Kervella, F. The´venin, D. Se´gransan, G. Berthomieu, B. Lopez,
P. Morel, and J. Provost, The diameters of α centauri a and b-a comparison
79
of the asteroseismic and vinci/vlti views, Astronomy & Astrophysics, 404 (2003),
pp. 1087–1097.
[28] P. Kustaanheimo and E. L. Stiefel, Perturbation theory of kepler motion based
on spinor regularzation, Journal fu¨r die Reine und Angewandte Mathematik, 218
(1965), pp. 204–219.
[29] J. Laskar, The limits of earth orbital calculations for geological time-scale use,
Philosophical Transactions of the Royal Society of London. Series A: Mathematical,
Physical and Engineering Sciences, 357 (1999), pp. 1735–1759.
[30] , Chaos in the solar system, in Annales Henri Poincare´, vol. 4 (Supp. 2),
Springer, 2003, pp. 693–705.
[31] J. Laskar, A. Fienga, M. Gastineau, and H. Manche, La2010: a new orbital
solution for the long-term motion of the earth, Astronomy & Astrophysics, 532
(2011), p. A89.
[32] J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A. Correia, and
B. Levrard, A long-term numerical solution for the insolation quantities of the
earth, Astronomy & Astrophysics, 428 (2004), pp. 261–285.
[33] T. Levi-Civita, Sur la Re´gularisation du Proble`me des Trois Corps, Acta Math,
42 (1920), pp. 99–144.
[34] F.-P. Luis and P.-C. Ernesto, The symmetric elliptic and hyperbolic restricted
3-body problem on the unit circle, Journal of Geometry and Physics, 99 (2016),
pp. 28–41.
80
[35] R. McGehee and C. Lehman, A paleoclimate model of ice-albedo feedback forced
by variations in earth’s orbit, SIAM Journal on Applied Dynamical Systems, 11
(2012), pp. 684–707.
[36] R. McGehee and E. Widiasih, A quadratic approximation to budyko’s ice-albedo
feedback model with ice line dynamics, SIAM Journal on Applied Dynamical Sys-
tems, 13 (2014), pp. 518–536.
[37] W. B. McKinnon, D. Prialnik, S. A. Stern, and A. Coradini, Structure
and evolution of kuiper belt objects and dwarf planets, The solar system beyond
Neptune, 1 (2008), pp. 213–241.
[38] A. Nadeau and R. McGehee, A simple formula for a planet’s mean annual
insolation by latitude, Icarus, 291 (2017), pp. 46–50.
[39] E. Noether, Invariante variationsprobleme, Ko¨niglich Gesellschaft der Wis-
senschaften Go¨ttingen Nachrichten Mathematik-physik Klasse, 2 (1918), pp. 235–
267.
[40] G. R. North, Theory of energy-balance climate models, Journal of the Atmo-
spheric Sciences, 32 (1975), pp. 2033–2043.
[41] G. R. North and T. J. Crowley, Application of a seasonal climate model to
cenozoic glaciation, Journal of the Geological Society, 142 (1985), pp. 475–482.
[42] G. R. North, J. G. Mengel, and D. A. Short, Simple energy balance model
resolving the seasons and the continents: Application to the astronomical theory of
the ice ages, Journal of Geophysical Research: Oceans, 88 (1983), pp. 6576–6586.
[43] G. R. North and M. J. Stevens, Energy-balance climate models, Cambridge
University Press Cambridge, 2006.
81
[44] P. E. Olsen, J. Laskar, D. V. Kent, S. T. Kinney, D. J. Reynolds, J. Sha,
and J. H. Whiteside, Mapping solar system chaos with the geological orrery,
Proceedings of the National Academy of Sciences, 116 (2019), pp. 10664–10673.
[45] H. Pa¨like, J. Laskar, and N. J. Shackleton, Geologic constraints on the
chaotic diffusion of the solar system, Geology, 32 (2004), pp. 929–932.
[46] H. Pollard, Mathematical Introduction to Celestial Mechanics, Prentice-Hall
Mathematics Series, Prentice-Hall, 1966.
[47] W. Sellers, A global climatic model based on the energy balance of the earth-
atmosphere system, Journal of Applied Meteorology, 8 (1969), pp. 392–400.
[48] D. A. Short, J. G. Mengel, T. J. Crowley, W. T. Hyde, and G. R.
North, Filtering of milankovitch cycles by earth’s geography, Quaternary Research,
35 (1991), pp. 157–173.
[49] A. B. Sorokovich, The Averaged Hyperbolic Restricted Three-Body Problem, So-
viet Astronomy, 26 (1982), pp. 721–726.
[50] D. Souami and J. Souchay, The solar system’s invariable plane, Astronomy &
Astrophysics, 543 (2012), p. A133.
[51] K. F. Sundman, Recherches sur le proble`me des trois corps, Ex Officina Typo-
graphica Societatis Litterariae Fennicae, 1907.
[52] K. F. Sundman et al., Me´moire sur le proble`me des trois corps, Acta mathe-
matica, 36 (1913), pp. 105–179.
[53] D. E. Trilling, M. Mommert, J. L. Hora, D. Farnocchia, P. Chodas,
J. Giorgini, H. A. Smith, S. Carey, C. M. Lisse, M. Werner, et al., Spitzer
82
observations of interstellar object 1i/‘oumuamua, The Astronomical Journal, 156
(2018), p. 261.
[54] K. K. Tung, Topics in mathematical modeling, Princeton University Press Prince-
ton, NJ, 2007.
[55] J. Walsh, Climate modeling in differential equations, The UMAP Journal, 36
(2015), pp. 325 – 363.
[56] E. Widiasih, Dynamics of the budyko energy balance model, SIAM Journal on
Applied Dynamical Systems, 2 (2013), pp. 2068 – 2092.
[57] P. Wiegert and S. Tremaine, The evolution of long-period comets, Icarus, 137
(1999), pp. 84–121.
[58] D. R. Williams, Jupiter Fact Sheet, 2018 (accessed June 09, 2020). https:
//nssdc.gsfc.nasa.gov/planetary/factsheet/jupiterfact.html.
[59] , Sun Fact Sheet, 2018 (accessed May 17, 2020). https://nssdc.gsfc.nasa.
gov/planetary/factsheet/sunfact.html.
[60] , Earth Fact Sheet, 2019 (accessed May 17, 2020). https://nssdc.gsfc.
nasa.gov/planetary/factsheet/earthfact.htmll.
[61] , Planetary Fact Sheet - Ratio to Earth Values, 2019 (accessed May 17,
2020). https://nssdc.gsfc.nasa.gov/planetary/factsheet/planet_table_
ratio.html.
[62] A. Winter, H. Ishioroshi, T. Watanabe, T. Oba, and J. Christy, Caribbean
sea surface temperatures: Two-to-three degrees cooler than present during the little
ice age, Geophysical Research Letters, 27 (2000), pp. 3365–3368.
83
[63] A. Wintner, The analytical foundations of celestial mechanics, Princeton, NJ,
Princeton university press; London, H. Milford, Oxford university press, 1941.,
(1941).
Appendix A
Mathematica Code
Mathematica was used to create all images and generate all numbers in this document.
The code is given in the section below.
A.1 Regularized HR3BP in 2D
ClearAll["Global‘*"]
SetAttributes[MyExport, HoldAll];
MyExport[picture_, type_] := (SetDirectory[NotebookDirectory[]];
Export[ToString[HoldForm[picture]] <> "." <> ToString[type], picture])
Laplace[{z1_, z2_, w1_, w2_},
m1_] := {(z1^2 - z2^2) - 2*(z1*w2 - z2*w1) (z1*w2 + z2*w1)/m1, (2*z1*z2) +
2*(z1*w2 - z2*w1) (z1*w1 - z2*w2)/m1}/(z1^2 + z2^2)
Eccentricity[{z1_, z2_, w1_, w2_}, m1_] := Module[{a, b},
{a, b} = Laplace[{z1, z2, w1, w2}, m1];
84
85
Sqrt[a^2 + b^2]]
SemiMajorAxis[{z1_, z2_, w1_, w2_}, m1_] :=
With[{e = Eccentricity[{z1, z2, w1, w2}, m1], L = 2*(z1*w2 - z2*w1)},
L^2/(m1*(1 - e^2))]
PerihelionAngle[{z1_, z2_, w1_, w2_}, m1_] := Module[{a, b},
{a, b} = Laplace[{z1, z2, w1, w2}, m1]; ArcTan[a, b]*180/Pi]
Period[{z1_, z2_, w1_, w2_}, m1_ ] :=
With[{A = SemiMajorAxis[{z1, z2, w1, w2}, m1]}, 2*Pi*Sqrt[A^3/m1]]
(*HYPERBOLIC ORBIT*)
a = 1000/5.2;
e = 2;
R[H_] := a*(e*Cosh[H] - 1)
tau[H_] := Sqrt[(e + 1)/(e - 1)] Tanh[H/2]
cos[H_] := (1 - tau[H]^2)/(1 + tau[H]^2)
sin[H_] := 2*tau[H]/(1 + tau[H]^2)
X[H_] := R[H]*cos[H]
Y[H_] := R[H]*sin[H]
ParametricPlot[{X[H], Y[H]}, {H, -2, 2}]
F10 = m2*(x - X[H])/((x - X[H])^2 + (y - Y[H])^2)^1.5 -
m2*X[H]/(X[H]^2 + Y[H]^2)^1.5 // Simplify;
F20 = m2*(y - Y[H])/((x - X[H])^2 + (y - Y[H])^2)^1.5 -
m2*Y[H]/(X[H]^2 + Y[H]^2)^1.5 // Simplify;
86
F1 = F10 /. {m2 -> 0.3550700522909202, x -> 2 (z1[s]^2 - z2[s]^2),
y -> 4*z1[s]*z2[s], H -> H[s]};
F2 = F20 /. {m2 -> 0.3550700522909202, x -> 2 (z1[s]^2 - z2[s]^2),
y -> 4*z1[s]*z2[s], H -> H[s]};
m1 = 0.6449299477090799; (*kg*)
m2 = 1 - m1;
ode = {z1’[s] == w1[s]/2, z2’[s] == w2[s]/2,
w1’[s] == h1[s]*z1[s] - 2*(z1[s]^2 + z2[s]^2)*(F1*z1[s] + F2*z2[s]),
w2’[s] == h1[s]*z2[s] - 2*(z1[s]^2 + z2[s]^2)*(F2*z1[s] - F1*z2[s]),
h1’[s] == -2*F1*(z1[s]*w1[s] - z2[s]*w2[s]) -
2*F2*(z1[s]*w2[s] + z2[s]*w1[s]),
H’[s] == 2*(z1[s]^2 + z2[s]^2)/(e*Cosh[H[s]] - 1)};
(* JUPITER ORBIT *)
smax = 500;
For[i = 0, i < 100, i++,
aJ = 1;
eJ = 0.0489;
z1j = 1/Sqrt[2] Sqrt[aJ (1 + eJ)]*Cos[Pi*i/100]; (*Sqrt[m1/(4*aJ)]*)
z2j = 1/Sqrt[2] Sqrt[aJ (1 - eJ)]*Sin[Pi*i/100];
w1j = -(1/Sqrt[2]) Sqrt[aJ (1 + eJ)] Sqrt[m1/aJ]*Sin[Pi*i/100];
w2j = 1/Sqrt[2] Sqrt[aJ (1 - eJ)] Sqrt[m1/aJ] Cos[Pi*i/100];
87
h1j = (w1j^2 + w2j^2 - m1)/(2*(z1j^2 + z2j^2));
H0j = -3;
icj = {z1[-smax] == z1j, z2[-smax] == z2j, w1[-smax] == w1j,
w2[-smax] == w2j, h1[-smax] == h1j, H[-smax] == H0j};
solj = NDSolve[
ode~Join~icj, {z1[s], z2[s], w1[s], w2[s], h1[s], H[s]}, {s, -smax,
smax}][[1]];
finaleccjup =
Eccentricity[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj /. s -> smax;
finalsmajup =
SemiMajorAxis[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj /. s -> smax;
finallaplacejup =
Laplace[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj /. s -> smax;
finalperijup =
PerihelionAngle[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj /. s -> smax;
perturbedorbitjupiter =
ParametricPlot[{2 (z1[s]^2 - z2[s]^2), 4*z1[s]*z2[s]} /.
solj, {s, -smax, smax},
AxesLabel -> {"q1", "q2"}, PlotStyle -> Gray];
starorbit =
ParametricPlot[{X[H[s]], Y[H[s]]} /. solj, {s, -smax, smax},
AxesLabel -> {"q1", "q2"}, PlotStyle -> Blue ];
icorbitjupiter =
ParametricPlot[{2 ((1/Sqrt[2] Sqrt[aJ (1 + eJ)]*
88
Cos[Sqrt[m1/(4*aJ)]*s])^2 - (1/Sqrt[2] Sqrt[aJ (1 - eJ)]*
Sin[Sqrt[m1/(4*aJ)]*s])^2),
4*(1/Sqrt[2] Sqrt[aJ (1 + eJ)]*Cos[Sqrt[m1/(4*aJ)]*s])*(1/Sqrt[2] Sqrt[
aJ (1 - eJ)]*Sin[Sqrt[m1/(4*aJ)]*s])}, {s, 0, Pi/Sqrt[m1/(4*aJ)]},
AxesLabel -> {"q1", "q2"}, PlotStyle -> {Gray, Dashed}];
finalorbitjupiter =
ParametricPlot[{2 ((1/Sqrt[2] Sqrt[finalsmajup (1 + finaleccjup)]*
Cos[Sqrt[m1/(4*finalsmajup)]*s])^2 - (1/Sqrt[2] Sqrt[
finalsmajup (1 - finaleccjup)]*Sin[Sqrt[m1/(4*finalsmajup)]*s])^2),
4*(1/Sqrt[2] Sqrt[finalsmajup (1 + finaleccjup)]*
Cos[Sqrt[m1/(4*finalsmajup)]*s])*(1/Sqrt[2] Sqrt[
finalsmajup (1 - finaleccjup)]*Sin[Sqrt[m1/(4*finalsmajup)]*s])}, {s,
0, 50}, AxesLabel -> {"q1", "q2"}, PlotStyle -> {Red}];
jupiterorbit =
Show[finalorbitjupiter, icorbitjupiter,
TicksStyle -> {{Medium, Black}, {Medium, Black}},
AxesStyle -> {{Medium, Black}, {Medium, Black}}];
jupandstar =
Show[starorbit, perturbedorbitjupiter, PlotRange -> {{-10, 10}, {-20, 20}},
AxesOrigin -> {0, 0}];
eccicjupiter =
Plot[Eccentricity[{1/Sqrt[2] Sqrt[aJ (1 + eJ)]*
Cos[Sqrt[m1/(4*aJ)]*s], 1/Sqrt[2] Sqrt[aJ (1 - eJ)]*
Sin[Sqrt[m1/(4*aJ)]*s], -(1/Sqrt[2]) Sqrt[
aJ (1 + eJ)] Sqrt[m1/aJ]*Sin[Sqrt[m1/(4*aJ)]*s],
89
1/Sqrt[2] Sqrt[aJ (1 - eJ)] Sqrt[m1/aJ] Cos[Sqrt[m1/(4*aJ)]*s]},
m1], {s, -smax, smax}, PlotRange -> All, AxesLabel -> {s, eccentricity},
PlotStyle -> {Gray, Dashed}];
eccorbjupiter =
Plot[Eccentricity[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj, {s, -smax,
smax}, PlotRange -> All, AxesLabel -> {s, eccentricity},
PlotStyle -> {Red}];
jupiterecc =
Show[eccicjupiter, eccorbjupiter, PlotRange -> {{-smax, smax}, {0, 0.1}},
AxesOrigin -> {0, 0}, TicksStyle -> {{Medium, Black}, {Medium, Black}},
AxesStyle -> {{Medium, Black}, {Medium, Black}}];
smaicjupiter =
Plot[SemiMajorAxis[{1/Sqrt[2] Sqrt[aJ (1 + eJ)]*
Cos[Sqrt[m1/(4*aJ)]*s], 1/Sqrt[2] Sqrt[aJ (1 - eJ)]*
Sin[Sqrt[m1/(4*aJ)]*s], -(1/Sqrt[2]) Sqrt[
aJ (1 + eJ)] Sqrt[m1/aJ]*Sin[Sqrt[m1/(4*aJ)]*s],
1/Sqrt[2] Sqrt[aJ (1 - eJ)] Sqrt[m1/aJ] Cos[Sqrt[m1/(4*aJ)]*s]},
m1], {s, -smax, smax}, PlotRange -> All, AxesLabel -> {s, semimajoraxis},
PlotStyle -> {Gray, Dashed}];
smaorbjupiter =
Plot[SemiMajorAxis[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj, {s, -smax,
smax}, PlotRange -> All, AxesLabel -> {s, semimajoraxis},
PlotStyle -> {Red}];
jupitersma =
Show[smaicjupiter, smaorbjupiter, PlotRange -> {{-smax, smax}, {0.5, 1.4}},
90
AxesOrigin -> {0, 0.5}, TicksStyle -> {{Medium, Black}, {Medium, Black}},
AxesStyle -> {{Medium, Black}, {Medium, Black}}];
laplaceicjupiter =
Plot[Laplace[{1/Sqrt[2] Sqrt[aJ (1 + eJ)]*
Cos[Sqrt[m1/(4*aJ)]*s], 1/Sqrt[2] Sqrt[aJ (1 - eJ)]*
Sin[Sqrt[m1/(4*aJ)]*s], -(1/Sqrt[2]) Sqrt[
aJ (1 + eJ)] Sqrt[m1/aJ]*Sin[Sqrt[m1/(4*aJ)]*s],
1/Sqrt[2] Sqrt[aJ (1 - eJ)] Sqrt[m1/aJ] Cos[Sqrt[m1/(4*aJ)]*s]},
m1], {s, -smax, smax}, PlotRange -> All,
AxesLabel -> {"s", "laplace vector"}, PlotStyle -> {Gray, Dashed}];
laplaceorbjupiter =
Plot[Laplace[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj, {s, -smax, smax},
PlotRange -> All, AxesLabel -> {"s", "laplace vector"}, PlotStyle -> {Red}];
jupiterlaplace =
Show[laplaceicjupiter, laplaceorbjupiter, PlotRange -> All,
AxesOrigin -> {0, 0}, TicksStyle -> {{Medium, Black}, {Medium, Black}},
AxesStyle -> {{Medium, Black}, {Medium, Black}}];
periicjupiter =
Plot[PerihelionAngle[{1/Sqrt[2] Sqrt[aJ (1 + eJ)]*Cos[Sqrt[m1/(4*aJ)]*s],
1/Sqrt[2] Sqrt[aJ (1 - eJ)]*Sin[Sqrt[m1/(4*aJ)]*s], -(1/Sqrt[2]) Sqrt[
aJ (1 + eJ)] Sqrt[m1/aJ]*Sin[Sqrt[m1/(4*aJ)]*s],
1/Sqrt[2] Sqrt[aJ (1 - eJ)] Sqrt[m1/aJ] Cos[Sqrt[m1/(4*aJ)]*s]},
m1], {s, -smax, smax}, PlotRange -> All,
AxesLabel -> {"s", "perihelion angle"}, PlotStyle -> {Gray, Dashed}];
91
periorbjupiter =
Plot[PerihelionAngle[{z1[s], z2[s], w1[s], w2[s]}, m1] /. solj, {s, -smax,
smax}, PlotRange -> All, AxesLabel -> {"s", "perihelion angle"},
PlotStyle -> {Red}];
jupiterperi =
Show[periicjupiter, periorbjupiter, PlotRange -> All, AxesOrigin -> {0, 0},
TicksStyle -> {{Medium, Black}, {Medium, Black}},
AxesStyle -> {{Medium, Black}, {Medium, Black}}];
Print["i = ", i, " ", "(", 2 (z1j^2 - z2j^2), ", ", 4*z1j*z2j, ") ",
" ecc=", finaleccjup, " %sma=", (finalsmajup - aJ)/aJ, " perihelion=",
finalperijup]]
(* EARTH ORBITS *)
smax=500;
For[i = 0, i <16, i++,
aE = 0.192;
eE = 0.0167;
z1e=1/Sqrt[2] Sqrt[aE (1+eE)]*Cos[Pi*i/16]; (*Sqrt[m1/(4*aJ)]*)
z2e=1/Sqrt[2] Sqrt[aE (1-eE)]*Sin[Pi*i/16];
w1e=-(1/Sqrt[2]) Sqrt[aE (1+eE)]Sqrt[m1/aE]*Sin[Pi*i/16];
w2e=1/Sqrt[2] Sqrt[aE (1-eE)]Sqrt[m1/aE]Cos[Pi*i/16];
h1e = (w1e^2+w2e^2-m1)/(2*(z1e^2+z2e^2));
H0e=-3;
ice = {z1[-smax]==z1e,z2[-smax]==z2e,w1[-smax]==w1e,
92
w2[-smax]==w2e,h1[-smax]==h1e,H[-smax]==H0e};
sole =NDSolve[ode~Join~ice,{z1[s],z2[s],w1[s],w2[s],h1[s],H[s]},
{s,-smax,smax}][[1]];
finaleccearth=Eccentricity[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole
/.s->smax;
finalsmaearth=SemiMajorAxis[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole
/.s->smax;
finallaplaceearth = Laplace[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole
/.s->smax;
finalperiearth = PerihelionAngle[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole
/.s->smax;
perturbedorbitearth=
ParametricPlot[{2(z1[s]^2-z2[s]^2),4*z1[s]*z2[s]}/.sole,{s,-smax,smax},
AxesLabel -> {x1[s], x2[s]}, PlotStyle -> Gray];
icorbitearth = ParametricPlot[{2((1/Sqrt[2] Sqrt[aE (1+eE)]*
Cos[Sqrt[m1/(4*aE)]*s])^2 - (1/Sqrt[2] Sqrt[aE (1-eE)]*
Sin[Sqrt[m1/(4*aE)]*s])^2),
4*(1/Sqrt[2] Sqrt[aE (1+eE)]*Cos[Sqrt[m1/(4*aE)]*s])*(1/Sqrt[2]
Sqrt[aE (1-eE)]*Sin[Sqrt[m1/(4*aE)]*s])}, {s, 0,Pi/Sqrt[m1/(4*aJ)]},
AxesLabel -> {"Subscript[q, 1]", "Subscript[q, 2]"},
PlotStyle -> {Gray,Dashed}];
finalorbitearth =
ParametricPlot[{2((1/Sqrt[2] Sqrt[finalsmaearth(1+finaleccearth)]*
Cos[Sqrt[m1/(4*finalsmaearth)]*s])^2 - (1/Sqrt[2]
93
Sqrt[finalsmaearth (1-finaleccearth)]*
Sin[Sqrt[m1/(4*finalsmaearth)]*s])^2),
4*(1/Sqrt[2] Sqrt[finalsmaearth(1+finaleccearth)]*C
os[Sqrt[m1/(4*finalsmaearth)]*s])*(1/Sqrt[2]
Sqrt[finalsmaearth (1-finaleccearth)]*
Sin[Sqrt[m1/(4*finalsmaearth)]*s])},
{s, 0,50}, AxesLabel -> {"q1", "q2"},
PlotStyle -> {Darker[Green]}];
earthorbit=Show[icorbitearth, finalorbitearth,
TicksStyle->{{Medium,Black}, {Medium, Black}}, AxesStyle->
{{Medium, Black}, {Medium, Black}}];
eccicearth = Plot[Eccentricity[{1/Sqrt[2] Sqrt[aE(1+eE)]*
Cos[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE (1-eE)]*
Sin[Sqrt[m1/(4*aE)]*s],-(1/Sqrt[2]) Sqrt[aE
(1+eE)]Sqrt[m1/aE]*Sin[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE
(1-eE)]Sqrt[m1/aE]Cos[Sqrt[m1/(4*aE)]*s]},m1],{s,-smax,smax},
PlotRange->All, AxesLabel -> {s, eccentricity}, PlotStyle -> {Gray, Dashed}];
eccorbearth =
Plot[Eccentricity[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole,{s,-smax,smax},
PlotRange->All, AxesLabel -> {s, eccentricity},
PlotStyle -> {Darker[Green]}];
earthecc=Show[eccicearth,eccorbearth, PlotRange -> {{-smax,smax},{0.01,0.025}},
AxesOrigin -> {0,0.01}, TicksStyle->{{Medium,Black}, {Medium, Black}},
AxesStyle-> {{Medium, Black}, {Medium, Black}}];
94
smaicearth = Plot[SemiMajorAxis[{1/Sqrt[2] Sqrt[aE(1+eE)]*
Cos[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE (1-eE)]*
Sin[Sqrt[m1/(4*aE)]*s],-(1/Sqrt[2]) Sqrt[aE
(1+eE)]Sqrt[m1/aE]*Sin[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE
(1-eE)]Sqrt[m1/aE]Cos[Sqrt[m1/(4*aE)]*s]},m1],{s,-smax,smax},
PlotRange->All, AxesLabel -> {s, semimajoraxis},
PlotStyle -> {Gray, Dashed}];
smaorbearth =
Plot[SemiMajorAxis[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole,{s,-smax,smax},
PlotRange->All, AxesLabel -> {s, semimajoraxis},
PlotStyle -> {Darker[Green]}];
earthsma=Show[smaorbearth,smaicearth, PlotRange ->{{-smax,smax},{0,0.3}},
AxesOrigin -> {0,0}, TicksStyle->{{Medium,Black}, {Medium, Black}},
AxesStyle-> {{Medium, Black}, {Medium, Black}}];
laplaceicearth = Plot[Laplace[{1/Sqrt[2] Sqrt[aE(1+eE)]*
Cos[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE (1-eE)]*
Sin[Sqrt[m1/(4*aE)]*s],-(1/Sqrt[2]) Sqrt[aE
(1+eE)]Sqrt[m1/aE]*Sin[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE
(1-eE)]Sqrt[m1/aE]Cos[Sqrt[m1/(4*aE)]*s]},m11],{s,-smax,smax},
PlotRange->All, AxesLabel -> {"s", "laplace vector"},
PlotStyle -> {Gray, Dashed}];
laplaceorbearth = \
Plot[Laplace[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole,{s,-smax,smax},
PlotRange->All, AxesLabel -> {"s", "laplace vector"},
PlotStyle -> {Darker[Green]}];
95
earthlaplace=Show[laplaceicearth,laplaceorbearth, PlotRange -> All,
AxesOrigin -> {0,0}, TicksStyle->{{Medium,Black}, {Medium, Black}},
AxesStyle-> {{Medium, Black}, {Medium, Black}}];
periicearth = Plot[PerihelionAngle[{1/Sqrt[2]Sqrt[aE(1+eE)]*
Cos[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE(1-eE)]*
Sin[Sqrt[m1/(4*aE)]*s],-(1/Sqrt[2]) Sqrt[aE
(1+eE)]Sqrt[m1/aE]*Sin[Sqrt[m1/(4*aE)]*s],1/Sqrt[2] Sqrt[aE
(1-eE)]Sqrt[m1/aE]Cos[Sqrt[m1/(4*aE)]*s]},m1],{s,-smax,smax},
PlotRange->All, AxesLabel -> {"s", "perihelion angle"},
PlotStyle -> {Gray, Dashed}];
periorbearth =
Plot[PerihelionAngle[{z1[s],z2[s],w1[s],w2[s]},m1]/.sole,
{s,-smax,smax}, PlotRange->All, AxesLabel -> {"s", "perihelion angle"},
PlotStyle -> {Darker[Green]}];
earthperi=Show[periicearth,periorbearth, PlotRange -> All,
AxesOrigin -> {0,0}, TicksStyle->{{Medium,Black}, {Medium, Black}},
AxesStyle-> {{Medium, Black}, {Medium, Black}}];
Print["i = ", i, " ","(",2(z1e^2-z2e^2), ", ",4*z1e*z2e,") ",
earthorbit, " ", finaleccearth, earthecc, finalsmaearth,
earthsma, finalperiearth, earthperi ]]
96
A.2 Budyko-Widiasih Model
Clear["Global‘*"]
SetAttributes[MyExport, HoldAll];
MyExport[picture_, type_] := (SetDirectory[NotebookDirectory[]];
Export[ToString[HoldForm[picture]] <> "." <> ToString[type],
picture])
Q[a_, e_] = 342.998*a^2/Sqrt[1 - e^2];
beta = (23.5/180)*Pi;
sapprox[y_] =
1 - (5/8)*LegendreP[2, Cos[beta]]*LegendreP[2, y] - (9/64)*
LegendreP[4, Cos[beta]]*LegendreP[4, y] - (65/1024)*
LegendreP[6, Cos[beta]]*LegendreP[6, y];
Integrate[sapprox[y], {y, 0, 1}];
sapproxplot =
Plot[sapprox[y], {y, 0, 1}, PlotStyle -> Black, AspectRatio -> 1,
AxesLabel -> {"sin(lat)", "Insolation"}];
sigma[y_, gamma_] :=
Sqrt[1 - (Sqrt[1 - y^2] * Sin[beta] * Cos[gamma] -
y * Cos[beta])^2];
s[y_] := (2/Pi^2) *
NIntegrate[sigma[y, gamma], {gamma, 0, 2 * Pi}];
sfirstprincplot =
Plot[s[y], {y, 0, 1}, PlotStyle -> {Black, Dashed},
AspectRatio -> 1, AxesLabel -> {"sin(lat)", "Insolation"}];
Show[sfirstprincplot, sapproxplot]
97
alpha1 = 0.32; (*land/water*)
alpha2 = 0.62; (*ice*)
alpha[eta_, y_] :=
Piecewise[{{Subscript[alpha, 1],
0 <= y < eta}, {Subscript[alpha,
2], eta < y <=
1}, {(Subscript[alpha, 1] + Subscript[alpha, 2])/2,
y == eta}}]
Manipulate[
Plot[alpha[eta, y], {y, 0, 1}, AspectRatio -> 1,
AxesLabel -> Automatic, PlotStyle -> Black], {eta, 0,
1}]
A = 202;
B = 1.9;
o[T_] := A + B*T
Plot[o[T], {T, 250, 300}, PlotStyle -> {Black}, AspectRatio -> 1,
AxesLabel -> {"T", "OLR"}]
alphabar2[eta_] :=
NIntegrate[Subscript[alpha, 1]*sapprox[y], {y, 0, eta}] +
NIntegrate[Subscript[alpha, 2]*sapprox[y], {y, eta, 1}];
98
alphabar[eta_] =
0.62 (1. - 1.2187411549050502 eta +
0.1644023830328752 eta^3 + 0.07139427615677992 eta^5 -
0.017055504284604825 eta^7) +
0.32 (0. + 1.2187411549050502 eta -
0.1644023830328752 eta^3 - 0.07139427615677992 eta^5 +
0.017055504284604825 eta^7);
Tbarstar[eta_, a_, e_] := (Q[a, e]*(1 - alphabar[eta]) - A)/
B; (*This represents the Global Mean Temperature *)
c = 3.04;
Plot[Tbarstar[eta, 1, 0.0167] ], {eta, 0, 1}, PlotStyle -> Black,
AxesLabel -> {"eta", "GMT"},
TicksStyle -> {{Medium, Black}, {Medium, Black}},
AxesStyle -> {{Medium, Black}, {Medium, Black}},
AxesOrigin -> {0, 0}]
Tstar[eta_, y_, a_,
e_] := (Q[a, e]*sapprox[y] (1 - alpha[eta, y]) - A +
c*Tbarstar[eta, a, e])/(B + c);
Plot[Tstar[Sin[66.5*Pi/180], y, 1, 0.0167], {y,
0, 1}, PlotStyle -> {Black, Thick}]