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 c2H which, for a given critical value c2Hcrit, 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, C0. 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 P1 and P4 on opposite sides of the centre of mass of the system O having mass m and the other two bodies P2 and P3 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 r4=−r1; r3=−r2; V4=−V1; V3=−V2. Note that V1 and V4 do not need to be coplanar with V2 and V3. 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.

Figure 1.

The initial configuration of the coplanar CSFBP.

Figure 1.

The initial configuration of the coplanar CSFBP.

Figure 2.

The parallelogram nature of the coplanar CSFBP after time t.

Figure 2.

The parallelogram nature of the coplanar CSFBP 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:  

(1)
formula
m1=m is the mass of P1 and m2=M is the mass of the more massive body P2. The force function U can be written as:  
(2)
formula
where ri is the length of the radius vector of Pi from the centre of mass for i= 1, 2 and r12 is the distance between P1 and P2.

The energy integral becomes  

(3)
formula

Let E0=−E. The moment of inertia I may be written as  

(4)
formula

Then Sundman's inequality,  

(5)
formula
which relates the moment of inertia I, the kinetic energy T and the magnitude of the angular momentum vector c, becomes:  
(6)
formula

For a given energy and angular momentum, equation (6) gives boundary surfaces in the r1r2r12 space which define the regions of real motion.

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

(7)
formula
Let us now introduce dimensionless variables ρ1, ρ2, and ρ12, and a new dimensionless constant C0 called the Szebehely constant:  
(8)
formula
In terms of these variables, Sundman's inequality (6) becomes  
(9)
formula
and the kinematic constraint becomes  
(10)
formula

Four types of collisions are possible in this problem.

  • ‘12’ collision (double binary (DB) collision). If ρ12= 0 then P1 collides with P2 and P3 collides with P4. The inequalities (10) are satisfied only if ρ21.

  • ‘13’ collision (DB collision). If graphic, then P1 collides with P3 and P2 collides with P4.

  • ‘23’ collision (single binary (SB) collision). If ρ2= 0 then P2 collides with P3. The inequalities are satisfied only if ρ121.

  • ‘14’ collision (SB collision). If ρ1= 0 then P1 collides with P4. The inequalities (10) are satisfied only if ρ122.

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 C0= 0. Each tube represents a particular hierarchical arrangement of the bodies as follows (see Fig. 3).

Figure 3.

The four possible hierarchical arrangements in the CSFBP.

Figure 3.

The four possible hierarchical arrangements in the CSFBP.

‘12’-type hierarchy (Fig. 3a). The P1P2 bodies and the P3P4 bodies form binaries which orbit around the baricenter of the four-body system. r12 is small, while r1 and r2 are of similar size. This is a double binary (DB) system.

‘13’-type hierarchy (Fig. 3b). The P1P3 bodies and P2P4 bodies form binaries which orbit around the baricenter of the four-body system. r12 is large, while r1 and r2 are of similar size. This is a DB system.

‘23’-type hierarchy (Fig. 3c). The P2P3 bodies form a central binary and the P1 and P4 bodies orbit around it. r2 is small, while r1 is large. This is a single binary (SB) system.

‘14’-type hierarchy (Fig. 3d). The P1P4 bodies form a central binary and the P2 and P3 bodies orbit around it. r1 is small, while r2 is large. This is a SB system.

For C0= 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 C0= 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.

Figure 4.

The projection of the phase space on to the ρ1–ρ2 plane for the case of a mass ratio of μ= 0.1 for increasing values of C0. The dotted regions are forbidden for real motion.

Figure 4.

The projection of the phase space on to the ρ1–ρ2 plane for the case of a mass ratio of μ= 0.1 for increasing values of C0. The dotted regions are forbidden for real motion.

When the Szebehely constant is allowed to increase from C0= 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 C0, the more extended these forbidden regions become, until for a critical value of C0=Ccrit 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 C0 reaches a second critical value of C0=Ccrit 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 C0 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, Ccrit 1=Ccrit 2= 46.841 and for μ= 0.1, Ccrit 1= 0.7792 and Ccrit 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 r1, r2 and r12 to r1, r2 and α, the angle formed by P1OP2. The corresponding conjugate momenta are p1, p2 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):  

(11)
formula
where c is the constant angular momentum.

The equations of motion can then be written as:  

(12)
formula
 
(13)
formula
 
(14)
formula
 
(15)
formula
 
(16)
formula
 
(17)
formula

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 p1=p2= 0 so that V1 and V2 are perpendicular to the line the bodies lie on.

We will vary the initial r1 and r2 variables with fixed energy (E0) and Szebehely constant (C0) values.

The initial value of pα and the angular momentum constant c are determined from E0 and C0, through equations (8) and (11). We recall that H=−E0.

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  

(18)
formula
where ξ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:  
(19)
formula
where Df(x) is the Jacobian matrix evaluated at x. If L1(x*) = 0, the orbit evolving from x* is regular. If L1(x*) > 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 x0:  
(20)
formula
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):  
(21)
formula
where δx*≪ 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 Δ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.

Figure 5.

The time evolution in time steps Δt of the RLI numbers for a– regular, b– weakly chaotic and c– chaotic orbits.

Figure 5.

The time evolution in time steps Δt of the RLI numbers for a– regular, b– weakly chaotic and c– chaotic orbits.

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*:  

(22)
formula
If x* is in a regular region, 〈ΔL(x*)〉t* is small. Otherwise, if x* is in a chaotic region, 〈Δ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:  

(23)
formula
while the antiparallel alignment index is:  
(24)
formula
The time-evolution of the deviation vectors can be calculated according to equation (19). From the definition of the alignment indices it follows that 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.

Figure 6.

The time evolution in time steps Δt of the SALI numbers for a– regular, b– weakly chaotic and c– chaotic orbits.

Figure 6.

The time evolution in time steps Δt of the SALI numbers for a– regular, b– weakly chaotic and c– chaotic orbits.

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 C0 which covered the complete range of different possibilities for global hierarchical stability. i.e. C0 < Ccrit 1 (unstable); Ccrit 1 < C0 < Ccrit 2 (‘23’ hierarchy stable) and Ccrit 2 < C0 (stable). The values of μ, E0 and C0 chosen are given in Table 1.

Table 1.

The μ, C0, and E0 values for which we investigated the phase space of the CSFBP.

Table 1.

The μ, C0, and E0 values for which we investigated the phase space of the CSFBP.

For any given μ, E0 and C0 each CSFBP orbit is uniquely defined by its initial values r1 and r2. The nature of each orbit can therefore be depicted in r1r2 space by colour code. Using this fact we can construct r1r2 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 r1 and r2 covering all regions of real motion in the r1r2 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.

Table 2.

Criteria and colour codes for the different categories of orbits.

Table 2.

Criteria and colour codes for the different categories of orbits.

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 r1r2 plane, where each point on the plane represents the starting conditions of an orbit. Each point in the r1r2 plane is colour coded according to the nature of the orbits as defined in Table 2.

For each μ–C0-E0 triple we produce two r1r2 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 r1r2 space forμ= 1.0

In this case E0 was fixed to be 7 and C0 was chosen to be 10, 40 and 60. Recall that Ccrit 1=Ccrit 2= 46.841. The initial values of r1 and r2 were varied from 0 to 4 with a step size of 0.01. We used units such that GM2= 1.

In the case where all bodies have equal mass, an interchange of the r1 and r2 produces the same physical orbit. Thus the categorization of the orbits is symmetric with respect to the r1=r2 line with the exception that the ‘23’ and ‘14’ collisions are reversed. That is, for r1r2 the initial ordering of the equal masses is 1234. For r2r1 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 r1r2 space to r2r1 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 r1r2 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 r1=r2= 1: see Fig. 7. The top line of the three graphs shows the RLI categorization of the orbits for C0= 10, 40 and 60, while the bottom line shows the SALI categorization.

Figure 7.

μ= 1, E0= 7. Graphs of the r1, r2 initial conditions for C0= 10, 40, 60 left to right. The top line of three graphs show the RLI categorization of the different (r1, r2) orbits, while the bottom line of three graphs show the SALI categorization for the same C0 values. The colours indicate collisions: red –‘12’ type, green –‘13’ type, yellow –‘14’ type, blue –‘23’ type. The light grey regions are regular, the dark grey regions are undetermined, the black regions are chaotic, and the white regions are forbidden for real motion.

Figure 7.

μ= 1, E0= 7. Graphs of the r1, r2 initial conditions for C0= 10, 40, 60 left to right. The top line of three graphs show the RLI categorization of the different (r1, r2) orbits, while the bottom line of three graphs show the SALI categorization for the same C0 values. The colours indicate collisions: red –‘12’ type, green –‘13’ type, yellow –‘14’ type, blue –‘23’ type. The light grey regions are regular, the dark grey regions are undetermined, the black regions are chaotic, and the white regions are forbidden for real motion.

The graphs can be separated into two regions: the double binary (DB) region around r2r1 and the two single binary regions around r2r1 (SB1) and r1r2 (SB2). In the DB region there always exists two massive regular regions, one on either side of the r2=r1 line, independent of the C0 value. For small C0 values, the SB regions contain primarily collision orbits which indicates chaotic behaviour. As C0 increases the number of the collision orbits decreases. For all values of C0, 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.

C0= 10 (Fig. 7, left-hand panel). In this case C0 is well below the critical value of Ccrit 1=Ccrit 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.

C0= 40 (Fig. 7, middle panel). In this case C0 is just below the critical value of Ccrit 1=Ccrit 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 r1 values in the SB1 region, the islands are not well determined; they are dark grey in the RLI graph. For large r1 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 r1 values and similarly in the SB2 region for greater r2 values.

C0= 60 (Fig. 7, right-hand panel). In this case C0 is much greater than the critical value Ccrit 1=Ccrit 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 r1r2 space forμ= 0.1

In this case E0 was fixed to be 1.2 and C0 was chosen to be 0.4, 0.8 and 0.9. Recall that for μ= 0.1 Ccrit 1= 0.7792 and Ccrit 2= 0.8886. r1 and r2 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 r1 and r2 values. Thus we need to integrate the whole range of the initial parameters. In this case we integrated graphic 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 r1= 2.2 and r2= 1.5 in order to be able to see more of the detail: see Fig. 8. As in Fig. 7 for each C0 value, we have two graphs showing the categorization of the orbits according to the RLI and SALI numbers.

Figure 8.

μ= 0.1, E0= 1.2. Graphs of the r1, r2 initial conditions for C0= 0.4, 0.8, 0.9 from left to right. The top line of graphs show the RLI categorizations and the bottom the SALI. The colour coding is as before in Fig. 7.

Figure 8.

μ= 0.1, E0= 1.2. Graphs of the r1, r2 initial conditions for C0= 0.4, 0.8, 0.9 from left to right. The top line of graphs show the RLI categorizations and the bottom the SALI. The colour coding is as before in Fig. 7.

The graphs can be separated into four regions: two double binary regions around r2r1 but r2 < r1 (DB1) and r2 > r1 (DB2), and two single binary regions, r2r1 (SB1) and r1r2 (SB2). In the DB1 region there always exists a massive regular region just below the r2=r1 line, independent of the C0 value. For small C0 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 C0 increases, the number of collision and chaotic orbits decreases. In the following we analyse each graph in detail.

C0= 0.4 (Fig. 8, left-hand panel). In this case C0 is well below the critical value Ccrit 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 r2≈ 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. r2= 0 and r2= 0.5), there exists islands of undetermined behaviour, dark grey in the RLI graph and light grey in the SALI graph. For large initial r1 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.

C0= 0.8 (Fig. 8, middle panel). In this case Ccrit 1 < C0 < Ccrit 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.

C0= 0.9 (Fig. 8, right-hand panel). In this case C0 is greater than Ccrit 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 C0 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 C0 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 Ccrit 2 value where hierarchical stability is guaranteed, the phase space is regular almost everywhere. For smaller C0 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 (c2H 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.

References

Froeschlé
C.
Lega
E.
Gonczi
E.
,
1997
,
Celest. Mech. & Dyn. Astr.
 ,
67
,
41
Marchal
C.
Saari
D. G.
,
1975
,
Celest. Mech.
 ,
12
,
115
Pál
A.
Sándor
Zs.
,
2002
, Publ. Astron. Dept. Eötvös Univ.,
12
,
27
Roy
A. E.
Ovenden
M. W.
,
1955
,
MNRAS
 ,
115
,
296
Roy
A. E.
Steves
B. A.
,
2001
,
Celest. Mech. & Dyn. Astr.
 ,
78
,
299
Sándor
Zs.
Érdi
B.
Efthymiopoulos
C.
,
2000
,
Celest. Mech. & Dyn. Astr.
 ,
78
,
113
Skokos
Ch.
,
2001
,
J. Phys. A (Math. Gen)
 ,
34
,
10029
Steves
B. A.
Roy
A. E.
,
2001
, in
The Restless Universe: Application of Gravitational N-Body Dynamics to Planetary, Stellar and Galactic Systems
 .
IOP Publishing
,
Bristol
, p.
301
Széll
A.
,
2003
,
PhD thesis
 . Glasgow Caledonian Univ.
Voglis
N.
Contopoulos
G.
Efthymiopoulos
C.
,
1998
,
Phys. Rev. E
 ,
57
,
372
Zare
K.
,
1976
,
Celest. Mech.
 ,
14
,
73
Zare
K.
,
1977
,
Celest Mech.
 ,
16
,
35