## Abstract

A special symmetric configuration of the general four-body problem called the Caledonian Symmetric Four-Body Problem is investigated. Recently, Steves & Roy have discovered a global stability parameter for the system, called the Szebehely constant. The connection between the chaotic behaviour of the phase space and the global stability given by the Szebehely constant is analysed by using the relative Lyapunov indicator (RLI) and smallest alignment indices (SALI) methods of determining regular and chaotic motion. It is found that as the Szebehely constant is increased, making the system hierarchically stable from a global point of view, the corresponding phase space becomes increasingly more regular and less chaotic.

## Introduction

Approximately two thirds of the stars in our Galaxy exist as part of multistellar systems. Around one fifth of these are part of triple systems, while a rough estimate suggests that a further one fifth of these triples belong to quadruple or higher systems, which can be modelled by the four-body problem. The four-body problem is also increasingly being used to study and explain many of the complex dynamical phenomena found in the Solar system and exoplanetary systems. Studying the stability and its relation to chaotic and regular regions of motion in four-body systems is therefore fundamental to understanding the evolution of quadruple stellar and exoplanetary systems.

It has long been known that the three or more body problem is non-integrable; however information on the global stability of a gravitational system can be found by studying the regions of the phase space where real motion is possible. In the general three-body problem, Zare (1976, 1977) and Marchal & Saari (1975) showed that the energy *H* and angular momentum *c* integrals can combine to form a stability criterion *c*^{2}*H* which, for a given critical value *c*^{2}*H*_{crit}, produces a phase space for the three-body system which contains topologically separate sub-regions. Thus the hierarchical arrangement of the system which exists in one sub-region of the real motion cannot physically evolve into the hierarchical arrangement of the system existing in another sub-region of the real motion. Or in other words, the closest two bodies which form a binary pair perturbed by the third outer body cannot rearrange their order, to form a closer binary with the third body. In this situation the system can be said to be hierarchically stable for all time.

The large number of degrees of freedom in the four or more few-body problems make stability criteria of this form difficult to formulate. Recently, Roy & Steves (2001) have developed the Caledonian Symmetric Four-Body Problem (CSFBP), a symmetrical four-body problem which utilizes all possible symmetries so that it can easily be studied in three-dimensional space. In this special symmetric problem, they have discovered a stability criterion (Steves & Roy 2001), based on a similar combination of the energy *E* and angular momentum *c* integrals, which they called the Szebehely constant, *C*_{0}. They showed that the hierarchical stability and possible hierarchical evolutions of the symmetric four-body system can be predicted solely from knowledge of the total energy and angular momentum of the system.

The recent development of fast and efficient methods of detecting chaotic behaviour in gravitational dynamical systems such as the method of the relative Lyapunov indicator (RLI) (Sándor, Érdi & Efthymiopoulos 2000) and the method of the smallest alignment indices (SALI) (Skokos 2001) is enabling detailed pictures of the chaotic behaviour of systems to be studied.

In this paper, we adapt the two chaos-detecting methods RLI and SALI to the coplanar CSFBP and investigate the chaotic and regular structure of the phase space. By studying the variation of the phase space for a variety of Szebehely constants, we can determine the relationship between the global hierarchical stability of the the four-body system and the chaotic/regular regions of motion.

## Description of The Coplanar CSFBP

The CSFBP is formulated by using all possible symmetries. The main feature of the model is its use of two types of symmetry: (i) past–future symmetry and (ii) dynamical symmetry. Past–future symmetry exists in an *N*-body system when the dynamical evolution of the system after *t*= 0 is a mirror image of the dynamical evolution of the system before *t*= 0. It occurs whenever the system passes through a mirror configuration, i.e. a configuration in which the velocity vectors of all the bodies are perpendicular to all the position vectors with respect to the centre of mass of the system (Roy & Ovenden 1955). Dynamical symmetry exists when the dynamical evolution of two bodies on one side of the centre of mass of the system is paralleled by the dynamical evolution of the two bodies on the other side of the centre of mass of the system. The resulting configuration is always a parallelogram, but of varying length, width and orientation.

The CSFBP is three-dimensional and involves initially two pairs of bodies, one pair on either side of the centre of mass, with each pair having components of unequal masses, but the same two mass values as for the other pair. To set up the model we make the following assumptions.

All four bodies are finite point masses, with two bodies

*P*_{1}and*P*_{4}on opposite sides of the centre of mass of the system*O*having mass*m*and the other two bodies*P*_{2}and*P*_{3}having mass*M*. We define μ=*m*/*M*so that 0 < μ≤ 1.At

*t*= 0, the bodies are collinear and their velocity vectors perpendicular to the line the bodies lie on. This ensures past–future symmetry.At

*t*= 0, the radius and velocity vectors of the bodies with respect to the four-body centre of mass of the system*O*are given by*r*_{4}=−*r*_{1};*r*_{3}=−*r*_{2};*V*_{4}=−*V*_{1};*V*_{3}=−*V*_{2}. Note that*V*_{1}and*V*_{4}do not need to be coplanar with*V*_{2}and*V*_{3}. This ensures that dynamical symmetry holds for all time and the configuration of the four bodies always forms a parallelogram.

In this paper, we will simplify the problem to focus solely on the coplanar CSFBP where the radius and velocity vectors are coplanar. See Fig. 1 for the initial configuration of the coplanar CSFBP and Fig. 2 for a possible configuration after time *t*.

Let the centre of mass of the system be at rest and at the origin of the coordinate system. Given the dynamical symmetry conditions of (c) valid at any time *t*, the equations of motion of the CSFBP simplify to:

*m*

_{1}=

*m*is the mass of

*P*

_{1}and

*m*

_{2}=

*M*is the mass of the more massive body

*P*

_{2}. The force function

*U*can be written as: where

*r*is the length of the radius vector of

_{i}*P*from the centre of mass for

_{i}*i*= 1, 2 and

*r*

_{12}is the distance between

*P*

_{1}and

*P*

_{2}.

Let *E*_{0}=−*E*. The moment of inertia *I* may be written as

*I*, the kinetic energy

*T*and the magnitude of the angular momentum vector

*c*, becomes:

For a given energy and angular momentum, equation (6) gives boundary surfaces in the *r*_{1}*r*_{2}*r*_{12} space which define the regions of real motion.

From the kinematic constraints of the problem, we also have the limitation

Let us now introduce dimensionless variables ρ_{1}, ρ

_{2}, and ρ

_{12}, and a new dimensionless constant

*C*

_{0}called the Szebehely constant: In terms of these variables, Sundman's inequality (6) becomes and the kinematic constraint becomes

Four types of collisions are possible in this problem.

‘12’ collision (double binary (DB) collision). If ρ

_{12}= 0 then*P*_{1}collides with*P*_{2}and*P*_{3}collides with*P*_{4}. The inequalities (10) are satisfied only if ρ_{2}=ρ_{1}.‘13’ collision (DB collision). If , then

*P*_{1}collides with*P*_{3}and*P*_{2}collides with*P*_{4}.‘23’ collision (single binary (SB) collision). If ρ

_{2}= 0 then*P*_{2}collides with*P*_{3}. The inequalities are satisfied only if ρ_{12}=ρ_{1}.‘14’ collision (SB collision). If ρ

_{1}= 0 then*P*_{1}collides with*P*_{4}. The inequalities (10) are satisfied only if ρ_{12}=ρ_{2}.

## Global Hierarchical Stability Criteria in The CSFBP

Setting the angular momentum to zero in equation (6) or the Szebehely constant to zero in equation (9) is equivalent to defining the surfaces of zero-velocity for the problem. The zero-velocity surfaces have been studied in detail in Roy & Steves (2001).

It can be shown that real motion can take place in four tube-like regions of the ρ_{1}ρ_{2}ρ_{12} space that are connected to each other near the origin for *C*_{0}= 0. Each tube represents a particular hierarchical arrangement of the bodies as follows (see Fig. 3).

‘12’-type hierarchy (Fig. 3a). The *P*_{1}*P*_{2} bodies and the *P*_{3}*P*_{4} bodies form binaries which orbit around the baricenter of the four-body system. *r*_{12} is small, while *r*_{1} and *r*_{2} are of similar size. This is a double binary (DB) system.

‘13’-type hierarchy (Fig. 3b). The *P*_{1}*P*_{3} bodies and *P*_{2}*P*_{4} bodies form binaries which orbit around the baricenter of the four-body system. *r*_{12} is large, while *r*_{1} and *r*_{2} are of similar size. This is a DB system.

‘23’-type hierarchy (Fig. 3c). The *P*_{2}*P*_{3} bodies form a central binary and the *P*_{1} and *P*_{4} bodies orbit around it. *r*_{2} is small, while *r*_{1} is large. This is a single binary (SB) system.

‘14’-type hierarchy (Fig. 3d). The *P*_{1}*P*_{4} bodies form a central binary and the *P*_{2} and *P*_{3} bodies orbit around it. *r*_{1} is small, while *r*_{2} is large. This is a SB system.

For *C*_{0}= 0 the bodies are free to travel from one region to another and thus it is possible for the hierarchical arrangement of the system to change to any of the four possibilities.

Combining equation (9) with the kinematic constraint of equation (10) which gives the extreme values of ρ_{12} for a given ρ_{1} and ρ_{2}, it is possible to capture information on the connectivity of the possible regions of motion in the ρ_{1}ρ_{2}ρ_{12} space by projecting the resulting boundary surfaces on to the ρ_{1}ρ_{2} plane. Examples of such projections can be seen in Fig. 4. Real motion is possible only in the regions shown without dots. For *C*_{0}= 0 the projection shows the allowed regions of motion are connected. The ρ_{2}≪ρ_{1} arm represents the ‘23’-type hierarchy, the ρ_{1}≈ρ_{2} arm represents both the ‘12’ and ‘13’-type hierarchies and the very narrow ρ_{1}≪ρ_{2} arm represents the ‘14’-type hierarchy. The connected projection indicates that change from one hierarchical arrangement to another is possible.

When the Szebehely constant is allowed to increase from *C*_{0}= 0 in equation (9), projections of the boundary surfaces on to the ρ_{1}ρ_{2} plane show regions near the origin where motion is not allowed. The greater the value of *C*_{0}, the more extended these forbidden regions become, until for a critical value of *C*_{0}=*C*_{crit 1} the lower allowed region, where ρ_{2}≪ρ_{1}, becomes totally disconnected from the other allowed regions. At this point, any ‘23’ hierarchy characterizing this lower allowed region will remain a ‘23’ hierarchy forever and any other type of hierarchy cannot evolve into a ‘23’ hierarchy. When the value of *C*_{0} reaches a second critical value of *C*_{0}=*C*_{crit 2}, all the allowed regions are totally disconnected and evolution from one hierarchy to another is impossible, ensuring that there is hierarchical stability for all time. Detailed derivation of these results can be found in Steves & Roy (2001).

The two critical values of *C*_{0} which give hierarchical stability first for the ‘23’ hierarchy and then for all four hierarchies for all time are functions only of the mass ratio μ. For example for μ= 1, *C*_{crit 1}=*C*_{crit 2}= 46.841 and for μ= 0.1, *C*_{crit 1}= 0.7792 and *C*_{crit 2}= 0.8886. In this paper we shall study these two μ values as representatives. μ= 1 corresponds to a quadruple stellar system. In the μ= 0.1 case, the CSFBP models a system with two very massive planets or brown dwarfs in a binary star system.

## Equations of Motion and Initial Conditions for The Coplanar CSFBP

### The equations of motion

It is convenient to change the position variables from *r*_{1}, *r*_{2} and *r*_{12} to *r*_{1}, *r*_{2} and α, the angle formed by *P*_{1}*OP*_{2}. The corresponding conjugate momenta are *p*_{1}, *p*_{2} and *p*_{α}.

By utilizing the symmetries and reducing the general equations of motion the Hamiltonian for the system can be written as (Széll 2003):

where*c*is the constant angular momentum.

The equations of motion can then be written as:

Following the description of the initial configuration as given in Section 2, the initial values of the canonical variables are given as follows: α= 0 so that initially all bodies are collinear; and *p*_{1}=*p*_{2}= 0 so that *V*_{1} and *V*_{2} are perpendicular to the line the bodies lie on.

We will vary the initial *r*_{1} and *r*_{2} variables with fixed energy (*E*_{0}) and Szebehely constant (*C*_{0}) values.

The initial value of *p*_{α} and the angular momentum constant *c* are determined from *E*_{0} and *C*_{0}, through equations (8) and (11). We recall that *H*=−*E*_{0}.

## RLI and Sali Methods of Chaos Detection

A traditional way to detect chaotic behaviour of orbits in dynamical systems is the calculation of the maximum Lyapunov characteristic exponent (LCE). The LCE of a trajectory emanating from an initial point ** x*** of the phase space is defined as the limit

*ξ*_{t}is the image of an initial infinitesimal small deviation vector

*ξ*_{0}after time

*t*. The evolution of

*ξ*_{t}in the tangent space can be calculated by numerical integration of the equations of motion together with their linearized equations: where

*D*

**is the Jacobian matrix evaluated at**

*f(x)***. If**

*x**L*

^{1}(

***) = 0, the orbit evolving from**

*x**** is regular. If**

*x**L*

^{1}(

***) > 0 it is chaotic. A serious disadvantage of the calculation of the LCE is that it can not be obtained after finite integration time. Thus its value can only be extrapolated, which makes the identification of weakly chaotic orbits very uncertain. Furthermore, weak chaos plays an essential role in understanding the long-term stability of the Solar system and other exoplanetary systems. In recent years there has been a growing interest in the development and application of fast chaos detection methods. Being aware of the incompleteness of the references below, we refer to the FLI method of Froeschlé, Lega & Gonczi (1997), to the method of spectral distance by Voglis, Contopoulos & Efthymiopoulos (1998), and to the SALI method of Skokos (2001). In Sándor et al. (2000) we introduced the concept of the relative Lyapunov indicator (RLI), which has also proven in our investigations to be an efficient tool of chaos detection. Before the definition of the RLI, we recall the definition of the finite-time Lyapunov indicator originating from an initial point**

*x*

*x*_{0}: The RLI method is based on the fact that the finite-time approximation of the Lyapunov indicator of two neighbouring orbits,

*L*(

***,**

*x**t*) and

*L*(

***+δ**

*x****,**

*x**t*) evolves similarly (as a function of

*t*) for regular orbits, and differently for chaotic orbits. In order to quantify the time-evolution of

*L*(

**,**

*x**t*) for neighbouring regular and chaotic orbits, we introduced the idea of the relative Lyapunov indicator (Sándor et al. 2000): where δ

***≪ 1. Our investigations on different dynamical systems such as the restricted three-body problem (Sándor et al. 2000); the 2D, 4D, and 6D symplectic mappings (yet unpublished), and the stability of certain exoplanetary systems (Pál & Sándor 2002) show that the curve Δ**

*x**L*(

***,**

*x**t*) (as a function of

*t*) exhibits typical behaviour for regular and chaotic orbits, which differ essentially from each other.

In Figs 5(a), (b) and (c), the behaviour of the curves Δ*L*(** x***,

*t*) (as a function of time) can be seen for regular, weakly chaotic and strongly chaotic orbits of the CSFBP. The curves separate well in magnitude, with the strongly chaotic and regular orbits having RLI values of 10

^{−6}and 10

^{−8}, respectively. Of particular note is that these distinguishing values are obtained very quickly, even after 1000 iterations or time-steps Δ

*t*.

Because the chaotic behaviour of orbits can be detected after a relatively short time of numerical integration, the method of the relative Lyapunov indicators enables us to study a large set of initial conditions, and thus allows us to explore the regular and chaotic regions of the phase space. In order to separate regular and chaotic regions of the phase space, we calculate the average value of Δ*L*(** x***,

*t*) for a given integration time

*t**:

*** is in a regular region, 〈Δ**

*x**L*(

***)〉**

*x*_{t*}is small. Otherwise, if

*** is in a chaotic region, 〈Δ**

*x**L*(

***)〉**

*x*_{t*}will be larger.

We compared the RLI method with the SALI method by Skokos (2001). The SALI method is based on the calculation of the evolution of two initial infinitesimal small deviation vectors ξ_{0} and η_{0}. The parallel alignment index is defined as:

*d*

_{−}= 0 implies that the deviation vectors coincide. If

*d*

_{+}= 0, then the initial vectors become opposite to each other. (The deviation vectors are normalized to 1 at each step of the calculation.) It is true also that

*d*

_{+},

*d*

_{−}∈[0, 2]. According to Skokos (2001), in the case of a chaotic orbit (in at least a 4D symplectic mapping) the two deviation vectors tend to coincide. Thus

*d*

_{−}→ 0 and

*d*

_{+}→ 2. In the case of a regular orbit, both alignment indices tend to positive values in the interval [0, 2]. Summarizing the behaviour of the alignment indices: in the chaotic case, the smaller alignment index (SALI) tends to zero, while in the regular case, it tends to a non-zero value. The different behaviour of the SALI number for ordered and chaotic orbits makes it an appreciable tool in detecting ordered and chaotic orbits of dynamical systems. In Figs 6(a), (b) and (c), the behaviour of the SALI number (as a function of time) can be seen for ordered, weakly chaotic, and strongly chaotic orbits of the CSFBP with the same initial conditions as given in the cases shown in Fig. 5.

## Results

In order to determine the relationship between the global hierarchical stability of the CSFBP and the chaotic and regular nature of its phase space, we studied two indicative CSFBP systems stipulated by μ= 1 (a quadruple stellar system) and μ= 0.1 (a binary star system with two very massive planets or brown dwarf stars). For each given mass ratio μ, we kept the energy of the system constant, effectively providing a relative scale of size for the system and then chose a range of increasing Szebehely constants *C*_{0} which covered the complete range of different possibilities for global hierarchical stability. i.e. *C*_{0} < *C*_{crit 1} (unstable); *C*_{crit 1} < *C*_{0} < *C*_{crit 2} (‘23’ hierarchy stable) and *C*_{crit 2} < *C*_{0} (stable). The values of μ, *E*_{0} and *C*_{0} chosen are given in Table 1.

For any given μ, *E*_{0} and *C*_{0} each CSFBP orbit is uniquely defined by its initial values *r*_{1} and *r*_{2}. The nature of each orbit can therefore be depicted in *r*_{1}*r*_{2} space by colour code. Using this fact we can construct *r*_{1}*r*_{2} graphs colouring each point black when it belongs to a chaotic orbit, light grey when it belongs to a regular orbit, and grey-scaled from light grey to black for orbits in between.

### Method

We chose a fine grid containing of the order of 150 000 pairs of initial values of *r*_{1} and *r*_{2} covering all regions of real motion in the *r*_{1}*r*_{2} plane. We numerically calculated the RLI and SALI numbers for each CSFBP orbit, using 10 000 iterations for each case.

The initial separation in the RLI method was chosen to be 10^{−11}. The chaotic nature of each orbit was categorized for both methods of fast chaos detection. See Table 2 for the categories and colour code.

Where the conservation of energy failed during the integration, the orbit was categorized as one of the four types of collisions according to the criteria given at the end of Section 2. The different behaviours of the orbits were plotted in the *r*_{1}*r*_{2} plane, where each point on the plane represents the starting conditions of an orbit. Each point in the *r*_{1}*r*_{2} plane is colour coded according to the nature of the orbits as defined in Table 2.

For each μ–*C*_{0}-*E*_{0} triple we produce two *r*_{1}*r*_{2} graphs: one for each of the methods of detecting chaotic RLI and SALI behaviour. The results are presented in the following sections.

* 6.1.1 The nature of the orbits in the initial r*_{1}r_{2} space for*μ= 1.0*

_{1}r

_{2}space for

In this case *E*_{0} was fixed to be 7 and *C*_{0} was chosen to be 10, 40 and 60. Recall that *C*_{crit 1}=*C*_{crit 2}= 46.841. The initial values of *r*_{1} and *r*_{2} were varied from 0 to 4 with a step size of 0.01. We used units such that *GM*^{2}= 1.

In the case where all bodies have equal mass, an interchange of the *r*_{1} and *r*_{2} produces the same physical orbit. Thus the categorization of the orbits is symmetric with respect to the *r*_{1}=*r*_{2} line with the exception that the ‘23’ and ‘14’ collisions are reversed. That is, for *r*_{1}≫*r*_{2} the initial ordering of the equal masses is 1234. For *r*_{2}≫*r*_{1} the initial ordering of the equal masses is 2143. These initial conditions represent the same physical orbit since the masses are equal. The resulting categorization will therefore be the same except that the type of collisions will be reordered according to the different order. Thus comparing the collision type from *r*_{1}≫*r*_{2} space to *r*_{2}≫*r*_{1} space; ‘23’ collisions (blue) become ‘14’ collisions (yellow) and vice versa. ‘12’ (red) and ‘13’ (green) collisions remain the same. We therefore need only to integrate half the *r*_{1}*r*_{2} space eg we integrated 4 × 100 × 4 × 100/2 = 80 000 orbits. In order to be able to see the detail near the origin where the most important features of the phase space can be found, we have plotted the graphs up to *r*_{1}=*r*_{2}= 1: see Fig. 7. The top line of the three graphs shows the RLI categorization of the orbits for *C*_{0}= 10, 40 and 60, while the bottom line shows the SALI categorization.

The graphs can be separated into two regions: the double binary (DB) region around *r*_{2}≈*r*_{1} and the two single binary regions around *r*_{2}≪*r*_{1} (SB1) and *r*_{1}≪*r*_{2} (SB2). In the DB region there always exists two massive regular regions, one on either side of the *r*_{2}=*r*_{1} line, independent of the *C*_{0} value. For small *C*_{0} values, the SB regions contain primarily collision orbits which indicates chaotic behaviour. As *C*_{0} increases the number of the collision orbits decreases. For all values of *C*_{0}, collision orbits exist near the boundaries of real motion. Collisions are more likely to occur near the boundaries of real motion, because at these boundaries by definition at least one pair of bodies has collided. The expected type of collision will be the type associated with the hierarchy that the system initially starts in, as such a configuration is more likely to lead to collision before it can change its hierarchy. Thus in the DB region the collision type near the boundaries is ‘12’ collision or red, as all orbits in this region begin in a ‘12’-type hierarchy. The SB2 region has ‘14’ collisions or yellow near its boundaries and the SB1 region has ‘23’ collisions or blue near its boundaries. In the following we analyse each of the graphs in detail.

*C*_{0}= 10 (Fig. 7, left-hand panel). In this case *C*_{0} is well below the critical value of *C*_{crit 1}=*C*_{crit 2}.

The regular part of the DB region is surrounded by a relatively thick chaotic region. The chaotic behaviour is characterized by the existence of all four kinds of collisions, denoted red, green, blue and yellow. The most outer part of the DB region is red, indicating the existence of ‘12’-type collisions near the boundaries of real motion. In the RLI case, there is a sizable dark grey region of transition between the ordered and chaotic regions. We recall that the RLI values here are between 10^{−8} and 10^{−6}. In the SALI graph there is a much smaller dark grey region (i.e. SALI values between 0.01 and 0.00001), indicating that in the SALI case the chaotic and regular regions are more clearly delineated. The most marked difference between the SALI and RLI graphs can be found close to the origin. In the SALI graph this area appears to be regular but in the RLI graph it is primarily a transition region with a small chaotic area. Thus the regularity of this region is doubtful.

The SB regions are well separated from the DB region by a region of ‘12’ collisions (red). It begins at around the line from point (0.5, 0) to (0.8, 0.3). This area contains only collision orbits indicating chaotic behaviour. A very interesting pattern can be recognized here. The green (‘13’ collision) and red (‘12’ collision) colours which dominate both SB regions alternate, while blue (‘23’ collision) and yellow (‘14’ collision) regions can be found between them. The whole SB1 region is surrounded by blue ‘23’ collisions, while SB2 is by yellow ‘14’ collisions.

It can be concluded that the picture is very chaotic except for the central region of the DB area. The collision orbits dominate the phase space. Real systems are likely to be found only in the central regular DB region.

*C*_{0}= 40 (Fig. 7, middle panel). In this case *C*_{0} is just below the critical value of *C*_{crit 1}=*C*_{crit 2}.

In the SB regions, islands of chaotic and regular behaviour appear. Most of the collisions can now be found nearer the DB area. For the smaller *r*_{1} values in the SB1 region, the islands are not well determined; they are dark grey in the RLI graph. For large *r*_{1} values, the orbits are more regular; they are light grey in both the RLI and SALI graphs.

In the SB1 region, real physical orbits are therefore likely to be found for greater *r*_{1} values and similarly in the SB2 region for greater *r*_{2} values.

*C*_{0}= 60 (Fig. 7, right-hand panel). In this case *C*_{0} is much greater than the critical value *C*_{crit 1}=*C*_{crit 2}.

The SB and DB regions are now disconnected. The phase space is regular except close to the boundaries of real motion where chaotic motion and collisions exist. The DB, SB1 and SB2 regions are surrounded by ‘12’-, ‘23’- and ‘14’-type collisions, respectively. Only one type of collision is seen in each disconnected region, corresponding to the one type of hierarchy state that is possible in that region.

The pictures show very regular behaviours.

*The nature of the orbits in the initial r*_{1}r_{2} space for*μ= 0.1*

_{1}r

_{2}space for

In this case *E*_{0} was fixed to be 1.2 and *C*_{0} was chosen to be 0.4, 0.8 and 0.9. Recall that for μ= 0.1 *C*_{crit 1}= 0.7792 and *C*_{crit 2}= 0.8886. *r*_{1} and *r*_{2} were varied from 0 to 3.34 with a step size of 0.0083.

Now the phase space is not invariant to the interchanging of the *r*_{1} and *r*_{2} values. Thus we need to integrate the whole range of the initial parameters. In this case we integrated orbits. As in the previous case the most important features of the phase space can be found close to the origin. Hence we have plotted the graphs up to *r*_{1}= 2.2 and *r*_{2}= 1.5 in order to be able to see more of the detail: see Fig. 8. As in Fig. 7 for each *C*_{0} value, we have two graphs showing the categorization of the orbits according to the RLI and SALI numbers.

The graphs can be separated into four regions: two double binary regions around *r*_{2}≈*r*_{1} but *r*_{2} < *r*_{1} (DB1) and *r*_{2} > *r*_{1} (DB2), and two single binary regions, *r*_{2}≪*r*_{1} (SB1) and *r*_{1}≪*r*_{2} (SB2). In the DB1 region there always exists a massive regular region just below the *r*_{2}=*r*_{1} line, independent of the *C*_{0} value. For small *C*_{0} values, the other regions contain almost solely collision orbits which indicates chaotic behaviour. The SB2 region is very narrow compared to the SB1 region. As *C*_{0} increases, the number of collision and chaotic orbits decreases. In the following we analyse each graph in detail.

*C*_{0}= 0.4 (Fig. 8, left-hand panel). In this case *C*_{0} is well below the critical value *C*_{crit 1}.

The SB2 region is too narrow to have much structure at this scale. ‘14’-type (yellow) collision orbits dominate. Some chaotic (black) orbits can be seen at around *r*_{2}≈ 0.9.

In the DB2 region the phase space has a very complicated structure. The large number of collision orbits suggests chaotic behaviour. The ‘12’- (red) and ‘13’- (green) type binary–binary collisions form a regular pattern. Some ‘14’ (yellow) collisions can be found between them. This dominance of collisions indicates chaotic behaviour. Between the collision regions, islands of indeterminate behaviour can be seen, i.e. they are light grey (regular) in the SALI graph and dark grey (undetermined) in the RLI graph. The centres of these islands are light grey in the RLI graphs, and therefore these orbits are likely to be regular.

The largest non-collision island can be found at (0.35, 0.4) and it is mostly regular. Slightly below, at (0.1, 0.2), another non-collision island can be found. In the RLI graph this region is black though it has some light grey regions, indicating that regular orbits exist there. This island is special in that it is the only area where the SALI method shows regular behaviour and the RLI shows chaotic behaviour. In the future we will investigate this island in detail.

The DB1 region has very similar behaviour to the μ= 1 case. It contains a massive regular region with a chaotic and collision region at the boundaries of real motion.

The SB1 region has a very spectacular pattern. Clear big non-collision regular islands (light grey) can be found between the predominantly ‘13’ (green) and ‘12’ (red) collision orbits in the RLI graph. Along the upper and lower boundaries of the SB1 region (i.e. *r*_{2}= 0 and *r*_{2}= 0.5), there exists islands of undetermined behaviour, dark grey in the RLI graph and light grey in the SALI graph. For large initial *r*_{1} values, the RLI shows these regions to be regular (light grey) as well.

It can be concluded that the phase space is quite chaotic with some regular islands. Real physical systems are most likely to be found in the regular regions of the DB1 and SB1 areas with some in the fewer regular regions of the DB2 area.

*C*_{0}= 0.8 (Fig. 8, middle panel). In this case *C*_{crit 1} < *C*_{0} < *C*_{crit 2}.

The ‘23’ hierarchy region is now totally disconnected from the rest. The most dramatic change from Fig. 8 (left-hand panel) can be seen in the SB1 region. It is now predominantly regular. Chaotic regions can be seen only at the boundaries. Within the DB1 region at the boundary that has just disconnected from the SB1 region, a small chaotic zone exists. Thus it can be concluded that the chaotic behaviour is the result of the overlapping of the DB1 and SB1 regions.

The DB1 region is otherwise similar to the previous case.

Some important changes can be found in the DB2 region as well. The non-collision area has now grown to the full length of the DB2 region and is more regular as is indicated by the white spots in the RLI graph. It is very interesting that the big regular island at (0.2, 0.3) has completely disappeared. The ‘13’ (green) collision region has disappeared as well from the (0, 0.4) area. The chaotic non-collision island at (0.1, 0.2) has not only survived, but increased in size. In the SB2 region no difference can be seen.

It can be concluded that the phase space shows more regular behaviour. Big chaotic regions can be found only in the DB2 region. The most likely areas to find real physical systems are in the SB1 and DB1 regions.

*C*_{0}= 0.9 (Fig. 8, right-hand panel). In this case *C*_{0} is greater than *C*_{crit 2} and all hierarchy regions are disconnected from each other. All CSFBP orbits are hierarchically stable.

A dramatic change can be found in the DB2 region. The collision orbits have mostly disappeared and most of the orbits are regular.

The phase space is very regular everywhere. Real systems are likely to be found in each region.

Overall, Figs 7 and 8 show that there is a connection between the global hierarchical stability criteria and the chaotic behaviour of the phase space; namely, as *C*_{0} increases and the system becomes more hierarchically stable, the phase space becomes more regular.

It seems that strong chaos exists where the DB and SB regions connect and overlap. As the phase space becomes disconnected most of the orbits become regular, with chaotic regions existing only near the boundaries of real motion where collisions occur.

## Conclusions

We investigated the stable and chaotic behaviour of the CSFBP. Two fast chaos detection methods were applied to the system, the RLI and SALI methods.

We analysed the phase space in detail in the μ= 1, 0.1, 0.01 and 0.001 cases (Széll 2003). For greater μ values the CSFBP is a model for quadruple stellar systems, while for smaller μ values the CSFBP is a model for planetary systems. In this paper we presented the most representative results in the μ= 1 and μ= 0.1 cases. For the μ values we drew graphs of the phase space with different *C*_{0} values. The points of the graphs belong to different initial conditions of the CSFBP. In each graph we denoted the different kind of collisions by different colours and the chaotic and regular orbits with black and light grey, while the intermediate regions were grey-scaled. The main behaviours of the graphs were discussed in detail.

We found that the regular and chaotic behaviour of the phase space is closely connected with the Szebehely constant. The larger the Szebehely constant, the more regular the phase space becomes. For a Szebehely constant greater than the *C*_{crit 2} value where hierarchical stability is guaranteed, the phase space is regular almost everywhere. For smaller *C*_{0} values, chaotic behaviour occurs primarily near the boundaries of real motion, i.e. where collisions occur and in regions where the different hierarchies connect and overlap. Thus it can be seen that there exists a relationship between the global stability parameter (Szebehely constant) and the local regular and chaotic behaviour of the system. To the best knowledge of the authors no such relationship has been shown so far. In the future we will analyse the three-body problem from this point of view, as we believe that a similar connection between the global stability criteria (*c*^{2}*H* criteria) and the regularity of the phase space must exist.

We can assume that observable physical systems are mostly regular. If they were chaotic then sooner or later close encounters would happen between the bodies and the systems would fall apart. Thus the regular regions of the graphs can be connected to real systems. For large μ values, quadruple stellar systems might be found in the regular regions of the phase space, while for small μ values, planetary systems might be found in the ordered regions of the phase space.

The identification of the chaotic regions is also important. Finding any real system in the chaotic region of the phase space would indicate that the system was young.

### Acknowledgments

This work was supported by the British–Hungarian Scientific and Technological Foundation under grant number GB-53/01.