T Manos

E Athanassoula

Regular and chaotic orbits in barred galaxies -I. Applying the SALI/GALI method to explore their distribution in several models

Keywords: galaxies, kinematics and dynamics -galaxies, structure

The distinction between chaotic and regular behavior of orbits in galactic models is an important issue and can help our understanding of galactic dynamical evolution. In this paper, we deal with this issue by applying the techniques of the Smaller (and Generalized) ALingment Indices, SALI (and GALI), to extensive samples of orbits obtained by integrating numerically the equations of motion in a barred galaxy potential. We estimate first the fraction of chaotic and regular orbits for the two-degree-of-freedom (DOF) case (where the galaxy extends only in the (x, y)-space) and show that it is a non-monotonic function of the energy. For the three DOF extension of this model (in the z-direction), we give similar estimates, both by exploring different sets of initial conditions and by varying the model parameters, like the mass, size and pattern speed of the bar. We find that regular motion is more abundant at small radial distances from the center of the galaxy, where the relative non-axisymmetric forcing is relatively weak, and at small distances from the equatorial plane, where trapping around the stable periodic orbits is important. We also find that the variation of the bar pattern speed, within a realistic range of values, does not affect much the phase space's fraction of regular and chaotic motions. Using different sets of initial conditions, we show that chaotic motion is dominant in galaxy models whose bar component is more massive, while models with a fatter or thicker bar present generally more regular behavior. Finally, we find that the fraction of orbits that are chaotic correlates strongly with the bar strength.

INTRODUCTION

Exploring the nature of orbits in galaxies constitutes a very important issue, not only because of the evident astronomical interest in classifying the types of orbits that exist in such systems, but also because orbits are needed for constructing self-consistent models of galaxies. In order to study and understand the structure and dynamics of a galaxy it is necessary to identify the type of dynamics characterizing the motion of stars (regular or chaotic) and estimate the percentages of each type in the phase space of the galaxy.

From previous works that describe galaxies and their star motion, it is well-known that the analysis of periodic or-⋆ E-mail:thanosm@master.math.upatras.gr (TM); lia@oamp.fr (EA) bits, and their stability, can provide very useful information about galaxy structure. Stable periodic orbits are associated with regular motion, since they are surrounded by tori of quasi-periodic motion. Thus regular orbits are trapped in the vicinity of the parent periodic orbit. On the other hand unstable periodic orbits breed chaos. If we pick an orbit in the immediate phase-space vicinity of a chaotic one, we find that the two orbits will diverge exponentially with time. Such chaotic orbits fill up all the phase space region that is available to them. Nevertheless, recent results in galactic dynamics show that there are chaotic orbits that can support galaxy features, even thin features like the outer parts of bars, spirals and rings (Kaufmann & Contopoulos 1996;Patsis et al. 1997;Patsis 2006;Romero-Gómez et al. 2006[START_REF] Pettini | Geometry and topology in Hamiltonian dynamics and statistical mechanics[END_REF]Voglis et al. 2006aVoglis et al. ,b, 2007;;[START_REF] Contopoulos | [END_REF]; Athanassoula et al. 2009aAthanassoula et al. ,b, 2010;;Harsoula & Kalapotharakos 2009;Tsoutsis et al. 2009). These orbits are often called "sticky" and their true chaotic nature takes very long to be revealed [START_REF] Contopoulos | Order and chaos in dynamical astronomy[END_REF]. Athanassoula et al. (2010) termed these orbits "confined chaotic orbits" because the structures they generate are well confined in configuration space. There are also several recent related results showing that strong local instability does not mean diffusion in phase space (e.g. Giordano & Cincotta 2004;Cincota & Giordano 2008;Cachucho et al. 2010).

In our next paper, paper II, we will focus on this special category of chaotic orbits discussing their significance from an observational point of view.

Ferrers' potentials (Ferrers 1877) have proven very efficient for studying the main properties of real bars. The detection of periodic orbits and their stability in these models have been studied in detail by many researchers (e.g. Athanassoula et al. 1983;Pfenniger 1984;Skokos et al. 2002a,b;Patsis et al. 2002Patsis et al. , 2003a,b),b). It is well-known that stable periodic orbits of low period possess in their neighborhood sizeable domains of quasiperiodic (or regular) motion. A number of questions arises, therefore, concerning these regions of stability: How far from the stable periodic orbit, can regular motion be sustained? Where are chaotic regions located in phase space and configuration space? What are the model's parameters that favor large islands of stability around the main stable periodic orbits? In this paper we plan to address such questions and to provide some tentative answers.

Chaos arises in many areas of dynamical systems and its effect is generally related to the type of nonlinear terms present in the model's equations of motion (see [START_REF] Contopoulos | Order and chaos in dynamical astronomy[END_REF] for a review). The dynamics of stars in galactic potentials is a particular case, where nonlinear terms play an important role. The distinction between regular and chaotic motion is not trivial and becomes more subtle in systems of many DOF. Thus, it is necessary to use fast and precise methods to identify the nature of orbits in such models.

Many methods have been developed over the years dealing with this problem. The inspection of successive intersections of an orbit with a Poincaré surface of section (PSS) [START_REF] Lieberman | Regular and Stochastic Motion[END_REF] has been particularly useful for two DOF Hamiltonian systems. One of the most popular methods of chaos detection is the computation of the maximal "Lyapunov Characteristic Exponent" (LCE) (Oseledec 1968;Benettin et al. 1980a,b;Skokos 2010). Along the same lines as LCE, several other methods of chaos detection have been proposed based on the study of the behavior of deviation vectors: We may mention, e.g. the "Fast Lyapunov Indicator" (Froeschlé et al. 1997;Froeschlé & Lega 1998) and the "Mean Exponential Growth of Nearby Orbits" (Cincota & Simó 2000;Cincota et al. 2003). On the other hand, there also exist methods based on the analysis of time series constructed by the coordinates of each orbit such as the "Frequency Map Analysis" (Laskar 1990;Laskar et al. 1992;Laskar 1993). More details about these and other relative methods can be found in the book of [START_REF] Contopoulos | Order and chaos in dynamical astronomy[END_REF] and in the review of Skokos (2010) .

In the present paper, we use a method called the "Smaller ALingment Index" (SALI), based on the properties of two deviation vectors of an orbit for the quick and efficient distinction of chaotic motion, originally introduced by Skokos (2001). It has been successfully applied to different dynamical systems (Skokos 2001;Skokos et al. 2002;Voglis et al. 2002;Skokos et al. 2003a,b;Kalapotharakos et al. 2004;[START_REF] Skokos | [END_REF]Panagopoulos et al. 2004;Széll et al. 2004;Antonopoulos & Bountis 2006;Antonopoulos et al. 2006;Capuzzo-Dolcetta et al. 2006;Voglis et al. 2006aVoglis et al. , 2007;;Kalapotharakos et al. 2008;Manos et al. 2008a;Manos & Athanassoula 2008;Stránský et al. 2009;[START_REF] Macek | [END_REF]. In every case, it was confirmed to be a fast and reliable indicator of the chaotic or ordered nature of orbits. We also use SALI's recent generalization, the so-called "Generalized ALignment Index" (GALI) introduced by Skokos et al. (2007), using a set of p initially linearly independent deviation vectors of the system. Thus, following more deviation vectors, we manage to acquire more extra information about the complexity of the regular motion, i.e. the dimensionality of the invariant torus (or simply torus from now on) on which the orbit lies (Skokos et al. 2008;Manos et al. 2008b,c;Bountis et al. 2009;[START_REF] Manos | [END_REF].

The paper is organized as follows: Sections 2 and 3 are devoted to the chaos detectors we use, i.e. the Lyapunov spectra and the SALI/GALI methods. In Section 4, we present in detail the model, which is composed by a bulge, a disk and a Ferrers bar (Ferrers 1877). In Section 5, we present our results for the two DOF restriction of the general model, calculating the regular and chaotic orbits and constructing also charts of the different types of motion in the associated phase space. Section 6 is dedicated to the study of the full three DOF model in terms of computing percentages of regular and chaotic orbits, selecting different sets of initial conditions and varying the parameters of the bar component. Thus, with the help of the SALI method, we first distinguish the true nature of the studied orbits in phase space (Section 7) and correlate the bar force with the presence of larger amount of chaotic motion in phase space. Finally, in Section 8, we summarize our conclusions.

LYAPUNOV EXPONENTS

The Lyapunov Characteristic Exponents (LCEs) or Lyapunov Characteristic Numbers (LCN) or more simple Characteristic Exponents are very important for the study of dynamical systems, for distinguishing between regular and chaotic behavior of orbits in phase space (Benettin et al. 1980a,b;[START_REF] Lieberman | Regular and Stochastic Motion[END_REF][START_REF] Pettini | Geometry and topology in Hamiltonian dynamics and statistical mechanics[END_REF]Skokos 2010). In practice, the LCEs describe the rate of separation of infinitesimally close trajectories. The mathematical definition of LCE relies on Oseledečs multiplicative theorem (Oseledec 1968).

A flow x(t) generated by an autonomous first-order system is given by:

dx(t) dt = F (x(t)), ( 1
)
where F is its velocity field. Let us consider a trajectory in M -dimensional phase space together with a nearby trajectory, with initial conditions x0 and x0 + ∆x0, respectively. These evolve with time yielding the tangent vector ∆x(x0, t) with its Euclidean norm:

d(x0) = ∆x(x0, t) .
(2)

Writing ∆x = (∆x1, ..., ∆xM ) ≡ w, the time evolution for w if found by linearizing (1), to obtain the variational equations:

dw dt = J(x(t))w, (3)
where J(x(t)) = ∂F/∂x is the Jacobian matrix of the F (x).

The mean exponential rate of divergence of two initially close trajectories is:

σ(x0, w) = lim t→∞ ( 1 t ) ln d(x0, t) d(x0, 0) . ( 4
)
It can be shown that σ exists and is finite. Furthermore, there is an M -dimensional basis êi of w such that for any w, σ takes one of the M (possibly non-distinct) values σi(x0) = σ(x0, êi), ∀i = 1, 2, ..., M, which are the Lyapunov characteristic exponents. These can be ordered by size σ1 σ2... σM .

THE METHOD OF THE SALI (AND GALI) SPECTRA

Let us consider a Hamiltonian flow of N DOF, an orbit in the 2N -dimensional phase space with initial condition

P (0) = (x1(0), x2(0), ..., x2N (0)
) and two deviation vectors w1(0), w2(0) from the initial point P (0). In order to compute the SALI for that orbit one has to follow the time evolution of the orbit itself as well as the two deviation vectors w1(t), w2(t) which initially point in two arbitrary directions. The evolution of the deviation vectors is given by the variational equations (3) of the flow. At every time step the two deviation vectors w1(t) and w2(t) are normalized by setting:

ŵi(t) = wi(t) wi(t) , i = 1, 2 (5)
and the SALI is computed as (Skokos 2001):

SALI(t) = min { ŵ1(t) + ŵ2(t) , ŵ1(t) -ŵ2(t) } . (6)
The properties of the time evolution of the SALI rapidly distinguish between ordered and chaotic motion as follows:

The SALI fluctuates around a non-zero value for ordered orbits, while it tends exponentially to zero for chaotic orbits (Skokos et al. 2003a,b). In general, two different initial deviation vectors become tangent to different directions on the torus, producing different sequences of vectors, so that the SALI always fluctuates around positive values. On the other hand, for chaotic orbits, any two initially different deviation vectors in time tend to align in the direction defined by the maximal LCE (mLCE). Hence, they either coincide with each other, or become opposite, which leads to the SALI falling exponentially to zero. Thus, this completely different behavior of the SALI helps us distinguish between ordered and chaotic motion in Hamiltonian systems of any dimensionality. An analytical study of SALI's behavior for such orbits was carried out in [START_REF] Skokos | [END_REF], where it was shown that SALI ∝ e -(σ 1 -σ 2 )t , σ1, σ2 being the two largest Lyapunov exponents.

The usual technique to decide whether an orbit can be called chaotic or regular is to check, after some time interval, if its SALI has become less than a very small threshold value. In the following we will take this value to be equal to 10 -8 . Depending on the location of the orbit, this limit can be reached more or less fast, as there are phenomena that can hold off the final characterization of the orbit, and certain orbits behave as regular for long times before finally drifting away from regular regions and starting to wander in a chaotic domain.

The generalized alignment index of order p (GALIp) is determined through the evolution of 2 p 2N initially linearly independent deviation vectors wi(0), i = 1, 2, ..., p, so it is related to the computation of many LCEs rather than just the maximal one. The evolved deviation vectors wi(t) are normalized every few time steps in order to avoid overflow problems, but their directions are left intact. Then, according to Skokos et al. (2007), GALIp is defined to be the volume of the p-parallelogram having as edges the p unitary deviation vectors ŵi(t), i = 1, 2, ..., p:

GALIp(t) = ŵ1(t) ∧ ŵ2(t) ∧ ... ∧ ŵp(t) . ( 7
)
From the definition of GALIp it becomes evident that if at least two of the deviation vectors become linearly dependent, the wedge product in Eq. ( 7) becomes zero and the GALIp vanishes.

In the case of a chaotic orbit, all deviation vectors tend to become linearly dependent, aligning in the direction defined by the maximal Lyapunov exponent and GALIp tends exponentially to zero following the law:

GALIp(t) ∼ e -[(σ 1 -σ 2 )+(σ 1 -σ 3 )+...+(σ 1 -σp)]t , (8)
where σ1 > . . . > σp are approximations of the first p largest Lyapunov exponents. In the case of regular motion, on the other hand, all deviation vectors tend to fall on the Ndimensional tangent space of the torus, where the motion is quasiperiodic. Thus, if we start with p N general deviation vectors, these will remain linearly independent on the Ndimensional tangent space of the torus, since there is no particular reason for them to become aligned. As a consequence, GALIp in this case remains practically constant for p N . On the other hand, for p > N , GALIp tends to zero, since some deviation vectors will eventually become linearly dependent, following power laws that depend on the dimensionality of the torus. According to Christodoulidi & Bountis (2006) and Skokos et al. (2008) one obtains the following formula for the GALIp, associated with quasiperiodic orbits lying on k-dimensional tori (where k is potentially equal to, or smaller than the system's size):

GALIp(t) ∼    constant, if 2 p k 1 t p-k , if k < p 2N -k 1 t 2(p-N ) , if 2N -k < p 2N . (9)
In the case where k = N , GALIp remains constant for 2 p N and decreases to zero as ∼ 1/t zi.

(10)

The method has been applied successfully in several Hamiltonian systems like the FPU lattice (Skokos et al. 2008) and coupled symplectic maps (Bountis et al. 2009) for the detection not only of regular and chaotic motion but also the dimensionality of the torus on which a regular trajectory lies on.

Practically, the SALI is equivalent to GALI2 and the distinction between regular and chaotic motion in the 2D Ferrers barred model (two DOF) can be done with either one of them. For the full 3D version of the model, we also use GALI3, depending on the properties of the model (if its phase space is dominated by large chaotic or regular areas) and on our goals. GALI3 generally demands more CPU time, since it is necessary to follow the evolution of three deviation vectors instead of two. For the case of chaotic trajectories, however, GALI3 decays much faster and the calculation can stop well before the end of total time interval. Summarizing therefore, if someone had to estimate the amount of chaotic vs. regular regions in phase space, GALI3 would be more efficient for models where chaos is dominant, while SALI would be preferable for models with large regions of order. In our runs, we have used both SALI (GALI2) for the general description of chaotic vs. ordered regions, and GALI3 to follow specific regular orbits and study the dimensionality of the torus on which they lie.

THE MODEL POTENTIAL

The motion of a test particle in a 3D rotating model of a barred galaxy is governed by the Hamiltonian:

H = 1 2 (p 2 x + p 2 y + p 2 z ) + V (x, y, z) -Ω b (xpy -ypx). ( 11
)
The bar rotates around its z-axis (short axis), while the xdirection is along the major axis and the y along the intermediate axis of the bar. The px, py and pz are the canonically conjugate momenta, V is the potential, Ω b represents the pattern speed of the bar and H is the total energy of the orbit in the rotating frame of reference (Jacobi constant).

The corresponding equations of motion are:

ẋ = px + Ω b y, ẏ = py -Ω b x, ż = pz, (12) ṗx
= - ∂V ∂x + Ω b py, ṗy = - ∂V ∂y -Ω b px, ṗz = - ∂V ∂z .
The potential V of our model consists of three components:

(i) A disc, represented by a Miyamoto-Nagai disc (Miyamoto & Nagai 1975):

VD = - GMD x 2 + y 2 + (A + √ z 2 + B 2 ) 2 , ( 13
)
where MD is the total mass of the disc, A and B are its horizontal and vertical scale-lengths, and G is the gravitational constant.

(ii) A bulge, which is modeled by a Plummer sphere (Plummer 1911) whose potential is:

VS = - GMS x 2 + y 2 + z 2 + ǫ 2 s , ( 14
)
where ǫs is the scale-length of the bulge and MS is its total mass.

(iii) A triaxial Ferrers bar (Ferrers 1877), the density ρ(x) of which is:

ρ(x) = ρc(1 -m 2 ) 2 , m < 1 0 , m 1 , ( 15
)
where ρc = 105 32π GM B abc is the central density, MB is the total mass of the bar and

m 2 = x 2 a 2 + y 2 b 2 + z 2 c 2 , a > b > c > 0, ( 16
)
with a, b and c being the semi-axes. The corresponding potential is:

VB = -πGabc ρc n + 1 ∞ λ du ∆(u) (1 -m 2 (u)) n+1 , (17)
where

m 2 (u) = x 2 a 2 + u + y 2 b 2 + u + z 2 c 2 + u , ( 18
) ∆ 2 (u) = (a 2 + u)(b 2 + u)(c 2 + u), (19)
n is a positive integer (with n = 2 for our model) and λ is the unique positive solution of:

m 2 (λ) = 1, (20)
outside of the bar (m 1), and λ = 0 inside the bar. The corresponding forces are given analytically by Pfenniger (1984).

This model has been used extensively for orbital studies by Pfenniger (1984), Skokos et al. (2002a,b) and by Patsis et al. (2002Patsis et al. ( , 2003a,b),b). Our so-called standard (S) model has the following values of parameters G=1, Ω b =0.054 (54 km • sec -1 • kpc -1 ), a=6, b =1.5, c=0.6, A=3, B=1, ǫs=0.4, MB=0.1, MS=0.08, MD=0.82, both for its two DOF and three DOF versions. The units used are: 1 kpc (length), 1000 km • sec -1 (velocity), 1 M yr (time), 2 × 10 11 M (mass). The total mass G(MS + MD + MB) is set equal to 1.

Our study is mainly focused on understanding the effect of the bar on the dynamics of the model's phase space, by varying the bar mass MB, its semiaxis lengths b, c and its pattern speed Ω b . In order to do this we have considered three more models: models C, B and M whose short z-semiaxis (c-parameter), intermediate y-semiaxis (b-parameter) and mass of the bar (MB-parameter) are twice those of the initial reference model. For each of these four (S, C, B and M ) models we launch three different sets of initial conditions (distributions I, II, III). Hence, we consider the three standard models (IS, IIS, IIIS) and their variations (IC, IIC, IIIC), (IB, IIB, IIIB) and (IM, IIM, IIIM ), depending on the varying parameter. These sets of initial conditions will be discussed in detail in Section 7. The maximal time of integration of the orbits through the equations of motion and the variational equations is set to be T=10,000M yr (10 billion yrs), that corresponds to a time less than, but of the order of one Hubble time. We chose to take a maximum integration time that is the same for all orbits, independent of their energy. As a result, some of the orbits -and particularly the ones in the innermost regions of the galaxy -may have been calculated for a larger multiple of their own dynamical times than others. This choice, however, is necessary because it allows comparison with, or applications to a galaxy or a simulation, where all information refers to one time only, the same for all orbits. It was therefore also adopted in previous orbital structure studies in galactic models (e.g. El-Zant & Shlosman 2002, 2003;Valluri et al. 2010). This choice is essential in the case of some sticky orbits which may appear regular if calculated till e.g. 10,000M yr, but show clear signs of chaoticity when calculated for much longer times, i.e. orbits whose diffusion time scale is equal to many Hubble times. Such times, however, are irrelevant for our applications, and thus these orbits should not be counted as chaotic, even if in the long run they may become so. These issues will be discussed further in Paper II.

The relative strength of the non-axisymmetric forces can be estimated by the quantity Qt, defined as:

Qt(r) = (∂Φ(r, θ)/∂θ)max/(r∂Φ0/∂r), (21)
The maximum of Qt over all radii shorter than the bar extent is often referred to as Q b , and used as a measure of the bar strength (e.g. Buta et al. 2003Buta et al. , 2004;;Laurikainen et al. 2004;Buta et al. 2005;Durbala et al. 2009). From now on, we will for simplicity and conciseness, often refer to the "strong non-axisymmetric forcings" simply as "strong bars".

In

THE TWO DOF MODEL POTENTIAL

The two DOF Ferrers model is described by the Hamiltonian Eq. ( 11) setting (z, pz) = (0, 0), i.e.:

H = 1 2 (p 2 x + p 2 y ) + V (x, y) -Ω b (xpy -ypx). (22)
In the two bottom inset figures of Fig. 2a, we present two typical orbits in the (x, y)-plane of the standard model. One is regular (bottom left inset), with initial condition: (x, y, px, py) = (0, -0.625, -0.201, -0.06) (orbit R) and the other chaotic (bottom right inset), with initial condition: (x, y, px, py) = (0, -0.625, -0.002, -0.24) (orbit C), both for the Hamiltonian value H = -0.360. Their qualitatively different behavior is also shown in the PSS in Fig. 2a ((y, py)-plane). The chaotic orbit C (gray points) tends to fill with scattered points the chaotic region, while the ordered orbit R (black points) creates a set of points that form a closed invariant curve, on the left part of the picture. In the top inset figure of Fig. 2a, we show the different morphologies of three regular orbits around the periodic orbits in the centers of the three largest islands of regular motion present in the associated PSS. The model's basic barred shape, in the (x, y)-plane is mainly provided by trajectories around the main family of periodic orbits, the so-called x1 family (following the nomenclature of [START_REF] Contopoulos | [END_REF] which are orbits elongated along the bar, i.e along the x-axis in our case, and inside corotation. The stability island on the right gives rise to elliptical-like orbits, around the x2 family, elongated perpendicular to the bar, while the island on the left contains orbits which are elliptical-like but retrograde, very slightly elongated perpendicular to the bar and around the family x4.

We then apply the SALI method to the orbits R and C and present their different typical evolution behaviors: For the chaotic orbit C (gray curve in Fig. 2b), the SALI tends to zero (≃ 10 -16 ) exponentially after some transient time, while for the regular orbit R, it fluctuates around a positive number (black curve in Fig. 2b). The corresponding mLCE and σ1, for these orbits is shown in Fig. 2c, where the σ1 for the regular orbit tends to zero, while for the chaotic orbit it tends to a positive number after a long integration time.

This comparison shows clearly an advantage of SALI (GALI2), namely its ability to detect the chaotic character of an orbit faster than by simply following the mLCE, which typically takes a long time to converge. Even in our two DOF model, SALI starts to decrease exponentially after a relatively short time (in Fig. 2b,t 10 3 ) while the corresponding mLCE σ1 starts to converge to some positive non-zero value for t 10 3 (see gray curve in Fig. 2c).

Exploiting the efficiency of SALI, we take initial conditions on the (y, py)-plane (with x = 0) and calculate the values of the index to detect very small regions of stability (or instability) more globally. We are thus able to construct a map of chaotic and regular regions, very similar to what is depicted in a PSS, but with more accuracy and higher resolution (albeit at a somewhat longer CPU time). Furthermore, as a byproduct of our application, we obtain an accurate estimate of the percentage of regular to chaotic orbits on a surface of section of the given energy.

For example, choosing different values of the energy and using a sample of 50,000 initial conditions equally spaced on the same (y, py)-plane, we plot on the left column of Fig. 3 the PSS for H = -0.360, -0.335, -0.300, -0.260 and on the right the corresponding final SALI values obtained from the selected grid of initial conditions, where each point is colored according to its SALI values at the end of the integration. In the SALI plots, the light gray color corresponds to regular orbits, the black color represents the chaotic orbits/regions, while the intermediate gray shades between the two represent orbits with small rate of local exponential divergence, and/or orbits whose chaotic nature is revealed after long times, like for example the so-called sticky orbits, i.e. orbits that "stick" onto quasiperiodic tori for long times. Note that the fraction of these orbits (which lie mainly around the borders of the islands of stability) is very small (few % of the total amount of initial conditions) and can be discerned by eye in Fig. 3 only if one focuses on these particular regions. Thus, calculating from all these initial conditions the percentage of regular orbits, we are able to follow how this fraction varies as a function of the total energy. This problem was already addressed for four models by Athanassoula et al. (1983) but with fewer orbits, i.e. larger error bars. In the present model, the chosen energy values start from a regime of total order, at H = -0.46, up to values near the escape energy, Hesc = -0.20. The distinction between chaotic and regular orbit is that SALI< 10 -8 for chaotic orbits and 10 -8 for regular ones. In the latter range one does, of course, include "sticky" chaotic orbits as well. As is clear from Fig. 4, although the percentages of the regular orbits decreases sharply as the energy grows above H = -0.4, this tendency is reversed at higher energy values. Clearly, for H > -0.3, despite the fact that many stability islands have vanished, the size of the ordered domain increases because the main big island on the left of the PSS (see Fig. 3f,h) which corresponds to retrograde orbits, has become larger.

THE THREE DOF MODEL POTENTIAL

We now turn to the three DOF model and begin a study of ordered and chaotic domains. Before discussing the more general results about the study of the phase space of the model, we first present in detail the behavior of our chaos indicators for some typical orbits. Let us choose a regular and a chaotic orbit and compare the efficiency of the GALI method with that of the mLCE's.

In Fig. 5a,b we show the behavior of a chaotic orbit (C1), with initial condition: (x, y, z, px, py, pz) = (0.5875, 0.0, 0.33333, 0.0, 0.2, 0.0) and of a regular orbit (R1), with initial condition: (x, y, z, px, py, pz) = (0.97917, 0, 0.04167, 0, -0.17778, 0). For the chaotic one, the GALI2 (or SALI) has become almost zero for t ≃ 10 4 . However, as it is clear in Fig. 5a, the higher order GALI k , k = 5, 6, indicate the chaoticity of the orbit already by t ≃ 10 3 , since they have reached very small values at that time. Note the good agreement between the predicted slopes related to its two largest Lyapunov exponents (σ1 ≈ 0.00910 and σ2 ≈ 0.00345) given by Eq. ( 8). In Fig. 5b we show the GALI indices for the regular orbit R1, where both GALI2,3 ∝ constant and the GALI4,5,6 decay following the power laws:

GALI4(t) ∝ 1 t 2 , GALI5(t) ∝ 1 t 4 , GALI6(t) ∝ 1 t 6 , ( 23
) obtained from Eq. ( 9) for m = 3, indicating regular motion which lies on a 3D torus.

One example of a low dimensional motion (lower than 3) is provided by the regular orbit (R2), with initial condition: (x, y, z, px, py, pz) = (0.5875, 0.0, 0.29770, 0.0, 0.33750, 0.0). In Fig. 5c, we show the behavior of the Log(GALI) indices and their slopes for the regular orbit R2. Note that only GALI2 remains constant while the GALI3,4,5,6 tend to zero following power laws predicted by Eq. ( 9):

GALI2(t) ∝ const., GALI3(t) ∝ 1 t , GALI4(t) ∝ 1 t 2 , GALI5(t) ∝ 1 t 4 , GALI6(t) ∝ 1 t 6 , (24)
for m = 2. This implies that the orbit's motion lies on a 2D torus even though in three DOF Hamiltonian systems one generally expects the dimension of the torus to be three. In Fig. 6 we show the corresponding (x, y) and (x, z) projections for the above C1 (1st column), R1 and R2 orbits (2nd and 3rd column respectively). The chaotic one fills up its available regime. Regarding the two regular ones, there is a clear difference reflected in their projections. The "complexity" of R1 (quasiperiodic motion on a 3D torus) is more pronounced than the one of the orbit R2 (quasiperiodic motion on a 2D torus). Thus, in such cases GALI3 offers an extra advantage in helping us detect different "degrees" of regularity and does not serve only to distinguish between chaotic and regular motion.

THE DISTRIBUTION OF ORBITS IN PHASE SPACE

In this section we focus on the detection and quantification of the different kinds of orbital motion (regular and chaotic) in phase (and configuration) space, as some parameters of the bar component vary. A similar study, in different potentials, was done in El-Zant & Shlosman (2002, 2003), using the Lyapunov spectra as a tool. Here we concentrate mainly in the bar component and in understanding how its different properties can affect the general stability of the system. We use a relatively large number of trajectories (50,000), more than one set of initial conditions for the survey of the phase space, and the SALI method as a chaos detector for the reasons discussed in the previous sections. We discuss the relative fraction of regular and chaotic motion not only as a function of its spatial location but also as a function of the total energy and to correlate it with the bar strength, i.e. the relative non-axisymmetric forcing.

Initial conditions

We now need to define the sample of orbits whose properties (chaos or regularity) we will examine. The best of course would have been to draw these from a distribution function of the model. This, however, is not available, so we choose initial conditions on different grids in phase space or from the density distribution and the available energy range in the model and we measure the corresponding fraction of different kinds of motion. These percentages can depend on the choice of the sets of initial conditions and a priori are not expected to be equal. For this reason, before varying any other parameter we should consider as carefully as possible what initial conditions to choose and how well this choice "spans" the allowed phase system of the system. We first considered sets of orbits in such a way as to favor the dynamics near the main family of periodic orbits x1 (building blocks of the bar). To this end, we launched initial conditions along the bar's major x-axis with positive momenta py and position (z) or momentum values (pz) in the z-direction (distributions I and II below). We also considered initial conditions which follow the model's total density ρ in positions (x, y, z) while their positions are arbitrary (distribution III below). More specifically, the three different classes (distributions) of initial conditions we use here are:

• distribution I: 50,000 orbits equally spaced in the space (x, z, py) with x ∈ [0.0, 7.0], z ∈ [0.0, 1.5], py ∈ [0.0, 0.45] and (y, px, pz) = (0, 0, 0).

• distribution II: 50,000 orbits equally spaced in the space (x, py, pz) with x ∈ [0.0, 7.0], py ∈ [0.0, 0.35], pz ∈ [0.0, 0.35] and (y, z, px) = (0, 0, 0).

• distribution III: 50,000 orbits whose spatial coordinates are chosen randomly from the density distribution of our model, using the rejection method (Press et al. 1986) within the rectangular box -a

x a, -b y b, -c z c. We then draw a random value of the Hamilto- nian H (total energy) in the range [0, Hmax], where Hmax is 1.1 times the escape energy Hesc = -0.201 defined by the zero velocity curve going through the L1 or L2 Lagrangian point. Then we set (py, pz) = (0, 0) and we calculate the pxmomenta by the relation: px = H(x, y, z, py, pz) (keeping the positive solution). Note that in this case instead of giving py-momenta we tried px. This choice allows us to check how a different way of populating the phase space may affect the relative fraction of regular and chaotic motion. This third distributions has a better representation of the density, but still is arbitrary with respect to the velocities. Details of our distributions/models can be found in Table 1. In the next sections we start by studying the distribution of regular and chaotic motion in phase space, for different choices of initial conditions, different bar parameters (sizes, mass and pattern speed). We also briefly study the distribution as a function of energy and of location in configuration space.

Percentages of regular orbits vs. energy

We next examine the variation of the percentages of regular and chaotic orbits as a function of the energy and we show the results in Fig. 7a. We first divided the energy interval into 30 subintervals, each containing an equal number of orbits. In every subinterval we calculated the percentages of regular and chaotic orbits. In Fig. 7a we plot the percentage of regular orbits for energies up to H ≃ -0.15, i.e. up to a value somewhat beyond the escape energy for all models under study. Generally, we observe that the fraction of regular orbits decreases as the energy increases. An important peak, however, occurs before the escape energy, so that for (H > -0.20) the fraction of regular orbits increases again! This non-monotonic behavior is related to the growth of the islands around some basic stable periodic orbits in the phase space and is similar to what was observed in the two DOF model, comparing the fraction of regular orbits for the common energy interval, i.e. up to H -0.24 (see Fig. 4).

Picking another set of initial conditions we find similar qualitative results with some small quantitative differences. Hence, we find the same behavior when studying the models of distribution II, i.e. regular motion is dominant for small generally energies, then a decrease and then a peak. Distribution III manifests the same trend for small energies while the peak is less intense compared to the other two distributions.

Figure 8. Slices of the (x, py)-plane for different z values. Dark regions correspond to chaotic initial conditions, light gray to regular ones and intermediate colors to orbits with a small rate of local exponential divergence, and/or orbits whose chaotic nature is revealed after long times, like for example the so-called sticky orbits. The gray-scale bar represents the values of the SALI in logarithmic scale.

Spread of regular orbits in configuration space

We also explored the way in which regular and chaotic orbits are distributed along the z-direction of the configuration space. Following the evolution of each orbit, we calculated the mean absolute value of their z-coordinate (< |z| >) and plot in Fig. 7b the fraction of regular orbits as a function of the < |z| > for the models of distribution I. It clearly follows from these results that the intervals relatively near the (x, y)-plane (< |z| > < 0.35) contain mainly regular orbits, while those with larger values of mean absolute distance z are mostly chaotic. Similar qualitative results are obtained for distribution II (not shown here). Checking distribution III we find good match for the relative fraction of regular motion for relatively small < |z| > (up to 0.4) with the other models. For intermediate values (0.4 < < |z| > < 1.5) the motion is mostly chaotic while we don't have orbits reaching larger values, which is due to the way the initial conditions are given.

We also looked at these percentages as a function of the initial spherical radius (R spherical =

x 2 0 + y 2 0 + z 2 0 ) and the mean spherical radius over the evolution of the orbits (< R spherical >). Again, dividing the total range of the R spherical values in 30 subintervals, we calculated the percentages of regular orbits as a function of the mean R spherical value in each subinterval. We find that the fraction of regular orbits for all models decreases sharply with increasing R spherical up to R spherical < 1.5, where it reaches a minimum (see Fig. 7c). For 1.5 < R spherical < 7 it starts to increase gradually almost monotonically for model IS and with some rather week fluctuations for IC. Note that the model IB possesses a relatively large fraction of regular motion within the interval 2.5 < R spherical < 4.5 as well. Model IM follows a similar trend, but with a different scale than IB, containing its larger fraction of regular motion in the same roughly intervals, while for larger R spherical the motion is mainly chaotic. This result is in good agreement with the results in Fig. 7d, where the horizontal axis corresponds to the value of the mean spherical radius over the evolution time of the orbits. As previously, distribution II confirms well this behavior while distribution III is in good agreement for relatively small radial distances from the center of the galaxy and diversifies for larger ones.

Fraction of regular motion as a function of distance from the equatorial plane

We then examine whether the orbital motion is more "regular" near the (x, y)-plane or far from it and how this relates to their behavior as a function of < |z| > discussed in the subsection 7.3 and Fig. 7b. For this, we took the set of initial conditions from distribution I (version IC in this example) with initial conditions given in the space (x, z, py) and create a mesh in the (x, py)-plane with different slices in the z-direction that start from z = 0 up to z = 1.25 with step = 0.25. In Fig. 8, we see that as the distance of the initial conditions from the equatorial plane increases (from top to bottom row and from left to right panels), large islands of stability in the (x, py)-plane start to shrink, when x is between 0 and 2. This trend is almost monotonic up to the z-slice z ≈ 1 (Fig. 8, third row, left panel). For z = 1.25 (Fig. 8, third row, right panel) the stability region, lying between x = 1.5 and x = 2 starts to grow. This result turns out to be in good agreement with the one in Fig. 7b, where we plot the percentages of regular orbits as a function of < |z| >. This trend, as expected, is similar for the other two different sets of initial conditions (distributions IIC and IIIC), since the model parameters remain the same.

Percentages of regular orbits vs pattern speed

Furthermore, choosing one class of initial conditions (distribution IS), we compare different models for several values of the bar's pattern speed Ω b , keeping the remaining parameters constant. In this study, we try to explore the effect of the value of the pattern speed on the percentages of globally ordered or chaotic behavior of the system. Based on the structure and crowding of the periodic orbits, Contopoulos (1980) showed that bars have to end before corotation, i.e. that rL > a, where rL the Lagrangian, or corotation, radius. A similar result was reached from the calculation of the self-consistent response to a bar forcing (Athanassoula 1980). These two results are very useful, since they set a lower limit to the corotation radius, which can thus not be smaller than the bar length, but give no information on an upper limit. Comparing the shape of the observed dust lanes along the leading edges of bars to that of the shock loci in hydrodynamic simulations of gas flow in barred galaxy potentials, Athanassoula (1992a,b) was able to set both a lower and an upper limit to the corotation radius, namely rL = (1.2 ± 0.2)a. This restricts the range of possible values of the pattern speed for our model, from Ω b = 0.0367, that corresponds to the Lagrangian radius rL = 1.4a, to Ω b = 0.0554, that corresponds to rL = 1.0a. Using these extreme values, and the three intermediate frequencies: Ω b = 0.0403, Ω b = 0.0444 and Ω b = 0.0494, that correspond to the Lagrangian radii rL = 1.3a, rL = 1.2a and rL = 1.1a, respectively, we investigated how the value of the pattern speed affects the system and found that as Ω b increases the percentage of the regular orbits decreases relatively weakly. The corresponding results are given in Table 2. It turns out that the variation of the pattern speed does not affect drastically the relative fraction of regular and chaotic motion. One would probably need to try more extreme values of Ω b , beyond the upper and lower limits of the corotation radius, in order to see a significant change. Such values, however, would be unrealistic for real galaxies.

Fraction of regular and chaotic trajectories

In this section and in Fig. 9, we investigate the amount of regular motion in the phase space for all the three sets of initial conditions (distributions I, II and III), varying the axial ratios (b/a or c/a parameters) and the mass (MB parameter) of the bar component, as described in Table 1. In all distributions we used 50,000 initial conditions and we employed the bootstrap-method (Press et al. 1986) on several subsamples to make sure that this number is sufficiently high to provide reliable information. Computing the distributions of regular and chaotic orbits for the above classes of initial conditions we find the following:

7.6.1 Percentages for distribution I Before starting with the full 3D model, we first measured the percentages of regular orbits in the 2D model. We use initial conditions equally spaced in (x, py), with (y, px) = (0, 0),

x ∈ [0.0, 7.0] and py ∈ [0.0, 0.45]. We find relative fractions of regular motion equal to 45.40% for model IM , 53.43%, for model IS, 92.43% for model IB and 64.43% for model IC. This shows that at least in the 2D case, the massive bar model (IM ) has the highest fraction of chaotic orbits, followed by the standard distribution IS. The two models with the extended short (IC), or intermediate axis (IB) have the least chaos.

Let us now turn to the full 3D orbital coverage. In Fig. 9a we present the percentages of the regular orbits for the various sets of initial conditions of the distribution I.

Comparing the percentages, we see that increasing the mass of the bar in model IM increases chaotic behavior. This confirms the results obtained above and in Athanassoula et al. (1983) for the two DOF case. On the other hand, when the bar is thicker (model IC), i.e. the length of the z-axis larger, the system becomes more regular. Finally, the corresponding results for a fatter bar, i.e. with larger y-axis (model IB) demonstrate that the increase of the intermediate axis of the bar also provides the system with more ordered behavior. Comparing the 2D and 3D cases presented above, we note that the trends we find are the same, but that the relative fractions of regular orbits are somewhat higher in the 2D case. We will discuss this further towards the end of the section.

Percentages for distribution II

Similarly, in Fig. 9b, we present percentages of the regular orbits for initial conditions selected for distribution II and models S, M , B and C. As before, the increase of the bar mass (model IIM ) causes more extensive chaos, while for a thicker bar or fatter bar (models IIB and IIC), the system is again more regular. It is thus clear that the results of distributions I and II are similar, presumably dues to the fact that they contain initial conditions that cover the phase space in a similar manner, i.e. they favor the region around the x1 family of periodic orbits.

Percentages for distribution III

Finally, we plot in Fig. 9c the percentages of regular orbits for several sets of initial conditions from distribution III. A first general observation is that the fraction of ordered motion is significantly smaller than is distributions I and II, for all models. The basic reason for this difference is related to the way that the momenta of the positions are given in this distribution of initial conditions. In particular, all the orbits in this case are launched with a px momentum instead of a py as in distributions I and II. As a consequence, the big island of stability around the main family of periodic orbits x1 is populate with fewer orbits and thus in general more chaotic orbits are launched and measured. Nevertheless, again the increase of the bar mass in the IIIM model results in more chaotic behavior. As for the IB model, a thicker bar in the y-direction turns out to make the system more regular. Regarding the case of larger z-semiaxis (model IIIC), we notice a small decrease in the percentages of regular orbits (9.74% now, from 10.46%). It turns out that this particular different distribution of initial conditions doesn't reveal quite the same trend as in the previous cases (IS → IC and IIS → IIC). Concluding this section dedicated to the dynamical study of regular and chaotic motion in the phase space is that their corresponding fraction is basically related to the way that one populates the main stable (or unstable) periodic orbits of the system. We find that motion near the plane is generally stable (non chaotic) and that trajectories spending the largest fraction of time at large z are usually chaotic (with large Lyapunov exponents). These results are in agreement with the results obtained by El-Zant & Shlosman (2002, 2003) for different models. The increase of the bar mass leads to more chaotic motion in the phase space as found in Athanassoula et al. (1983) for a two DOF model, while in general the increase of the bar y or z-semiaxis leads to more regularity.

The effect of bar strength

Reviewing the above percentages of regular (and chaotic) motion in phase space, we see that they change in the same way when the basic parameters and properties of the model are varied. In order to understand this better we examine next the effect of the bar strength, measuring it as discussed in Section 4, i.e. from the relative strength of the non-axisymmetric forcing. In Fig. 10 we plot the relative fraction of chaotic motion as a function of bar strength Q b . Black filled squares correspond to distribution I, filled triangles to distribution II and filled circles to distribution III. This figure reveals that, provided the choice of initial conditions for the orbits is the same or similar, there is a tight correlation between the two quantities. the amount of chaos increases with increasing bar strength, as intuitively expected. What is unexpected, however, is how tight all these correlations are, with rank correlation coefficients ∼ -0.936 for all distributions.

Fig. 10 also shows that, except for the bar strength, the distribution of initial conditions chosen is also important. Distributions I and II differ little and their best fitting straight lines are roughly parallel and little displaced from each other. This is in good agreement with our previous statement (subsect. 7.1) that they both cover well the regular orbits trapped around the stable orbits of the x1 family, which are in fact the backbones of the bar (Athanassoula et al. 1983). On the other hand, distribution III has a much higher fraction of chaotic orbits, as could also be seen from Fig. 9c, attesting that the corresponding distribution of initial conditions is more chaotic, for the reasons we already described.

CONCLUSIONS

In this paper, paper I, we studied the detection and distribution of regular and chaotic motion in the phase space of barred galaxy models. Our results referring to the dynamical properties of the 2D/3D model with a Ferrers bar can be summarized as follows:

(i) We applied successfully the GALI method to distinguish both qualitatively and quantitatively between regular and chaotic orbits. We showed its efficiency and its advantage in detecting fast and accurately the chaotic nature of a trajectory.

(ii) In 2D systems, using the SALI (GALI2) method, we were able to identify efficiently and fast tiny regions of regular or chaotic motion, which are not clearly visible on PSSs.

(iii) Using GALI2,3 we detected regular motion in low dimensional tori, i.e. examples of regular orbits of the three DOF Ferrers' model that lie on a 2D torus, while the torus' expected dimension is generally three. Concerning chaotic orbits, the GALI3,4,5,6 indices decay faster than SALI and are able to detect chaos at early times well before this is evident from the calculation of the mLCE.

(iv) We also tried different models and different distributions of initial conditions for the orbits to test the fraction of regular and chaotic motion in various cases. We found that regular orbits are generally dominant at relatively small radial distances from the center of the galaxy and at small distances from the equatorial plane.

(v) We tested the effect of varying the bar pattern speed for our standard model and found that, within a realistic range of values, this parameter does not affect much the phase space dynamics. Within these limits, we found that the fraction of regular motion varies only by about 3 per cent, in the sense that the bars at the slow limit of the realistic range have more regular motion than the bars in the fast limit.

(vi) We varied the values of the main parameters of our models, such or the mass and the axial ratio of the bar component (in the 3D model). Chaos turns out to be dominant in galaxy models whose bar component is more massive, while models with a thicker or fatter bar present generally more regular behavior although the initial conditions are given can in general affect somewhat the relative percentages.

(vii) We found a very strong correlation (rank correlation coefficient at ∼0.94) between the fraction of chaotic orbits and the relative strength of the non-axisymmetric forcing. This holds for all our three initial conditions distributions taken individually, even though they were specifically chosen so as to represent different types of orbits. We can thus conclude that strong non-axisymmetric forcings (i.e. strong bars) are the main cause of the presence of large amount of chaotic motion in phase space.

In our next paper of this series (paper II) we focus on the significance of the "different degrees" of chaotic motion, like the strong and weak chaos, and their significance from an astronomical point of view. The fact that many "sticky"/chaotic orbits often persist and behave in a regular manner for very large time intervals, before showing their chaoticity, makes them astronomically important for supporting the galaxy structure. We present there a first approach in filtering these weak chaotic orbits from the strongly chaotic ones and classify them observationally as effectively ordered motion.

  Fig. 1 we show the above quantity for the 4 models we are going to study and discuss in this paper. From the quantity Qt it becomes clear that the increase of the mass of the bar component (model M : dotted line) gives rise to a "stronger bar" compared to the initial standard model (model S: solid line). On the other hand, an increase of the short z-semiaxis (c-parameter -model C: dashed line) or intermediate ysemiaxis (b-parameter -model B: dot-dashed line) has as a result "weaker bars". Hence, how do "strong bars" affect the relative regular and chaotic motion in the phase space of the full three DOF model? Before answering this question, let us first study and comprehend the dynamics of the two DOF model.

Figure 1 .

 1 Figure 1. The quantity Qt as a function of radius. Its maximum value corresponds to Q b , which gives an estimate of the relative strength of the non-axisymmetric forces. We see that the increase of the mass of the bar component (model M : dotted line) corresponds to "stronger bar" compared to the initial standard model (model S: solid line). On the other hand, when we increase the short z-semiaxis (c-parameter -model C: dashed line), or the intermediate y-semiaxis (b-parameter model B: dash-dotted line) the bar becomes weaker.

Figure 2 .

 2 Figure2. Analysis of orbits in the two DOF case: a) Poincaré surface of section in the (y, py)-plane for a regular (R-black color) and a chaotic (C-gray color) orbit (H = -0.36). The points of the successive intersections of the regular orbit R with this plane create a closed curve, while the points of the chaotic orbit C fill with scattered points all its available region of motion. The orbits can be seen in the two inset figures in the bottom of the panel. In the top inset figure, we show examples of the three different morphologies of the regular orbits around the three main periodic orbits in the centers of the three main islands of stability. Note that the basic barred shape in the (x, y)-plane is mainly provided by trajectories from the central stable region, around the family of periodic orbits x1 (orbits elongated along the bar's major axis). The other two stability islands give either elliptical-like orbits elongated perpendicular to the bar orbits (right island -family x2), or orbits which are elliptical-like but retrograde with respect the bar's pattern speed (left island -family x4). b) The corresponding evolution of the SALI for the R and C orbits of panel a. For the chaotic one C (gray line), the SALI decays exponentially fast to zero, while for the ordered one R (black line) it fluctuates around a non-zero number. c) The corresponding evolution of the maximal Lyapunov exponents σ 1 for the orbits R and C. For the chaotic orbit C (gray line) this tends after some transient time to a non-zero value, while for the regular orbit R (black color) it tends linearly to zero. Note that the Log(SALI) in y-axis of panel b ranges from -16 to 2, while the Log(σ 1 ) in panel c ranges from -5 to 0.

Figure 4 .

 4 Figure 4. Percentages of regular orbits for several values of the energy (H = -0.460,-0.435,-0.410,-0.360,-0.335,-0.300,-0.260,-0.240) in the two DOF Ferrers model. Note that while the fraction of the regular orbits decreases sharply as the energy grows above H = -0.4, this tendency is reversed at higher energy values.

Figure 3 .

 3 Figure 3. The PSS (left column) and SALI (right column) methods are in good agreement when surveying the same projection of the phase space. In the left column we show the PSS for the two DOF model potential with H = -0.360 (panel a), H = -0.335 (panel c), H = -0.300 (panel e) and H = -0.260 (panel g). In the right column we present regions of different values of the SALI for 50,000 initial conditions on the (y, py)-plane for the same values of the Hamiltonian (panels b,d,f and h, respectively). The light gray colored areas correspond to regular orbits, while the dark black ones to chaotic. Note the excellent agreement between the two methods as far as the gross features are concerned, as well as the fact that the SALI can easily pick out small regions of stability which the PSS has difficulties detecting, like those around the central islands in panels b,d,f and h.

Figure 5 .

 5 Figure5. Distinguishing chaotic from regular motion with GALI method. a) Exponential decay of all GALIs for the chaotic orbit C1 (lin-log scale) with the predicted slopes related to its two largest Lyapunov exponents (σ 1 ≈ 0.00910 and σ 2 ≈ 0.00345). b) Slopes of GALIs for the regular orbit R1 showing that its motion lies on a 3D torus, where GALI 2,3 ∝ constant while the GALI 4,5,6 decay following a power law. c) Slopes of GALIs for the regular orbit R2 lying on a 2D torus with only GALI 2 ∝ constant.

Figure 6 .

 6 Figure6. Two projections ((x, y) and (x, z)-planes) of the chaotic orbit C1 (1st column) and the quasiperiodic orbits R1 (2nd column) and R2 (3rd column). Note the different "complexity" between R1, whose regular motion lies on a 3D torus, and R2, which lies on a 2D as GALI indices detected.

Figure 7 .

 7 Figure 7. Percentages of regular orbits as a function of (a) the energy, (b) the averaged < |z| > value, (c) the initial spherical radius R spherical , and (d) the averaged < R spherical > radius of the orbit's evolution, for models IS, IB, IC and IM .

Figure 9 .

 9 Figure 9. Percentages of regular orbits (SALI 10 -8 ) for distribution I (left panel), II (middle panel) and III (right panel) of initial conditions for all the different model versions, where the bar mass and the length of its intermediate and minor (vertical) semi-axes are varied. Generally, when the bar mass increases the model tends to be more chaotic, while when the size of the bar semi-axes (b or c parameters) increases the fraction of regular motion grows. The small difference in the trends between models IS → IC and IIIS → IIIC is due to the different way in which the sample of initial conditions in phase space is chosen.

Figure 10 .

 10 Figure 10. The relative fraction of chaotic motion as a function of bar strength Q b increases (see Fig.1). Filled squares ( ) correspond to the models M , S, C and B (from right to left, respectively) for distribution I. while filled triangles ( ) are for distribution II and filled circles (•) for distribution III.

Table 1 .

 1 Varying parameters for all sets of initial conditions.

	distribution/model	M B	b	c
	IS, IIS, IIIS	0.1	1.5	0.6
	IC, IIC, IIIC	0.1	1.5	1.2
	IB, IIB, IIIB	0.1	3.0	0.6
	IM, IIM, IIIM	0.2	1.5	0.6

Table 2 .

 2 Models with different pattern speeds Ω b .

	model	Ω b	r L	% of Regular orbits
	IΩ 1	0.0367032	1.4a	28.39
	IΩ 2	0.0403014	1.3a	27.63
	IΩ 2	0.0444365	1.2a	27.26
	IΩ 4	0.0493654	1.1a	26.55
	IΩ 5	0.0554349	1.0a	25.55

The reason we choose an extended range of energies is because we wish to study some orbits that are in this range but don't necessarily escape during our adopted integration time.

ACKNOWLEDGEMENTS

We would like to thank Ch. Skokos, T. Bountis, A. Bosma, M. Romero-Gómez and P.A. Patsis for their fruitful comments and discussions on this work. We also thank the referee, P.M. Cincotta, for his report, which helped us improve the clarity of the paper. T. Manos acknowledges the "Karatheodory" graduate student fellowship No B395 of the University of Patras and funds from the Marie Curie fellowship No HPMT-CT-2001-00338 to Marseille observatory. Visits between the Marseille and Patras team members were partially funded by the region PACA in France. This work was partially supported by the grant ANR-06-BLAN-0172.


Claire Michel

Valérie Doya
email: valerie.doya@unice.fr

Sorin Tascu

Wilfried Blanc

Olivier Legrand

Fabrice Mortessagne

CONTROLLED EXCITATION OF SCAR MODES IN PASSIVE AND ACTIVE MULTIMODE CHAOTIC FIBERS

A multimode optical fiber with a D-shaped cross section is an experimental paradigm of a wave system with chaotic rays' dynamics. We show that seldom but usable modes, called scar modes, localized along some particular direction of the geometrical trajectories can be selectively excited. We report numerical simulations that demonstrate the importance of the so-called self focal point in the scar modes selection process. We use a localized illumination in a passive fiber, or, a localized gain in an Ytterbium-doped fiber, located in the vicinity of this special point to control the scar modes selection.

Introduction

The behavior of waves in a 2-dimensional bounded domain (billiard) is determined by the nature of the rays' dynamics in the geometrical limit. For simple shapes as the circle or the square, the rays' dynamics is regular but tends to be complex as the shape becomes more complicated. When the rays' dynamics presents the property of high sensitivity to initial conditions, the dynamics is chaotic. The influence of the chaotic nature of the dynamics on the properties of waves constitutes the active field of investigations of Wave Chaos [START_REF] Stöckmann | Quantum Chaos: an introduction[END_REF]. From vibrating plates [START_REF] Schaadt | Experiments on elastomechanical wave functions in chaotic plates and their statistical features[END_REF], to microwave cavities [START_REF] Stöckmann | Quantum Chaos" in billiard studied by microwave absorption[END_REF][START_REF] Sridhar | Experimental observation of scarred eigenfunctions of chaotic microwave cavities[END_REF] or optical fiber [START_REF] Doya | Speckle statistics in a chaotic multimode fiber[END_REF], all these various systems present generic features well-described by theoretical predictions that demonstrate the universality of Wave Chaos. Investigating how the chaotic geometrical limit influences the wave properties is not only a fundamental issue but reveals also applications in various domains ranging from room acoustics to reverberant electromagnetic cavities. Optics in particular is getting advantages of wave chaos. Chaotic fibers [START_REF] Doya | Optimized absorption in a chaotic doubleclad fiber amplifier[END_REF] are used to optimize the pump power absorption in double clad fiber amplifier and chaotic shapes provide improved features to microlasers [START_REF] Harayama | Lasing on scar modes in fully chaotic microcavities[END_REF][START_REF] Fang | Wave interference effect on polymer microstadium laser[END_REF][START_REF] Fang | Analysis of high-quality modes in open chaotic microcavities[END_REF].

As mentionned above, a multimode optical fiber with a truncated (D-shaped) cross section (see fig. 1) and large dimensions compared to the wavelength has been proved to be an efficient experimental tool to investigate the spatial properties of waves that propagate in a chaotic system [START_REF] Doya | Speckle statistics in a chaotic multimode fiber[END_REF]. Indeed, as the optical fiber is invariant along the direction of propagation z, the projection of the z-evolution of a ray over its transverse cross section can be seen as the time-evolution of a ray in the D-shaped billiard (see fig. 1(b)). The chaotic nature of the rays' dynamics in the D-shaped billiard controls the properties of the propagating modes in the wave-analogous system. Generic modes of chaotic systems are ergodic modes: their intensity is statistically uniformly distributed all over the cross section (fig. 1(c)). This behavior is the result of multiple interferences encountered by the waves that reflect on the chaotic boundary [START_REF] Berry | Regular and irregular semiclassical wavefunctions[END_REF]. Nevertheless, some surprising modes may appear: they present strong localization of intensity along the direction of some periodic trajectories of the associated billiard. For these non-ergodic modes, called scars [START_REF] Heller | Bound-State Eigenfunctions of Classically Chaotic Hamiltonian Systems: Scars of Periodic Orbits[END_REF], the spreading of intensity is "frozen" by constructive interferences along a periodic orbit. The intensity of the field of different scar modes along the short 2-bounce periodic orbit (p.o.) is presented in fig. 2. These singular modes have been the object of numerous theoretical investigations [START_REF] Berry | Quantum scars of classical closed orbits in phase space[END_REF][START_REF] Kaplan | Scars in quantum chaotic wavefunctions[END_REF][START_REF] Li | Numerical study of scars in a chaotic billiard[END_REF] and, recently, have been used in the field of applied optics. Their localized enhancement of intensity offers interesting properties: scar modes are currently employed in microlaser cavities to increase the directivity of laser emission and decrease the laser threshold [START_REF] Harayama | Lasing on scar modes in fully chaotic microcavities[END_REF][START_REF] Fang | Wave interference effect on polymer microstadium laser[END_REF][START_REF] Fang | Analysis of high-quality modes in open chaotic microcavities[END_REF]. These scar modes reveal some signatures of the underlying geometrical limit. They mark some periodic trajectory and concentrate their intensity over a specific point. In this paper, we utilize this "geometrical point" to excite selectively the scar modes which are difficult to exhibit in a multimode fiber. We show that in a passive multimode fiber, a localized excitation can select some scar modes. In order to be less dependent on the illumination condition, we propose, in a second part, to perform a "scar mode amplifier" by introducing a localized gain to control wave chaos. The conclusive part of the present paper is dedicated to the description of experimental feasibility of the scar modes amplifier.

Scar selection by localized excitation

The system we use is a silica fibre with n core =1.451 surrounded by a low refractive index silicone cladding with n cladding =1.41. The fibre is truncated at half its radius R=60µm. Scar modes can exist along a periodic orbit if µL << 2π where µ is the Lyapunov coefficient (i.e. the rate of divergence of two close initial conditions measured in the phase space) and L the length of the orbit. When this condition is fulfilled, constructive interferences along the periodic orbit freeze the development of the chaotic dynamics leading to localized scar modes. We consider the 2-bounce p.o. which is the shortest periodic orbit of the D-shaped billiard, with L = 3R and µ = 0.4/R, so that scar modes along the 2-bounce p.o. can exist. The condition of constructive interferences along the 2-bounce p.o. defines the transverse wavenumber κ p associated to a scar mode of order p (fig. 2(a)-(f)):

κ p L -∆ϕ - π 2 = 2πp (1)
∆ϕ is the phaseshift induced by the reflection at the core/cladding interface that depends on κ p , n core and n cladding and p is an integer. The additional π/2 phaseshift is due to the unique self-focal point [START_REF] Bogomolny | Smoothed wave functions of chaotic quantum systems[END_REF] on the 2-bounce p.o. A pencil of rays emitted from a self-focal point will focus on this same point after one period. Moreover, we note that the D-shaped cross section mimicks an unstable Fabry-Perot resonator: the self-focal point is close to the focal point of the curved dielectric mirror. This special point concentrates rays and also corresponds to a maximum of intensity for the scars that can be exploited for selective excitation. As seen in fig. 2, most of the scar modes along the 2-bounce p.o. exhibit a maximum of intensity in the vicinity of the self-focal point. This distinctive feature gives rise to the technique that we propose for achieving a selective excitation of scar modes: use an initial illumination with a Gaussian shape centered on this intensity spot. At the input of the fiber, the field of the Gaussian beam is decomposed over a finite number of guided modes whose weights are given by their spatial overlap with the incoming light. In order to preferentially excite the scar modes, the width of the Gaussian is set to the mean diameter of the bright spot which characterizes these modes, in practice σ = 5.64 µm. To describe the spectral contents of the field evolving along the fiber, we introduce the κz-spectrum C(κ, z) defined as the Fourier transform of the correlation function of the field ψ( r) with its initial condition [START_REF] Feit | Computation of mode properties in optical fiber waveguides by a propagating beam method[END_REF][START_REF] Michel | Selective Amplification of Scars in a Chaotic Optical Fiber[END_REF]:

C(κ, z) = z+(Lz/2) z-(Lz/2) dz ′ d rψ * ( r, 0)ψ( r, z ′ )e -iβ(κ)z ′ ( 2
)
where β is the longitudinal wavenumber and L z a propagation length. We use an algorithm implementing the Beam Propagation Method to simulate the propagation of the initial field along the fiber [START_REF] Doya | Speckle statistics in a chaotic multimode fiber[END_REF][START_REF] Feit | Light propagation in graded-index optical fibers[END_REF]. The figure 3(a) shows the κ-z-spectrum calculated after a propagation length that ensures that all the modes are well established. This figure has to be compared to the calculated spatial overlap between the initial Gaussian beam and the scar modes shown in fig. 3(b). We note that only a few modes dominate the spectrum: the scar modes from order 1 to 5. The excited modes correspond to the one having a maximal spatial overlap with the Gaussian input field. Thus, as the Gaussian beam is launched on the high intensity area of the scar modes, its spatial overlap is larger for scar modes than for the generic speckle modes. In a "passive" multimode chaotic fiber, we use the complex spatial structure of the modes to select some of them. One stigma of the geometrical dynamics, the self-focal point, arises in the spatial distribution of intensity of the scar modes and can be exploited to control the selection of scar modes. A perfect localization of excitation on the self-focal point can act as a scar modes selector but this method relies on a perfect command on the location and characteristic of the initial beam.

Scar selection using gain

The challenge is to overcome the intrinsic limitation of a passive system by introducing an active medium inside the fiber. The main question is: how to determine the optimized distribution of the gain ? A hint to answer may be found in a closely related domain. Indeed, some recent investigations in disordered media with gain have shown that a given localized mode can be selectively excited through a pumping confined to its spatial extension [START_REF] Sebbah | Random laser in the localized regime[END_REF].

In our case, due to technical constraints, there is no hope that the active medium precisely match the spatial pattern of a given scar mode. But here again the high intensity area close to the self-focal point will be helpful: we consider an Ytterbium-doped multimode fiber with the gain region localized on the position of the bright spot. To prevent light from being guided in the doped area, we consider a small but realistic refractive index difference between the silica core and the doped area of 5.10 -4 . To pump the system, a laser at 980 nm is focused in the large chaotic truncated core. The large number of ergodic modes excited ensures an uniform absorption of the pump along the propagation [START_REF] Doya | Optimized absorption in a chaotic doubleclad fiber amplifier[END_REF][START_REF] Leproux | Experimental study of pump power absorption along rare-earth-doped double clad optical fibres[END_REF]. The field of the pump is absorbed by the active ions and acts as a power reservoir for the signal. In order to evaluate the efficiency of the modes selection by the localized gain, a spatially incoherent signal (speckle field) at 1064 nm is coupled into the D-shaped core. The initial coupled energy is therefore spread over a great number of modes.

In figure 4(a), we report the far-field intensity of the signal at the input of the fiber integrated over the radial coordinate φ κ , I(φ κ ). The inset presents the initial far-field intensity I(r κ , φ κ ). All the wavenumber directions are explored. We let the field evolve along the fiber: amplification occurs and after 40 m of propagation along the exotic fiber amplifier, the maximum amplification factor is achieved. We note z max this length of maximum amplification. We plot the same quantities as (fig. 4(a)) in (fig. 4(b)) for z = z max . We can observe that a strong transverse wavenumber filtering process occurs during the propagation. The integrated far-field intensity tends to concentrate over two privileged angles that correspond to values of φ κ equal to π/2 and 3π/2. These angles are associated to the direction of the trans-verse wavenumber of the waves that bounces back along the direction of the 2-bounce p.o.

In the far-field pattern, we can also notice that some peaks concentrate the intensity along the 2-bounce p.o. exhibiting a selection on the transverse wavenumber not only on direction but also on the modulus. With the aim to get more informations on the mode selection process induced by the gain, we analyze the transverse wavenumber spectrum calculated for the length z max . We notice that few selected modes appear in the wavenumber spectrum: only some modes get amplified along the propagation. In table 6, the theoretical values of the scar mode wavenumbers deduced from Eq. ( 1) are compared to the wavenumbers corresponding to the peaks in the numerical spectrum. The comparison clearly establishes that each amplified wavenumber is associated to a scar mode wavenumber. This result and the particular structure of the far-field pattern shown in fig. 4(b) conspire to demonstrate that the amplification mechanism applies only to the scar modes. From the spatial overlap of a given scar mode with the doped region, one can even obtain its amplification rate thus accounting for the relative amplitudes of the different amplified scar modes. The diamonds in fig. 5 show the values taken by these overlaps. They are compared to the amplification rates directly deduced from the evolutions of the scar mode powers. The typical evolution of the power along the fiber is exemplified in fig. 6 for the scar mode of order p = 4. In a semi-logarithmic scale, the power exhibits a linear growth before reaching a plateau which corresponds to the saturation of the amplification mechanism. The amplification rate is then obtained through a linear fitting procedure by assuming a uniform amplification process of the form:

P p (z) = P p (0)exp(α p z), (3)
The constant slope of the curve gives the amplification rate α p . The crosses in fig. 5 represent the values of α p for the different scar mode orders. The good agreement between the two quantities emphasizes the importance of the overlap of the field with the doped area on the differential modes amplification process. The location of the gain controls the scar mode amplification: to get an efficient amplification, the doped area should have a good interaction with the scar modes.

First step to experimental realization

The experimental investigation of the selective amplification process put numerically into evidence requires the realization of a non-conventional fiber amplifier. This has been achieved and the experiments is under progress. In a first step, we performed a silica preform doped with 1250 ppm of Ytterbium ions by means of a standard MCVD (Modified Chemical Vapor Deposition) method. Ytterbium ions are incorporated without other ions as Aluminium usually used to ensure a good incorporation of rare-earth ions in the silica matrix. The index profile measured on the preform confirms the presence of Ytterbium ions by the way of a refractive index increase of 5.10 -4 . Due to this small refractive index difference between the doped area and the core, the fundamental mode of the weakly guiding doped area is spatially extended and suffers from important losses. As a consequence, this mode is not supported along the propagation. The MCVD method imposes the active ions to be located in the center of the preform. In order to localize the doped area in the vicinity of the self-focal point and also to get a D-shaped silica core, the circular preform is machined: a smaller preform with an off-centered doped area is cut into the preform [START_REF]Optical Development Workshop in the Institute of Plasma Physics[END_REF]. Then, the preform is pulled by controlling the temperature and the drawing velocity to get the fiber with a well-preserved D-shaped geometry and diameter over hundred of meters. We use a low refractive index polymer (Luvantix UVF PV-409R04AP) as the optical cladding. The figure 7 is the absorption spectrum of the Ytterbium-doped fiber measured experimentally.

The general behavior corresponds to the generic absorption spectrum of Ytterbium in silica matrix. From this measure, we can deduce an estimation of the absorption coefficient at the pump wavelength, a p = 6.2dB/m. To perform the scar mode amplifier, we realize the following experimental setup. A pump laser diode at 978 nm is focused at the fiber input. The signal from a Yag laser at 1064 nm is collimated and launched at the fiber input. A dichroic mirror is used to couple both the signal and the pump at the fiber input. For the signal illumination, we can either used a focused beam or a diffuse beam (for a noisy field illumination). We depict in the figure 8 the evolution of the gain curve along the propagation. The gain is obtained from the measured signal at the fiber output in presence of the pump over the measured signal without pump. The curve increases up to a critical length z max ≃ 15 m associated to the saturation length of the signal. We note that the value of the gain is small compared to usual optical amplifier. Indeed, this amplifier is unusual from different aspects: as instance, the quantity of incorporated Ytterbium is small (and without aluminium co-doping) and the signal propagates in a large D-shaped core reducing the interaction with the doped area is reduced. Nevertheless, the gain curve clearly shows that an amplification process occurs along the propagation: optical amplification is achieved with this exotic active fiber. The next step is to control the initial illumination so that the amplification occurs on the scar modes. This step is under process and will constitute the purpose of a forthcoming publication.

Conclusion

In systems that exhibiting deterministic chaotic in the geometrical limit, some modes show surprising enhancement of intensity along geometrical periodic orbit of the corresponding billiard. These modes are of particular interest because they freeze the chaotic nature of the system and tend to concentrate intensity in finite regions of the system. Nevertheless, they are isolated and sparse modes. We demonstrate that the self-focal point proper to the 2-bounce periodic orbit can be used to select the scar modes in a billard-analogous wave system: a D-shaped optical fiber. Intensity of the scar modes is concentrated on the selffocal point: we prove that using a localized excitation in a passive fiber, or introducing a localized gain on this specific point, we can excite selectively the scar modes. Thus, a pure geometrical quantity, the self-focal point, may play a crucial role on the wave properties. These numerical investigations have been done using realistic datas from measurements and the experiment of scar amplifier is currently in progress.

Table caption

Tab1 The table reports the theoretical value of the transverse wavenumber for each scar mode of order p along the 2-bounce p.o. compared to the measured value of the transverse wavenumber associated to the peak (amplified modes) on the spectrum.tab1.ps

Fig1Fig2Fig3Fig5Fig6Fig7

  Fig1 Transverse cross section of the D-shaped multimode fiber (a), typical rays' trajectory in the associated D-shaped billiard (b) and generic speckle mode of a D-shaped chaotic optical fiber (c).fig1.eps

FIGURE 8

 8 FIGURE 8

TABLE 1

 1

	p κ p from eq.1 measured κ p
	1	4.69	4.70
	2	6.78	7.06
	3	8.86	9.05
	4	10.95	10.99
	5	13.03	13.10
	6	15.12	15.29
	7	17.20	17.23
	8	19.29	19.59
	9	21.37	21.53
	10	23.46	23.68
	11	25.02	25.80
	12	27.10	27.93
	13	29.19	29.80

Acknowledgement

We want to address grateful thanks to S. Trezien and M. Udé for their technical supports in the preform and fiber manufacture. We also acknowledge fruitful discussion with P. Aschieri. This work is part of the projet n o 05-JCJC-0099-01 supported by the Agence Nationale de la Recherche.


Sho Fujimura

Shuji Shiraishi

A remark on the generation of all Euler trails for an Eulerian graph with 2 n edges

Keywords:

For an Eulerian graph G with 2 n edges., we give a new method of generating all Euler trails in G in parallel efficiently by using 2-chain cover of G.

A trail is a walk with no repeated edges. An Eulerian trail in a graph G is a walk in G passing through every edge exactly once. A graph G has an Eulerian trail if and only if it has either 0 or 2 vertices of odd degree. A graph is eulerian if it admits an Euler trail.

Since we deal with only graphs G with 2 n edges, we can consider an Euler trail a sequence of 2-chains of G. A 2-chain is a simple path of length two. A 2-chain cover of G is a set of edge-disjoint 2-chains that covers all edges of G.

A 2-chain cover of G is obtained efficiently in parallel by the following manner [1]: let T denote extended spanning tree of G and let f denote the number of descendants function of T . Let v be a vertex in T

(1) If f (v) is even, we can get the edge-disjoint 2-chains by using all edges in {(v, v i )|f (v i )} is even.

(2) If f (v) is odd and z is the parent of vertex v, we can get the edge-disjoint 2-chains by using all edges in {(z, v)|(v, v i ), f (v i )} is even.

We use the following fact mentioned in [START_REF] Fujimura | On the number of perfect matchings of line graphs[END_REF]:

Fact. If G has an even number of edges, all 2-chain covers of G are obtained.

By the above Fact an Euler trail corresponds to an either one of 2-chain covers in G. Now consider the converse. Let C be a 2-chain cover of G. We regard each 2-chain c = {n, m} in C (n, m are edge labels) as one new edge of which endpoints are the same of the 2-chain c. The new edge is then labeled by nm.

We call here such reduction 2-chain contraction of G and denote the resulting graph by G/2. We repeat above 2-chain contraction process until the edge size of G/2 i is two or G/2 j is disconnected.

If the edge size of G/2 i is two, we obtain an Euler trail by writing down the labels of edges of G/2 i . If G/2 j is disconnected, the 2-chain cover does not correspond to an Euler trail.

We summarize this result as follows:

Let G be an eulerian graph with 2 n edges. Then all Euler trails in G are generated in parallel efficiently.

About the complexity of 2-chain cover of G, see [START_REF] Shiraishi | A parallel algorithm for the maximum 2-chain edge packing problem[END_REF]. We give the above procedure by considering a detailed example. Let G = (V, E) , where V = {a, b, c, d, e, f }, E = {1, 2, 3, 4, 5, 6, 7, 8} (Fig. 1a) and its spanning tree S with root a. (red line in Fig. 1b). Note that we can choose an arbitrary spanning tree. The construction of an extended spanning tree based on S using non-tree edges labelling 1,4, 6 of G\S is shown in Fig. 2b T Below are the lists of all 2-chain covers of G:

{1,3}, {2,5}, {4,8}, {6,7} … C 1 {1,3}, {2,7}, {4,8}, {5,6} … C 2 {1,3}, {2,4}, {5,8}, {6,7} … C 3 {1,3}, {2,7}, {4,6}, {5,8} … C 4 {1,2}, {3,4}, {5,8}, {6,7} … C 5 {1,6}, {2,7}, {3,4}, {5,8} … C 6 {1,2}, {3,8}, {4,5}, {6,7} … C 7-1 {1,4}, {2,5}, {3,8}, {6,7} … C 7-2 {1,5}, {2,4}, {3,8}, {6,7} … C 7-3 {1,4}, {2,7}, {3,8}, {5,6} … C 8-1 {1,5}, {2,7}, {3,8}, {4,6} … C 8-2 {1,6}, {2,7}, {3,8}, {4,5} … C 8-3
Case C 1 : The procedure of generating Euler trails by 2-chain cover C 1 is the following: First, we construct the graph Next for the 2-chain contracted graph G ′ and its spanning treeS ′ , we obtain the 2-chain covers

G ′ = G/2 by 2-chain contraction of G based on C 1 . G ′ = (V, E), where V = {b, c, d, e}, E = {13,
C ′ 1 = {{31, 48}, {25, 67}}, C ′ 2 = {{31, 67}, {25, 84}} of G ′ (
Like the Case C 1 , 2-chain covers C ′ 1 = {{31, 48}, {27, 65}}, C ′ 2 = {{27, 31}, {65, 84}} of G ′ and 2-chain covers C ′′ 1 = {3148, 5672}, C ′′ 2 = {3127, 6584} of G ′′ are obtained.
C 1 '' C 2 '' Fig. 8: 2-chain covers C ′ 1 , C ′ 2 of G ′ and 2-chain covers C ′′ 1 , C ′′ 2 of G ′′

S'

Fig. 9 shows 2-chain covers C 4 '' Fig. 12 shows 2-chain covers Fig. 14 shows the 2-chain cover

C ′ 1 = {{31, 58}, {24, 67}}, C ′ 2 = {{31, 67}, {42, 85}} ,C ′ 3 = {{31, 85}, {24, 67}}, C ′ 4 = {{31, 42}, {58, 67}} of G ′ and 2-chain covers C ′′ 1 = {3148, 5672}, C ′′ 2 = {3127, 6584}, C ′′ 3 = {3148, 5672}, C ′′ 4 = {3127, 6584} of G ′′ .
Fig. 9: 2-chain covers C ′ 1 , C ′ 2 , C ′ 3 , C ′ 4 of G ′ and 2-chain covers C ′′ 1 , C ′′ 2 ,C ′′ 3 , C ′′ 4 of G ′′
C ′ 1 = {{31, 58}, {27, 64}}, C ′ 2 = {{27, 31}, {46, 85}} ,C ′ 3 = {{27, 64}, {31, 85}}, C ′ 4 = {{27, 58}, {31, 46}} of G ′ and 2-chain covers C ′′ 1 = {3158, 4672}, C ′′ 2 = {5846, 7213} ,C ′′ 3 = {1385, 2764}, C ′′ 4 = {1346, 7258} of G ′′ . Note that C ′′ 1 , C ′′ 3 have loops.
C 4 '' Fig. 10: 2-chain covers C ′ 1 , C ′ 2 , C ′ 3 , C ′ 4 of G ′ and 2-chain covers C ′′ 1 , C ′′ 2 ,C ′′ 3 , C ′′ 4 of G ′′
C ′ 1 = {{43, 85}, {21, 67}}, C ′ 2 = {{67, 85}, {21, 34}} of G ′ and 2-chain covers C ′′ 1 = {8543, 1276}, C ′′ 2 = {8567, 2134} of G ′′ .
C 2 '' Fig. 11: 2-chain covers C ′ 1 , C ′ 2 of G ′ and 2-chain covers C ′′ 1 , C ′′ 2 of G ′′
C ′ 1 = {{43, 85}, {27, 61}}, C ′ 2 = {{27, 85}, {34, 61}} of G ′ and 2-chain covers C ′′ 1 = {8543, 1672}, C ′′ 2 = {8527, 6134} of G ′′ .
C ′ 1 = {{45, 83}, {12, 76}} of G ′ and 2-chain cover C ′′ 1 = {4583, 1276} of G ′′ .
′ 1 = {{38, 41}, {52, 76}} of G ′ and 2-chain cover C ′′ 1 = {4138, 5276} of G ′′ .

Fig. 1 :

 1 Fig. 1: G and its spanning tree S with root a

Fig. 2 :Fig. 3 :

 23 Fig. 2: Vertex split and extended spanning tree T1 based on S

Fig. 4 :

 4 Fig. 4: All 2-chain covers of G

  Fig. 5: 2-chain contracted graph G ′ based on C1

Fig. 6 :Fig. 7 :

 67 Fig.6:Spanning tree S ′ of G ′ and 2-chain covers C ′ 1 , C ′ 2 of G ′

3 :

 3 The figure below shows that G ′ is the 2-chain contracted graph of G based on C 3 and its spanning tree S ′ .

  Fig.10 shows 2-chain covers C′ 1 = {{31, 58}, {27, 64}}, C ′ 2 = {{27, 31}, {46, 85}} ,C ′ 3 = {{27, 64}, {31, 85}}, C ′ 4 = {{27, 58}, {31, 46}} of G ′ and 2-chain covers C ′′ 1 = {3158, 4672}, C ′′ 2 = {5846, 7213} ,C ′′ 3 = {1385, 2764}, C ′′ 4 = {1346, 7258} of G ′′ . Note that C ′′ 1 , C ′′ 3 have loops.

Case C 5 :

 5 Fig.11shows 2-chain coversC ′ 1 = {{43, 85}, {21, 67}}, C ′ 2 = {{67, 85}, {21, 34}} of G ′ and 2-chain covers C ′′ 1 = {8543, 1276}, C ′′ 2 = {8567, 2134} of G ′′ .

6 :

 6 The figure below shows that G ′ is the 2-chain contracted graph of G based on C 6 and its spanning tree S ′ .

Fig. 12 :Fig. 13 :

 1213 Fig. 12: 2-chain covers C ′ 1 , C ′ 2 of G ′ and 2-chain covers C ′′ 1 , C ′′ 2 of G ′′

  Fig. 14: 2-chain cover C ′ 1 of G ′ and 2-chain cover C ′′ 1 of G ′′

Fig. 16 :

 16 Fig. 15: 2-chain cover C ′ 1 of G ′ and 2-chain cover C ′′ 1 of G ′′