HD4208 - First results:

Observations and Orbital Parameters of HD4208b:

HD4208 is a Sun like star of spectral type G5V and mass M = 0.93 M_sol. Unfortunately, within the litterature, no stella age has been determined by using spectroscopic measurements. Therefore, an age of 5 x 10^9 years will be assumed in the following discussion. The first announcement of an orbiting planet was given by (Vogt et al., 2001, ApJ, 568, p. 352). They derived the following orbital parameters by fitting a synthetic 2-body Kepler orbit to the observed radial velocity measurements:

  • a = 1.69 AU
  • e = 0.04 +/- 0.12
  • omega = 301 deg (argument of pericenter)
  • P = 829 +/- 36 days and
  • Msin(i) = 0.80 M_jup.

  • Motivation and Questions:

    The very first question to be addressed is, whether or not planets within HD4208 remain confined within the habitable zone and especially within the continues habitable zone? Since HD4208b is a "relative" close-in Jupiter-class giant planet, we can expect strong perturbations acting on any object within the terrestial region of HD4208. What kind of perturbations we will encounter and observe, is a question of dynamically following the time evolution of the Keplerian orbital elements of test-particles, distributed randomly throughout the habitable zone. Our motivation is, of course, the quest and search for a second Earth planet in orbit around a Sun-like star. This kind of theoretical work should give an understanding of conditions whether or not Telluric planets are dynamically capable to remain within the habitable zone of specific extrasolar planetary systems and provide some statistics of the occurrence frequency of Earth-like planets.

    The (Continues) Habitable Zone of HD4208:

    With the general consensus among biologists that carbon-based life requires water for its self-sustaining chemical reactions, the search for habitable planets has therefore focused on identifying environments in which liquid water is stable over billions of years. The habitable zone of a planetary system is defined to be the circumstellar region (or annulus) in which liquid water can exist on a planet with certain atmospheric and geologic properties.

    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:

  • R_in .le. a(1-e) .le. R_out and
  • R_in .le. a(1+e) .le. R_out,

    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
    initial conditions.
    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.

  • Variable-step, time-transformed symplectic integrator developed by (S. Mikkola (1997) - Celes. Mech. and
    Dyn. Astro. 67, p. 145 and Mikkola and Palmer (2000) - Cele. Mech. and Dyn. Astron. 77, p. 305).
  • "Pseudo"-variable step, hybrid mixed variable symplectic integrator (J. E. Chambers (1999) - Mon. Not. R.
    Astr. Socm 304, p.793) as it is implemented within the mercury6 integration package - MERCURY 6 Integrator .
  • Variable-step interpolation method, alias Bulirsch-Stoer-method. This method is well-documentet within
    (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.
  • Variable step interpolation method, known as the Gragg-Bulirsch-Stoer method and is implementet within the
    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 conditions:

    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
    dynamics over a 10^6 year period using 50000 year snapshots. The simulation is showing very interesting
    features of particle dynamics within the habitable zone of HD4208. In particular, regions of mean motion
    resonances display interesting dynamical features. Particles are beeing removed at around 1.25 AU (probably
    due to the 3:2 mean motion resonance located at 1.27 AU). Also the habitable zone is highly dominated by a
    high-eccentricity excitation at the 2:1 mean motion resonance - clearly visible at around a = 1.07 AU.
    Also high-order mean motion resonances (5:2 and 3:1 corresponding to semi-major axis at 0.80 AU and 0.91 AU,
    respectivly) are clearly visible. This is in close agreement with the MEGNO plots (as shown above for
    e_Giant = 0.2), and we can asses these regions to be high-chaotic regions in which mean motion chaotic dynamics
    are causing eccentricity excitation of test particles. An interesting point to remark is, that some high-chaotic
    regions do not have the ability of clearing particles - like mean motion resonances, as discussed above. An example
    for this statement is the region of particles with semi-major axis somewhat greater than 1.1 AU (also dominated by
    the 9:5 mean motion resonance). Within the MEGNO plot, particles starting at this location (and around this location)
    exhibit a high-chaotic time evolution within the systems phase space. From the animated motion in (a,e)-space,
    particles remain confined below an upper eccentricity. Also some "very"-high-order resonances are visible at around
    0.65 AU (best observed within the 400-particle simulation) within the first 5 x 10^4 - 1.5 x 10^5 years, but this is
    not for sure and some detailed study should be done in future runs. In addition, some wave-structure is apparent.
    Whether this means particle transport or something else, is not clearly understood. More runs need to be adressed. The
    time output interval of 50000 years (data sampling frequency) seemed to be to long in order to capture the dynamics in
    detail. Intermediate snapshots of the ongoing dynamics should give more insight. The wave-structure have some kind of
    transient and is dying out at 800000 years or so, leaving the region within the habitable zone completely stirred up
    in a random distribution. After 10^6 years some linear relation is seen between the max-e and the semi-major axis of
    test-mass particles. Regarding the accuracy of the simulation: The total system energy and angular momentum error for
    the simulation are dE/E = 5x10^-09 and dL/L = ~10^-12. Additional runs need to be performed in order study the dynamics
    for the initial time period ...

    The following jpg-encoded (using mencoder) movies are produced using different frame rates (fps) in order
    to control the playback-speed. Thus the simulation is the same, but the speed is increasing from 1 --> 5

    1-st simulation using 400 particles:
  • (a,e)-Movie: Test-Particle Simulation, slow1, 1Myrs, N=400, 0.97 MB
  • (a,e)-Movie: Test-Particle Simulation, slow3, 1Myrs, N=400, 0.97 MB
  • (a,e)-Movie: Test-Particle Simulation, slow5, 1Myrs, N=400, 0.97 MB
  • 2-nd simulation using 2000 particles:
  • (a,e)-Movie: Test-Particle Simulation, slow1, 1Myrs, N=2000, 1.53 MB
  • (a,e)-Movie: Test-Particle Simulation, slow3, 1Myrs, N=2000, 1.53 MB
  • (a,e)-Movie: Test-Particle Simulation, slow5, 1Myrs, N=2000, 1.53 MB
  • 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 ...

  • (a,e)-Movie: Test-Particle Simulation, 150000 years, 15.5 MB
  • (x,y)-Movie: Test-Particle Simulation, 150000 years, crosses, 23.1 MB
  • (x,y)-Movie: Test-Particle Simulation, 150000 years, dots, 6.6 MB
  • 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
    "conservative" estimate.

    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
    habitable zone at time t = 0 years and t = 5 x 10^9 years, respectively. The iso-apocenter contour lines (at 1.26 AU and 1.34 AU) measure
    the outer radius of the habitable zone at the two instance of times. The "(a,e)-triangle" formed by e = 0, iso-pericenter = 0.74 AU and
    iso-apocenter = 1.26 AU is the continouse habitable zone for HD4208 for an assumed age of 5 x 10^9 years (no spectroscopic age
    determination exists within the litterature). The continous HZ is the favorable region for an Earth-like planet to reside in, in order to
    develop life - since here the believed nessecary condition for liquid water is satisfied throughout the host stars lifetime. Follow the
    dynamics by watching the animated movie of the simulation spanning the first 150000 years. The following figure (movie) has the habitable
    zone o-plottet within the same figure to ease the comparison (at the very first time watching the animated simulation, I was stunned and
    needed to rerun it again and again in order to capture each and every detail ...).

    Further simulations, using e_giant = 0.0, 0.05, 0.1 and 0.15 are currently running and will be published in this section.

    Snapshots for particle dynamics (without self-gravitation) within the habitable zone of HD4208
    particle dynamics
    Klick on the image to start download (18.6 MB)

    Furhter Simulations for e = 0.15, 0.10, 0.05 and 0.0 of the observed outer giant planet (with mass 1Mjup):

    e = 0.15 of HD4208b (5.9 MB)

    e = 0.10 of HD4208b (19.3 MB)

    e = 0.05 of HD4208b (6.7 MB)

    e = 0.0 of HD4208b (5.6 MB)

    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.