HD4208 - First results:
Observations and Orbital Parameters of HD4208b:
Motivation and Questions:
The (Continues) Habitable Zone of HD4208:
Within the litterature (Kasting et al. (1993) - Icarus, 101, p. 108) provides detailed calculations for the extend of the habitable zone for main-sequence stars under the assumption of various atmosphere model parameters. The ability of a planet for supporting life, or habitability, is based on our knowledge of life on Earth and on Earth's geometric location within the Solar system. In general, this subject lacks constrains obtained by observations, in order to test the models. However, we will use the results obtained by Kasting et al. as a starting point.
The long-term occurrence of liquid water on a planet is mainly determined by the following key factors: a) the luminosity and surface temperature of the host star and b) geological and atmospheric properties of the terrestial planet. In combination, they determine a unique range from the host star for a planet to sustain liquid water.
Because of the process of stellar evolution, the luminosity of a main sequence star is an increasing function of time. An increasing luminosity results in an incease of total power out put of energy. This again, results in a shift of the habitable zone outwards as the host star evolves along the main sequence within the HR-diagram. Therefore, we define the zero-age-main-sequence (ZAMS) habitable zone (the habitable zone just after the emergence of the star on the main sequence) and the current (now) habitable zone defined by the current age of the host star. The intersection of these regions defines the continues habitable zone, in which planets are habitable throughout the total life-time of the host star.
From a geometrical point of view, the inner and outer radius of the habitable zone constrains the Keplerian orbit of a planet to certain bounds in eccentricity and semi-major axis, by the following conditional relations:
where R_in and R_out are the inner and outer radius of the habitable zone, respectively. The periastron distance is given by a(1-e) and apoastron distance by a(1+e). These relations imposes an constrain on possible allowed values for the eccentricity at a given semi-major axis. The following figure represents the above constraining relationship given R_in(ZAMS) = 0.63 AU, R_out(ZAMS) = 1.26 AU and R_in(now) = 0.74 AU and R_out(now) = 1.34 AU, as determined from Kasting et al. (1993), for the habitable zone within HD4208. Here, a stellar age of 5 x 10^9 years is assumed. As an example for the constrain on the eccentricity: The maximum eccentricity a planet is allowed to obtain for an orbit with semi-major axis a = 1.0 AU, is e_max = 0.26 - in order to remain within the habitable zone. The maximum eccentricity a planet is allowed to reach with a = 0.7 AU, is e_max = 0.1 - which is within the range of the ZAMS habitable zone. The horizontal line at e = 0.2 indicates the upper bound of the outer giant orbit eccentricity.
Numerical Methods:In order to study the time evolution of Keplerian orbital elements (a, e, i, omega, Omega, M) for a given
planet, several numerical solution methods are used to solve the equation of motion for a given set of
Recently, a new class of numerical algorithm, so-called symplectic methods, found a wide-spread use
within numerical celestial mechanics and dynamical astronomy. These methods, as long as the system under
study is "well-behaved", exhibits very good conservation properties of total system energy and angular
momentum and allows the use of a large step size for the numerical discretisation of time within the
integration process. Thus, numerical symplectic methods are fast (by a factor of 10) compared to classic
shemes as the Runge-Kutta or Burlirsch-Stoer methods. In general, symplectic methods are developed using
the principle of "operator-splitting" based on Lie-series, in which the total system Hamiltonian is split
in two components, each measuring the Keplerian perturbation (or indirect two-body perturbation) and the
interaction perturbation (or direct perturbations). Unfortunately, symplectic methods are generelly only
applicable to conservative Hamiltonian dynamical systems. That is, no systems in which energy-dissipation,
due to frictional processes, can be studied using a direct symplectic method. In addition, symplectic methods
are restricted to use a constant step size, allowing to study planetary systems with low-eccentric orbits in
order to capture the correct physics.
However, very recent developements within numerical celestial mechanics introduces numerical symplectic
methods that can follow both non-conservative and high-eccentric dynamical system. The usual problems are
solved by introducing time-transformation and hybrid symplectic methods. Within this study 4 major numerical
methods are used to follow the dynamics of test-particles. Two symplectic and two classic integrators.
The details of constructing symplectic integration methods is a research area of its own, and here we only
state some references.
Dyn. Astro. 67, p. 145 and Mikkola and Palmer (2000) - Cele. Mech. and Dyn. Astron. 77, p. 305).
Astr. Socm 304, p.793) as it is implemented within the mercury6 integration package - MERCURY 6 Integrator .
(Press et al. - Numerical Recipies in Fortran 77). This method is a super-slow but very accurate integrator and
is usually used to test symplectic results in order to check numerical accuracy.
ODEX-integrator by (Hairer, Norsett and Wanner (1993) - Solving Ordinary Differential Equations. Nonstiff
Problems, 2nd edition. Springer Ser. in Comp. Math. vol. 8). See also - Hairers Fortran and Matlab Codes . This
method is used to solve the system equation of motion and the corresponding variational equations, in order
to calculate the MEGNO indicator.
MEGNO-Scans within the Habitable Zone:Since some of the orbital parameters of HD4208b are uncertain and spanning over a given
parameter range, we need to explore the dynamics on test-particles as a function of these
parameters. The system considered is the so called restricted three-body problem. Since the
work of Poincare, it is known that this kind of system do exhibit chaotic behavior. Here, we
denote a chaotic system to show a sensitiv dependence on initial conditions. Within the
litterature, several "chaos-diagnosting" tools exists and can broadly be divded into to
branches: 1) Spectral Fourier Analysis Tools and 2) Lyapunov Tools. Within this study, I use a
Lyapunov derivate chaos indicator - the so-called MEGNO, which stands for mean exponential
growth of nearby orbits. It has been introduced by (Cincotta and Simo, Astr.Astroph.Suppl.
Ser. 147, p.205-228 (2000) and (Cincotta et al., Physica D (2003), 182, p. 151-178) and applied
to problems within dynamical astronomy and celestial mechanics by (K. Gozdziewski et al. (2001) -
A&A 378, p. 569) and (K. Gozdziewski (2002) A&A 393, p.997). However, within the litterature, this
method has never been applied to dynamical analysis problems of stability regions within the
habitable zones of extrasolar planets.
MEGNO is very efficient in distinguishing quasi-periodic and chaotic orbit time evolution on
a timescale that allows a large parameter survey for a given system. Following the theory,
describing and introducing the MEGNO indicator, the MEGNO-value approaches the limit --> 2 for
a quasi-periodic orbit, as time --> infinity. Since the convergence property of MEGNO is fast
in time, the dynamical properties of certain initial conditions can be studied for a large parameter
The following figure, show the result of a (a_EarthPlanet,e_Giant)-parameter scan. Here, a_EarthPlanet
deontes the semi-major axis of an Earth planet (here modelled as a test-particle) and e_Giant is the
outer planets eccentricity. The MEGNO value is color-coded and represented on a log-scale color bar.
Thus, a MEGNO of value 2 (quasi-periodic or regular time evolution) is represented by orange
(log_e(2) = 0.7). The scan covers a (N_x,N_y)=(100,60) grid with a total of 6000 initial conditions
within the habitable zone. All particles have coplanar orbits and started on initially circular orbits.
The semi-major axis and minimum mass of HD4208b are fixed to 1.7 AU and 1.0M_jup, respectively.
For high-eccentric orbits of the outer planet, several lines (indicating chaotic dynamics) become apparent
at 0.8 AU and 0.9 AU. It turns out that these lines represents high-order mean-motion resonances between
the two objects. A very prominent chaotic line, is clearly visible at a_EarthPlanet = 1.06 AU. This line
corresponds to the 2:1 mean-motion resonance. In addition, the 3:2 mean-motion resonance, at a_EarthPlanet =
1.27 AU, acts to be stabilizing through most of the giant planet's eccentricity range - rendering chaotic
orbits for e_Giant = 0.2. Also a "three-finger"-resonance area is apparent around a = 1.2 AU, corresponding to
the 9:5, 17:10 and 8:5 mean motion resonances. Surprising results! This is, where science is fun ... . What
this means regarding planetesimal accretion and subsequent planet formation, is still an open question ...
Resonance Fine Structure:The following figures represents high resolution MEGNO scans around the 2:1 mean motion
resonance line. The scan covers the semi-major axis of the Telluric planet over the
interval [1.0, 1.1]. A grid of (N_x, N_y) = (200, 120) have been used, corresponding to
a total of 24000 initial conditions. The mass of HD4208b is still set to 1 M_jup. The 2:1
resonance line is indicated by a vertical line at 1.06 AU. In addition, the iso-apocenter
contour line for the outer boundary of the ZAMS habitable zone is shown within the figures.
Also visible is the emergence of "overtone"-resonance lines on both the right and left side
of the nominal 2:1-line. These additional resonance lines are much weaker in rendering the
dynamical evolution into a chaotic system. It appears that the relative spacing between these
additional lines follow a fixed numerical relationship.
The last figure within this section is a high resolution of scan of the "three-finger"-resonance
area around a = 1.2 AU. This time the scan only covers an eccentricity range of e_giant within
[0.0; 0.1], since we learned from the initial MEGNO scan that the dynamics is rendered highly
chaotic for giant planet eccentricities larger than 0.1. Again, the single resonances are indicated
by vertical lines and the two iso-apocenter contour lines are superimposed within the figure.
Apparently, the 17:10 line turns out to be a "ghost" or "partial" line, rendering the dynamics
to be quasi-periodic within the line itself!!. Regarding future particle simulations, this result
predicts quasi-periodic or regular orbits of particles perturbed by a giant outer planet with
e_ginat = 0.1 and starting with semi-major axis around 1.28 AU.
Test-particle simulations within the habitable zone:
Initial condition for the test particles are circular orbit with random semi-major axis in the range
a = [0.61, 1.35] AU and M = [0, 360] degrees. A total of 2000 (the very first run used 400) particles
were used for the simulations. The outer giant planet (HD4208b) have e = 0.2 and m = 1 M_jup with
omega = 301 degrees.
Total integration time = 10^6 years :The very first simulation containes 400 test-particles simulated over a time period of 10^6 years. The
second run to be performed, using all CPUs (16) within the local linux cluster (dual- and single intel
CPU linux-boxes) at the AO, considers a total of 2000 particles. For this run, a total time of integration
of 10^8 years were choosen. It turns out that 6 (some machines lasts for 8 days of dedicated CPU time) days
are to long in order to do a proper parameter survey. The data of this run were used to develop several
scripts to extract specific data from the raw data files such as particle survival time, (a,e)- and
(x,y)-movies and (a_0,max(T))-plots.
From the 2000-particle data a (a,e)- and (x,y)-movie (see below) has been produced to show particle
The following jpg-encoded (using mencoder) movies are produced using different frame rates (fps) in order
1-st simulation using 400 particles:
2-nd simulation using 2000 particles:
Total integration time = 150000 years :In order to study the initial (time less than 150000 years) particle dynamics, a 150000 years run have be
performed using 2000 particles with random initial mean anomaly and semi-major axis with each particle
starting on an e = 0 circular orbit. To capture details, the data sampling frequency have been increased and
data batches of the orbital elements (Keplerian and Cartesian) are stored with an output interval of every 200
years. From the data the following animated simulation have been mpeg-encoded using a total of 750 *.jpg frames.
The orbital parameters for the giant planet are (e = 0.2, a = 1.69 AU and omega = 300 deg.) - which could be
considered as a "worst-case scenario" regarding the upper observed giant planets orbit eccentricity. Similar
simulations for different orbit parameters of the outer observed giant planet should be done in order to study the
dependence of particle dynamics as a function of giant planet parameters. Currently, such simulations are in
progress and results will be included and discussed within the thesis. The discussion above reflects much of the
same dynamical features, although in a much greater detail in what is going on ...
HD4208 - habitable or not ?! ...Now, we come to the all encompassing question: Is it possible for a planet (here modelled as a test-particle) to remain bounded
within the habitable zone of HD4208? Since the host star HD4208 is a sun-like star, it evolves along the main sequence within the
HR-diagram. According to stellar evolution models the star's total energy output, or luminosity, is a steady increasing function
of time. Hence, the habitable zone is an annulus with increasing inner and outer radius. Therefore, we define the ZAMS- or t=0 habitable
zone and current habitable zone and the resulting (overlapping) continous habitable zone. At a given time the habitable annulus
constrains the possible orbital parameters (a,e) for a particle given by the periastron- and apoastron relations r_per=a(1-e) and
r_apo=a(1+e). Example: for a given semi-major axis, a particles orbital eccentricity is only allowed to an upper maximum value subject
to the constrain imposed by the model parameters of the habitable zone. Confer (Kasting et al. 1993) for any details of calculating the
habitable region for various assumption on atmosphere and planet parameters. The following zones are calculated using Kasting et al.'s
Within the following figure (or animation) the iso-pericenter contour lines (at 0.63 AU and 0.74 AU) measure the inner radius of the
Furhter Simulations for e = 0.15, 0.10, 0.05 and 0.0 of the observed outer giant planet (with mass 1Mjup):
Summary of current results:The following plot represents a summary of the current results obtained from testparticle simulations using m_HD4208b = 1M_jup as a function
of outer giant planet eccentricity. The maximum eccentricity of testparticles, for a given e-parameter of the outer giant planet, follows
roughly a straight line (except at resonance locations), as a function of test-particle semi-major axis. The lines could formally be
quantified by a linear least-square fit to the maximum eccentricity value at the end of each 150000 yrs simulation. For a circular
(e_HD4208b = 0) orbit, the max-e line almost coincides with e = 0 (again, except at resonance locations) for the test particles and has
therefore been omitted within the plot.
From these short-term simulations, we learn that for a outer giant planet orbit eccentricity of 0.2 would render any planets (particles)
as un-habitable (from a dynamical point of view and assuming only one-planet perturbations), if we consider the region within the
continues habitable zone (and the current HZ). Particle orbits crosses both the pericenter (inner HZ boundary) and apocenter (outer HZ
boundary) contour lines. A periodic crossing of the iso-pericenter line represents the temporal visiting of a planet (particle) to
high-temperature regions. Correspondingly, the crossing of the iso-apocenter line represents temporal ice-ages. In general, the number
of habitable particles is increasing as the outer giant planet's eccentricity is decreasing. This result conforms with our intuition,
since for decreasing orbit eccentricity of the outer giant perturber, the pericenter distance increases with the result of weaker
perturbation action averaged over time.