
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 mainsequence 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 longterm 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 HRdiagram. Therefore, we define the zeroagemainsequence (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 lifetime 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 semimajor 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(1e) and apoastron distance by a(1+e). These relations imposes an constrain on possible allowed values for the eccentricity at a given semimajor 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 semimajor 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 givenplanet, 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, socalled symplectic methods, found a widespread use within numerical celestial mechanics and dynamical astronomy. These methods, as long as the system under study is "wellbehaved", 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 RungeKutta or BurlirschStoer methods. In general, symplectic methods are developed using the principle of "operatorsplitting" based on Lieseries, in which the total system Hamiltonian is split in two components, each measuring the Keplerian perturbation (or indirect twobody 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 energydissipation, 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 loweccentric orbits in order to capture the correct physics. However, very recent developements within numerical celestial mechanics introduces numerical symplectic methods that can follow both nonconservative and higheccentric dynamical system. The usual problems are solved by introducing timetransformation and hybrid symplectic methods. Within this study 4 major numerical methods are used to follow the dynamics of testparticles. 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 superslow but very accurate integrator and is usually used to test symplectic results in order to check numerical accuracy. ODEXintegrator 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. MEGNOScans within the Habitable Zone:Since some of the orbital parameters of HD4208b are uncertain and spanning over a givenparameter range, we need to explore the dynamics on testparticles as a function of these parameters. The system considered is the so called restricted threebody 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 "chaosdiagnosting" 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 socalled MEGNO, which stands for mean exponential growth of nearby orbits. It has been introduced by (Cincotta and Simo, Astr.Astroph.Suppl. Ser. 147, p.205228 (2000) and (Cincotta et al., Physica D (2003), 182, p. 151178) 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 quasiperiodic 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 MEGNOvalue approaches the limit > 2 for a quasiperiodic 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 space. The following figure, show the result of a (a_EarthPlanet,e_Giant)parameter scan. Here, a_EarthPlanet deontes the semimajor axis of an Earth planet (here modelled as a testparticle) and e_Giant is the outer planets eccentricity. The MEGNO value is colorcoded and represented on a logscale color bar. Thus, a MEGNO of value 2 (quasiperiodic 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 semimajor axis and minimum mass of HD4208b are fixed to 1.7 AU and 1.0M_jup, respectively. For higheccentric 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 highorder meanmotion 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 meanmotion resonance. In addition, the 3:2 meanmotion 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 "threefinger"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 motionresonance line. The scan covers the semimajor 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 isoapocenter 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:1line. 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 "threefinger"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 isoapocenter 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 quasiperiodic within the line itself!!. Regarding future particle simulations, this result predicts quasiperiodic or regular orbits of particles perturbed by a giant outer planet with e_ginat = 0.1 and starting with semimajor axis around 1.28 AU.
Testparticle simulations within the habitable zone:Initial conditions:Initial condition for the test particles are circular orbit with random semimajor 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 testparticles simulated over a time period of 10^6 years. Thesecond run to be performed, using all CPUs (16) within the local linux cluster (dual and single intel CPU linuxboxes) 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 2000particle data a (a,e) and (x,y)movie (see below) has been produced to show particle
The following jpgencoded (using mencoder) movies are produced using different frame rates (fps) in order 1st simulation using 400 particles:2nd 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 beperformed using 2000 particles with random initial mean anomaly and semimajor 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 mpegencoded 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 "worstcase 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 testparticle) to remain boundedwithin the habitable zone of HD4208? Since the host star HD4208 is a sunlike star, it evolves along the main sequence within the HRdiagram. 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(1e) and r_apo=a(1+e). Example: for a given semimajor 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 isopericenter 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 functionof outer giant planet eccentricity. The maximum eccentricity of testparticles, for a given eparameter of the outer giant planet, follows roughly a straight line (except at resonance locations), as a function of testparticle semimajor axis. The lines could formally be quantified by a linear leastsquare fit to the maximum eccentricity value at the end of each 150000 yrs simulation. For a circular (e_HD4208b = 0) orbit, the maxe 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 shortterm simulations, we learn that for a outer giant planet orbit eccentricity of 0.2 would render any planets (particles) as unhabitable (from a dynamical point of view and assuming only oneplanet 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 isopericenter line represents the temporal visiting of a planet (particle) to hightemperature regions. Correspondingly, the crossing of the isoapocenter line represents temporal iceages. 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.
