For computer simulations of dynamical systems one has to choose an appropriate dynamical model, that will give a good approximation for the system under consideration. In a first step we will apply the simplest model which is most probably the fastest one for the computations. This is, in our case, the elliptic restricted three-body problem (ER3BP): where the motion of a mass-less body (m3) is studied in the gravitational field of two massive bodies, i.e., the so-called primaries (mi and m2). Since m3 does not influence the motion of the primaries (i.e., in our case the two stars), they move on Keplerian orbits around their center of mass. This model is commonly used in studies of celestial mechanics and gives quite reasonable results if the mass of the third body is small compared to the other two.
If the third body is quite massive it is advisable to use the three-body problem (TBP), where we are faced with gravitational interactions between all three bodies.
For binaries, where a giant planet has already been detected, one has to use the following models instead of the ER3BP and the TBP, namely:
• the restricted four-body problem (R4BP): where the motion of a mass-less body is studied in the gravitational field of three massive bodies; and
• the four-body problem: where all bodies are massive, so that they interact gravitationally. Its application is also used to verify the results of the R4BP.
For the computations we used two integration methods which are known to work correctly in the case of close encounters, namely the Bulirsch-Stoer method  - which applies the Richardson extrapolation, where a large interval H is divided step by step into n finer subintervals h = H/n - and the Lie integration method [30, 31] - where a special linear differential operator (called the Lie operator) produces a Lie series which is used to solve the equations of motion; the adaptive step-size control of the program allows for the computation of close encounters. With both methods one can obtain high-accuracy solutions for ordinary differential equations.
To determine the state of motion of the computed orbits, one has to use a chaos indicator, on the one hand, or long-term orbital computations and analysis on the other hand.
As a chaos indicator the fast Lyapunov indicator (FLI) (see ) was used, which is quite a fast tool to distinguish between regular and chaotic motion. According to the definition - which is the length of the largest tangent vector:
(n denotes the dimension of the phase space) - it is obvious that chaotic orbits can be found very quickly because of the exponential growth of this vector in the chaotic region. For most chaotic orbits only a small number of primary revolutions is needed to determine the orbital behavior. In order to distinguish between stable and chaotic motion we define a critical value for the FLIs, which is usually set between 107 and 109 depending on the computation time. This method has often been applied to studies of Extrasolar Planetary Systems (see, e.g., [18-20, 22, 23, 33-36] and others).
A simple method used to analyze the orbital parameters of an extensive computation is, e.g., the sup method, where the variation of the sup of an action variable is examined. This method was introduced by Laskar  and later by Froeschle and Lega . For the planetary motion it is useful to apply this method to the eccentricity, since the variation of sup-e determines the orbital behavior quite well - as was shown by Laskar in his long-term studies of our solar system (see ).
Another good characterization of a region in a planetary system is achieved by the max-e, which is the maximum eccentricity of an orbit calculated over the whole integration time. The max-e stability maps are very useful - especially for studies in the so-called habitable zone of a sun-like star (see Dvorak, this volume) - in a sense that one can easily determine the stable regions. This method has been successfully applied to studies of extrasolar planetary systems (see, e.g., [18-20, 36,40-42] and others.)
Was this article helpful?