The phasespace of boxypeanut and Xshaped bulges in galaxies I. Properties of nonperiodic orbits
Abstract
The investigation of the phasespace properties of structures encountered in a dynamical system is essential for understanding their formation and enhancement. In the present paper we explore the phase space in energy intervals where we have orbits that act as building blocks for boxypeanut (b/p) and “Xshaped” structures in rotating potentials of galactic type. We underline the significance of the rotational tori around the 3D families x1v1 and x1v1 that have been bifurcated from the planar x1 family. These tori play a multiple role: (i) They belong to quasiperiodic orbits that reinforce the local density. (ii) They act as obstacles for the diffusion of chaotic orbits and (iii) they attract a large number of chaotic orbits that become sticky to them. There are also bifurcations of unstable families (x1v2, x1v2). Their unstable asymptotic curves wind around the x1v1 and x1v1 tori generating orbits with hybrid morphologies between that of x1v1 and x1v2. In addition, a new family of multiplicity 2, called x1mul2, is found to be important for the peanut construction. This family produces stickiness phenomena in the critical area of the radial and vertical inner Lindblad resonances (ILRs) and reinforces b/p bulges. Our work shows also that there are peanutsupporting orbits before the vertical ILR. Nonperiodic orbits associated with the x1 family secure this contribution as well as the support of b/p structures at several other energy intervals. Nonlinear phenomena associated with complex instability of single and double multiplicity families of periodic orbits show that these structures are not interrupted in regions where such orbits prevail. Depending on the main mechanism behind their formation, boxy bulges exhibit different morphological features. Finally our analysis indicates that “X” features shaped by orbits in the neighbourhood of x1v1 and x1v1 periodic orbits are pronounced only in sideon or nearly endon views of the bar.
keywords:
Galaxies: kinematics and dynamics – chaos – diffusion structure1 Introduction
The skeleton of the boxy and peanutlike (b/p) structures observed in edgeon disc galaxies is considered to be the families of periodic orbits (hereafter p.o.), which constitute the x1tree. These are the families bifurcating from the planar x1 family (Contopoulos, 1981) which lies on the equatorial plane, at the vertical resonances (Skokos et al., 2002a, b; Patsis et al., 2002). Despite the fact that the periodic orbits determine their environment in the phasespace, the main body of the structures are nonperiodic orbits. It is known that even in 2D Hamiltonian systems quasiperiodic orbits trapped around stable periodic orbits do not necessarily support a structure of similar morphology. A morphology similar to the one of the p.o. is secured mainly by quasiperiodic orbits on invariant curves close to it. Contrarily, quasiperiodic orbits on invariant curves at the outer parts of stability islands may support more complicated structures. For example quasiperiodic orbits trapped around elliptical x1 orbits in 2D Ferrers bars can support an ansae type bar, if orbits closer to the borders of the x1 stability islands are populated (Patsis, 2005). This independence between the morphologies of the p.o. and of the structures around them is expected to be more pronounced in 3D systems. Also the extent and the significance of the phenomenon of stickiness (Contopoulos & Harsoula, 2010; Bountis et al., 2012) in building a 3D bar and a b/p bulge has not been investigated up to now.
In a recent series of papers (Katsanikas & Patsis, 2011; Katsanikas et al., 2011a, b, 2013) we have presented the dynamical fingerprints of nonperiodic orbits close to p.o. in 3D rotating Hamiltonian potentials of galactic type. For this purpose we used the colourrotation method proposed by Patsis & Zachilas (1994) in order to visualise the 4D spaces of sections. In this way we have associated typical structures in the 4D space with the neighbourhood of stable, as well a with simple, double and complex unstable periodic orbits (Broucke, 1969; Hadjidemetriou, 1975; Contopoulos & Magnenat, 1985). Studies of other Hamiltonian systems of galactic type (Contopoulos & Harsoula, 2013), or even studies of maps (Lange et al., 2014; Richter et al., 2013; Zachilas et al., 2013), in which the 4D space is visualised, indicate that the main structures depend basically on the corresponding kind of instability, or stability, of the p.o. In the present paper we take advantage of the colourrotation technique in order to investigate the phase space of the b/p and Xshaped structures in our model. This method has been proven to be especially efficient in estimating the role of sticky orbits (Katsanikas et al., 2013).
Referring to b/pshaped bulges we mean boxy structures, the morphology of which can be schematically described either with the symbol “” or with the symbol “”. Bureau et al. (2006) and Aronica (2006) characterise these two qualitatively different profiles of edgeon disc galaxies as “OX” and “CX” respectively. In many cases wings of an “X” feature delineate the peanut shape of the bulge, so the terms “OX” and “CX” refer to an offcentred or centred “X” respectively.
We assume that the X feature traces the extrema of the b/p bulge. Let us assume that it reaches a maximum radius along the major axis and a maximum height vertical to the equatorial plane . As we concluded in Patsis et al. (2002) the b/p bulge is part of the 3D bar, which is in agreement with the results of studies with different approaches to the same problem (see e.g. Lütticke et al., 2000; Bureau & Athanassoula, 2005; Bureau et al., 2006). Observationally, the estimation of the extent of the bar in edgeon galaxies is usually based on identifying “bumps” in the surface brightness profiles in cuts along or parallel to their major axes as indicated originally by Wakamatsu & Hamabe (1984). The qualitative and quantitative analysis of such profiles allowed Lütticke et al. (2000) to estimate the ratio of the length of the bar to the length of the b/p distortion to be about . Aronica (2006) (his appendix C) presented the unsharpmasked profiles of a sample of 25 edgeon galaxies with b/p bulges, which reveal a clear X feature (see also Bureau et al., 2006). By measuring approximately the dimensions of the X features in the unsharpmasked images of this sample, one can empirically estimate roughly a variation of the ratio, to be . We estimated a mode value close to . This quantity depends primarily on the orientation of the bar with respect to the line of sight and is expected to be small for bars viewed nearly endon and large for bars viewed nearly sideon, since will not vary substantially. This is important information for deciding which orbits populate a b/p structure. The goal of our study is to identify orbits which could reinforce these structures and reproduce the observed ratios of their geometrical dimensions.
The study of the body of the orbital content of b/p bulges can provide valuable information about the properties of the stellar component of standard galactic bars similar to the bar of the Milky Way (Saito et al, 2011; Li et al., 2012; Wegg & Gerhard, 2013; Gardner et al., 2014). Nonperiodic orbits reflect the dynamical mechanism that determines the volume and the regions of the phase space, which they occupy. Moreover they provide kinematic imprints that can be compared with observed quantities (Kregel et al., 2004; Williams et al., 2011; Vásquez et al., 2013) or with results of body simulations (Athanassoula & Misiriotis, 2002; Bureau & Athanassoula, 2005; Saha & Gerhard, 2013; Gardner et al., 2014; Quillen et al., 2014). This comparison can lead in favouring or excluding dynamical mechanisms and help us understand the internal structure of the peanuts. The first step towards this comparison is to register all types of orbital behaviour, regular and chaotic, that leads in profiles observed in edgeon disc galaxies.
The paper is structured as follows: The model used and a brief review on the stability of p.o. in 3D Hamiltonian systems is given in section 2. In section 3 we discuss the origin of the b/p bulge, investigating the lowest energy orbits that could support a b/p feature. We present the dynamical phenomena in the neighbourhood of the main families of p.o. in the region of the radial and vertical ILRs in section 4. In section 5 we investigate the role of chaotic orbits in the sideon profiles based on the dynamics of the x1v1 family, while the different mechanisms that could lead to an ’” or “”type profile are considered in section 6. In section 7 we deal with projection effects that are important for the profiles we examine and finally we give a summary of the paper by enumerating our conclusions in section 8.
2 The model
Since clear b/p bulges and especially X features are usually associated with strong bars we have chosen the “strong bar case” model of Skokos et al. (2002b) for our calculations. This model has the most massive bar with respect to the disc from all other models investigated in that paper. The general model consists of a Miyamoto disc, a Plummer bulge and a 3D Ferrers bar. It is a popular and extensively studied model for 3D bars, initially used by Pfenniger (1984). Here, we briefly describe the three components of the model to facilitate the reader of the current paper.
The potential of a Miyamoto disc (Miyamoto & Nagai 1975) reads:
(1) 
where is the total mass of the disc, and are the horizontal and vertical scale lengths, and is the gravitational constant.
For the bulge we used a Plummer sphere (Plummer, 1911) with potential:
(2) 
where is the scale length of the bulge and is its total mass.
Finally, the density of a triaxial Ferrers bar is:
(3) 
where
(4) 
, , are the semiaxes and is the mass of the bar component. The corresponding potential and the forces are given in Pfenniger (1984).
For the Miyamoto disc we use A=3 and B=1, and for the axes of the Ferrers bar we set , as in Pfenniger (1984) and Skokos et al. (2002b). The masses of the three components satisfy . The length unit is taken as 1 kpc, the time unit as 1 Myr and the mass unit as . If not otherwise indicated, the lengths mentioned all over the text and in the figures are in kpc. Units will be not repeated on the axes of the figures.
The total potential of our model is given by
(5) 
The 3D bar is rotating around its short axis. The axis is the intermediate and the axis the long one. The system is rotating with an angular speed and the Hamiltonian governing the motion of a testparticle can be written in the form:
(6) 
where and are the canonically conjugate momenta. We will hereafter denote the numerical value of the Hamiltonian by E and refer to it as the Jacobi constant or, more loosely, as the ‘energy’. The corresponding equations of motion are:
(7) 
In Table 1 we summarise the parameters of the specific model we will use in our paper.
GM  GM  GM  E(rILR)  E(vILR)  

0.72  0.2  0.08  0.4  0.054  0.467878  0.438317  6.31 

2.1 Periodic orbits and their stability
Most information about the dynamics of our Hamiltonian system is given by the study of the main families of p.o. that belong to the x1 tree (Skokos et al., 2002a). The phase space we study is structured by the p.o. and the way it is structured depends on their stability.
The space of section in the case of a 3D system is 4D. In order to find a p.o., the equations of motion are solved for a given value of the Hamiltonian, starting with initial conditions in the plane =0, for . The exact initial conditions for the periodic orbit are calculated using a Newton iterative method. A periodic orbit is found when the initial and final coordinates coincide with an accuracy at least 10. For the integration of the orbits we used a fourth order RungeKutta scheme. The relative error in the energy remains in our calculations less than 10.
For the estimation of the linear stability of a periodic orbit we first consider small deviations from its initial conditions, and then integrate the perturbed orbit to the next upward intersection. In this way a transformation is established, which relates the initial with the final point. The relation of the final deviations of this neighbouring orbit from the periodic one, with the initially introduced deviations can be written in vector form as: . Here is the final deviation, is the initial deviation and is a matrix, called the monodromy matrix. It can be shown that the characteristic equation is written in the form . Its solutions obey the relations and and for each pair we can write:
(8) 
where and
(9) 
The quantities and are called the stability indices. One of them is associated with radial and the other one with vertical perturbations. If , and , the 4 eigenvalues are on the unit circle and the periodic orbit is called ‘stable’. If , and , , or , , two eigenvalues are on the real axis and two on the unit circle, and the periodic orbit is called ‘simple unstable’. If , , and , all four eigenvalues are on the real axis, and the periodic orbit is called ‘double unstable’. Finally, means that all four eigenvalues are complex numbers but off the unit circle. In this case the orbit is characterised as ‘complex unstable’ (Contopoulos & Magnenat, 1985; Heggie, 1985; Pfenniger, 1985a, b; Zachilas, 1993). We use the symbols S, U, DU, to refer to , , and periodic orbits respectively.
The distinction of the various types of instability has been initially presented by Broucke (1969) and Hadjidemetriou (1975), and has been used in studies of the stability of periodic orbits in systems of three degrees of freedom. For an extended description the reader may refer to Pfenniger (1984) and Contopoulos & Magnenat (1985). For a general discussion of the kinds of instability encountered in Hamiltonian systems of N degrees of freedom the reader may refer to Skokos (2001).
The ‘stability diagram’ is a diagram that describes the stability of a family of periodic orbits in a given potential when one of the parameters of the system varies (e.g. the numerical value of the Hamiltonian E), while all other parameters remain constant (Contopoulos & Barbanis, 1985; Contopoulos & Magnenat, 1985). In other words it gives the variation of the stability indices and as the energy, or any other parameter, varies. The resulting curves are called the “stability curves”. Such a diagram helps in finding the transitions from stability to instability or from one to another kind of instability. Intersections, or tangencies, of a stability index with the axis are of special importance for the dynamics of the system. When this happens a new family is generated by bifurcation from the initial one and has the same multiplicity with it. The new family of p.o. inherits also the kind of stability that characterises the parent family. Intersections, or tangencies, of the stability curves with the axis, also generate new families but are accompanied by period doubling. The central family of the system is the well known x1 family (see e.g Contopoulos & Grosbøl, 1989), which is reduced to circular orbits on the equatorial plane in the axisymmetric part of the potential. In the full potential the members of this family are p.o. of elliptical shape, elongated along the major axis of the bar, remaining always on the equatorial plane. The “x1tree” is composed by the x1 family itself together with the 3D families, which are generated at specific SU or US transitions. At these transitions the stability index associated with the vertical perturbations has two intersections, or a tangency, with the axis. They happen at the regions of the vertical resonances (see e.g Patsis & Grosbøl, 1996). The “x1tree” is the backbone of b/p bulges and X features (Patsis et al., 2002). In standard models the first SU and US transitions due to an intersection of with the axis are encountered near the vertical 2:1 resonance and the bifurcated simple periodic 3D families of p.o. of the x1tree are introduced in pairs. Following the nomenclature of the fiducial case in Skokos et al. (2002a) these families are called x1v1 and x1v2 (associated with the vertical 2:1 resonance), x1v3 and x1v4 (associated with the vertical 3:1 resonance) etc. Finally we have to remark that at the transition to complex instability we have in general no bifurcating families of periodic orbits.
The stability diagram of our model that describes the evolution of stability for all families of the x1tree as E varies is given in Fig. 1a. The two stability curves, that exist already for E (left side of the figure), belong to x1. The stability index is associated with the radial, while with the vertical perturbation. The families x1v1, x1v2 etc. are introduced at the intersections of with the axis. In Fig. 1b we zoom into the ILR region, which is of special importance for the b/p dynamics. We will refer to it in detail in section 4. The straight line segments in Fig. 1a and b at for E , labelled “(x1v1)” is drawn to mark the area over which the x1v1 family is complex unstable. In this region the inidices and in Eq. 8 cannot be defined since we have in Eq. 9. The location of the four eigenvalues of the x1v1 p.o. with respect to the unit circle is as shown in Fig. 1c. We will refer to Fig. 1 throughout the text, describing the dynamical behaviour in the neighbourhood of p.o. participating in the building of the b/p bulge.
3 Where does the peanut start?
The radial ILR region of the model is found from the energy interval over which the families x2 and x3 (Contopoulos & Grosbøl, 1989) exist. E(rILR) in Table 1 is taken as the left limit of the x2x3 loop in the characteristic or in the stability diagram. E(vILR) is taken at the energy where we have the SU transition of x1 and the introduction in the system of the x1v1 family. The first (vertical) instability strip in the energies of x1 occurs at E (Fig. 1a). One more 3D family (x1v2) associated with the vertical ILR starts at the right border of this interval where we have the US transition (Fig. 1a).
It is generally believed that the peanut starts at the vertical 2:1 (vILR) resonance, where we encounter 3D p.o. with elliptical projections on the equatorial plane and “smiles” or “frowns” in their sideon views, i.e. p.o. of the x1v1 family (Pfenniger, 1984; Combes et al., 1990; Patsis et al., 2002; Bureau et al., 2006; Patsis & Xilouris, 2006; MartinezValpuesta et al., 2006). We note that Combes et al. (1990) emphasise that at the end of their body models the radial and vertical 2:1 resonances coincide.
The bar is elongated along the yaxis and its backbone is the x1tree. The members of the x1v1 family are not the first 3D orbits in the system in general. There are other 3D families of p.o. in the very central part of the model associated with the bulge component of the potential, as well as 3D nonperiodic orbits that support the disc in the same region. Among them there are 3D quasiperiodic orbits with E associated with the planar x1 family. We investigate here their contribution to the 3D shape of the bar in the central region of the model by examining their effect in the resulting sideon profiles.
For E , the only existing barsupporting family of p.o. belonging to the x1tree is x1, which in this energy range is stable everywhere. Albeit planar, x1 contributes to the building of the bar away from the equatorial plane. Due to the stability of x1, quasiperiodic orbits around them will be trapped in “rotational” or “tube” tori (Vrahatis et al., 1997; Katsanikas & Patsis, 2011), which are 3D objects embedded in the 6D space. Perturbations of the initial conditions of the x1 p.o. in the or direction will result to 3D orbits extending in the 6D phase space.
Interesting is the shape of the associated with these tori orbits in the configuration space. We followed their morphological evolution as the perturbations of the initial conditions of the p.o. at a certain energy increase and also as the energy increases approaching the critical value, where we have the SU transition of x1 at the vertical 2:1 resonance (E ).
In Fig. 2a we depict the projections of two quasi periodic orbits around x1 at E . Both of them have the same initial condition as the x1 p.o. at this energy . The narrow, darker one, is perturbed in the z direction by , while the thicker one by . The projections of these two orbits have similar shapes. In Fig. 2a the thick orbit has a maximum thickness about 0.2. The energy of the two orbits is the same, thus the orbit with the higher has a shorter extent along the direction. In the 3D configuration space these orbits have a toroidal morphology. In the (y,z) projections they are characterised by a minimum for and a maximum z at maximum as in Fig. 2a. This is a typical, general morphology of the quasiperiodic orbits around x1 away from the ILRs of the system.
However, as energy increases we observe that initially for small perturbations, the minimum for becomes a local maximum, as in the case of the (y,z) projection of the orbit in Fig. 2b. In this particular case we have an orbit with E having initial conditions that deviate from those of the x1 p.o by =0.01. In order to clearly observe the formation of the local maximum, the scale of the zaxis is 0.02 times that of the yaxis.
By studying the quasiperiodic orbits around x1 for energies just before the SU transition of x1 at the vertical 2:1 resonance, we realise that the local maximum at becomes prominent for quasiperiodic orbits in the neighbourhood of x1. The orbit in Fig. 2c is for E = and occurs by perturbing the x1 p.o. with initial conditions =(0.0930694109,0,0,0) by =0.01. The projection of the orbit has taken a shape like an “” symbol. The p.o. x1 at this energy is still stable (, ). Despite the fact that we are before the SU transition at which the x1v1 family is bifurcated, the shape of this orbit is typical of an orbital structure expected to appear in the neighbourhood of the stable x1v1 family and its symmetric with respect to the equatorial plane (x1v1) (Patsis et al., 2002). There is a connection between the evolution of the border of the quasiperiodic orbits and the shape of the orbits of the family of p.o. to be bifurcated. In order to show this we plot in Fig. 2c with red the projections of the x1v1 and x1v1 p.o. for a slightly larger energy (E ) when they have been introduced in the system as stable p.o. the SU transition of x1.
On account of the stability of the x1 family, 3D quasiperiodic orbits are expected to be trapped around their p.o. in invariant tori called “tube” or “rotational”. (Vrahatis et al., 1997; Katsanikas & Patsis, 2011). By means of the colourrotation method (Patsis & Zachilas, 1994), these tori appear in the case of family x1, before the vILR, as typical “tube tori”. This happens because there is a 3D projection, , in which the torus appears as folded and selfintersecting (Vrahatis et al., 1997). In Fig. 3a we depict the perturbed x1 orbit by at E . For constructing this 4D object the orbit is integrated for 500 periods of the x1 p.o. at the same energy. We use the projection and the consequents are coloured according to their values, normalised in this case in the interval [0,1]. The colour bar, at the right side of the figure, gives the colour variation. The point of view in spherical coordinates is defined by the angles and the distance of the observer from the centre of the figure equals the largest dimension of the box surrounding the drawn orbit (for more details see Katsanikas & Patsis (2011)). The shape of the object resembles an “8” or “” symbol. At the region of the intersection of the two meeting branches we encounter colour values at the two opposite edges of the spectrum, which means that it is a pseudointersection. The figure gives a typical 4D Poincaré surface of section for all 3D barsupporting quasiperiodic orbits before the first SU transition of x1. A 3D view of this orbit in the configuration space is given in Fig. 3b. The point of view in this figure is . The quasiperiodic orbit is drawn with black. The red curve is the x1 p.o. on the plane at the same energy. One can easily understand that by viewing this orbit sideon, the shape will be as the one depicted in Fig. 2c.
The study of successive projections shows that the shape of the quasiperiodic orbits with the same deviation from the initial conditions of x1 changes monotonically as the energy varies. In all examined models we found particles in orbits supporting a b/p structure at the region of the ILRs and before the introduction of the x1v1 family in the system. Based on the appearance of the local maximum at in the projections, we can say that we have peanutsupporting orbits already at a distance of about 0.6 kpc along the major axis of the bar.
4 The region of the radial and vertical ILRs
4.1 Nonperiodic orbits in the neighbourhood of x1
Thus far we explored the contribution of the quasiperiodic orbits trapped by x1 to the formation of a peanut structure before the introduction in the system of the x1v1 family. The x1 family is introduced in the system at the 2:1 vertical resonance. The SU transition of x1, at E marks the beginning of the vertical ILR region, while in our model the radial ILR already exists at this energy. Beyond E the x1 family becomes simple unstable, with , for E . Let us call this energy interval E. At the end of the E interval we have the US transition of x1. Within E the family x1 is simple unstable (U), while x1v1 is stable (S). We find also the two known planar families of p.o., x2 as stable and x3 as simple unstable. The energy range over which the x2, x3 families extend is significantly larger than the E region. We get an overview of the interconnections of the various families of p.o. in the region of the radial and vertical ILR by zooming into the evolution of the stability indices of all implicated families (Fig. 1b). We note that in the scale we use for presenting Fig. 1b (let alone Fig. 1a) the index of x1v1 appears coinciding with the axis. In reality just remains always very close to the axis, having two tangencies at E and at E . At E meets and the x1v1 family becomes complex unstable. The stability indices associated with the x2x3 loop are plotted with red.
In E x1 is unstable relative to vertical perturbations, because . However, for perturbations on the equatorial plane x1 remains stable, because . A cross section in E with orbits deviating from x1 by or , and thus remaining on the equatorial plane, will be as the one in the example depicted in Fig. 4 for E . We note that despite the fact that for this E both x1 and x3 p.o. are characterised as simple unstable (U), the phase space structure in their neighbourhood, as expected, has different properties, because the two eigenvalues on the unit circle in the former case are associated with the radial perturbations, while in the latter with the vertical ones (see section 4.2 below).
As regards the dynamics in the neighbourhood of x1 away from the equatorial plane, for energies close beyond the bifurcating point, as is E (Fig. 5), we find the following: The vertically perturbed x1 orbits are characterised by slow diffusion in phase space. The reason is the presence of the invariant tori belonging to the stable family x1v1, which has been already bifurcated from x1 at E . Despite the fact that the presence of rotational or tube tori cannot restrict a region in the 4D space of section, their proximity to the initial conditions of the U x1 p.o. hinders considerably the diffusion of the consequents of perturbed x1 orbits in the available phase space. There is a range of small deviations from the initial conditions of x1, in and in directions, for which the integrated orbits will be orbits on the invariant tori of the bifurcated stable x1v1 p.o. at the same energy. This is due to the fact that close beyond the critical energy of the bifurcation of x1v1 from x1, the bifurcated family has almost the same as x1 and a small . In Fig. 5a we give the projection of the consequents of the x1 orbits perturbed successively by (red curves) and (blue curve), at E =0.438225. The location of x1 is at . The orbits on the red curves, for to , surround the initial conditions of the p.o. x1v1 (right) and x1v1 (left). The location of these p.o. is indicated with arrows (the plotted dots correspond to tiny tori around them). The orbit with surrounds in the plane all three p.o. In the 3D projection we observe (Fig. 5b) that all these orbits are on tori, the thickness of which (dimension) increases with increasing . The sideon morphology of the orbits on the red curves is as the one in Fig. 5c () for the closed curves around x1v1 and symmetric to them with respect to the equatorial plane for the closed curves around x1v1. The tori which surround all three of the p.o., like the one with , have a clear peanut structure, similar to the one we have seen in Fig. 2c, at an energy before the SU transition of x1. This is given in Fig. 5d.
If we include in the vertical (by , ) perturbations of x1 also or displacements, i.e. if we essentially perturb vertically quasiperiodic orbits on closed curves around x1 in Fig. 4, the situation does not change considerably. For the region inside the fourth closed curve around x1 in Fig. 4, there is again a range of vertical perturbations for which our orbits abut on x1v1 invariant tori and the orbits in the configuration space have a morphology similar to the orbits in Figs. 2b,c. For example, orbits perturbed by starting in this region belong to x1v1 invariant tori. Differentiations from this behaviour is observed when we start integrating orbits with or starting at and beyond the fourth closed curve around x1 in the x1 “stability island” of Fig. 4. There, for the same range of or perturbations, for which we reached x1v1 tori before, we encounter orbits with complicated morphologies that are not similar among themselves as the perturbation varies. This indicates that we reached the limits of the phasespace region occupied by the x1v1 tori in the 4D space of section.
For instance, in Fig. 6 we show two such orbits with initial conditions corresponding to perturbations of the orbit on the fourth invariant curve around x1 in Fig. 4, perturbed by (a) and by (b). The sideon profiles of these two orbits are not similar among themselves. However, both of them could support a b/p morphology of a central bulge. On one hand the sideon profile of the orbit in Fig. 6a reproduces a morphology resembling orbits trapped by either branches of x1v1 (x1v1 and x1v1), which are symmetric with respect to the equatorial plane. However, it also differs from the orbit in Fig. 5d. On the other hand the orbit in Fig. 6b has a x1 quasiperiodic orbit morphology similar to the one in Fig. 2a even if we continue its integration for dynamical times. In Patsis & Katsanikas (2014) (paper II), we will return to these orbits discussing orbits that combine boxy sideon and faceon profiles. For the time being we keep in mind that there is an extended volume around x1 orbits, even in the energy interval where it is simple unstable, that contributes to the boxiness of the sideon profile of the model.
Beyond the fourth drawn elliptical curve around x1 in Fig. 4, we enter a region where we have chains of islands of stability. This region on the plane has been proven to be of special importance for the dynamics of the b/p structures (see also paper II). Orbits with origin in this region may have similar sideon morphologies with the orbit in Fig. 6b for thousands of dynamical times. Their sticky character is revealed very illustrative by means of the colourrotation method. The 3D projections of the space of section show us that these orbits remain confined for large times (of the order of to dynamical times) in certain regions of the phase space. We remind that usually in Galactic Dynamics we are interested about timescales of tens or in extreme cases a few hundreds of dynamical times.
As an example, let us consider the orbit with initial conditions and integrate it for dynamical times. In Fig. 7a we plot this orbit with blue points in the projection on the surface of section presented in Fig. 4. Fig. 7a includes the invariant curves around x1 and x2 and the chaotic zone around both of them (black points). The points coloured red belong to the corresponding projection of the orbit in Fig. 6b. By careful inspection of Fig. 7a one can observe that the blue consequents build at least two thin, relatively denser layers with respect to their surrounding consequents in phase space. One of them is found immediately after the “red” orbit and an even thinner one further out. We indicate these structures with black arrows. The orbit is trapped for some time on these two layers before it expands and fills a larger region that seems to be confined by a chain of islands just before we enter the chaotic zone. By considering the projection (Fig. 7b) we realise that the two orbits (“red” and “blue”) build in fact cylindrical structures that remain confined also in the zdirection. This is more clear for the consequents of the “red” orbit, which need an extremely large integration time to diffuse in phase space (of the order of several Hubble times). Nevertheless the chaotic character of both orbits is evident in their 4D representations by means of the colourrotation method. They are given in Fig. 7c and Fig. 7d for the “red” and “blue” orbits respectively. In Fig. 7c we use as spatial coordinates and we colour the consequents according to their values for the “red” orbit. In Fig. 7d we use as spatial coordinates and we colour the consequents according to their values for the “blue” orbit. The colour bars giving the colours of the consequents are given on the right side of these figures. Thus, the blue and red consequents are not related to the colours we use in Figs. 7a, b for distinguishing the two orbits. The first feature to be observed is that in both Figs. 7c,d we have mixing of colours on the structures formed in the 3D projections. All 4D combinations of spatial coordinates and colours give qualitatively similar figures with those presented in Figs. 7c,d for each one of the two orbits. In all 3D projections we have namely a toroidal object with colour mixing. In Fig. 7d we observe that the “blue” orbit forms an inner toroidal structure before it expands, after about 5500 dynamical times, occupying a larger volume. The inner torus is tangent to the one depicted in Fig. 7c as their projections in Figs. 7a indicate. The part of the “blue” orbit corresponding to the inner torus is morphologically very close to the one presented in Fig. 6b.
In conclusion we can say that there is a plethora of sticky orbits in the energy interval in which x1 is simple unstable, that can reinforce boxy sideon profiles. In 3D projections of the 4D spaces of section, which include the plane, these orbits are found trapped in cylindrical structures between the closed curves around x1. This is an important dynamical mechanism that has to be taken into account for the building of the peanut.
4.2 The x2, x3 region
Apart from x1, the other simple unstable 2D family in the E region is x3. Unlike x1 it is radially unstable. On the plane x3 is found inside a chaotic zone, which surrounds the closed curves around x1 and the stability island around x2 (Fig. 4). If we perturb by the initial conditions of x3 we find in its neighbourhood orbits confined in the zdirection of the configurations space. However, they are not as sharply defined as the quasiperiodic orbit around x2. For comparison with x2, we give the morphology of two perturbed by orbits, i.e. of the stable x2 and of the simple unstable x3 orbit within 100 dynamical times in Fig. 8a,b respectively. Note that in order to understand the internal details of these orbits the scales of the axes in the (x,z) as well as in the (y,z) projections are not equal. The perturbed x3 orbits spend more time following a x2flow (elongated perpendicularly to the bar) but during their excursions in phase space away from the x2 island of stability their projections on the equatorial plane precess. As a result they fill a disky region around a box that extends along the xaxis (Fig. 8b, left panel). At any rate these orbits, like those in Fig. 6, remain close to the equatorial plane and support a boxy structure in the ILR region of the model. Despite the difference in their stability character, the contribution of nonperiodic orbits around x3 to the local density away from the equatorial plane is comparable with that of quasiperiodic orbits around x2. However, they have a clearly different morphology in their faceon views.
Unlike in the case of the x1 family, diffusion from the immediate neighbourhood of x3 is not hindered by x1v1 tori. This is expected, because x3 is not the parent family of x1v1 and so the tori of x1v1 are always at a distance from the x3 p.o.
The different dynamical behaviours in the neighbourhood of the two 2D simple unstable p.o. x1 and x3 is reflected in the different structures observed in the 4D surfaces of section of the perturbed orbits presented in Fig. 8. In order to get a clear view of these structures in the 4D cross section of the phase space we have integrated them for a period 80 times larger than the one used to depict them in the configuration space. Nevertheless, no new dynamical features are introduced in the system during this much larger period. In Fig. 9a we apply the colourrotation method to visualise the 4D space. In this particular case we use the projection and we colour the consequents according to their variation, normalised in the [0,1] interval. The x1 orbit is perturbed by . We observe that its 4D cross section is of the form of a ribbon on which colour has a smooth variation according to the colour bar that is given at the right side of the figure. This is a typical fingerprint of the 4D space of section in the neighbourhood of a p.o. with this type of instability for times larger even than a Hubble time (cf with figure 3 in Katsanikas et al., 2013) On the other hand the corresponding 4D space of section of the perturbed x3 simple unstable orbit is characterised in its projection by two conglomerations of consequents that are of toroidal form, as well as by a number of consequents between them that form a bridge. This can be observed in Figs. 9b,c and d. In Fig. 9b we can see that the colour on the front (left) toroidal structure has almost a single shade, reflecting the proximity of the values in this region. The same is true for the toroidal formation of the consequents at the back (right) side of the cross section. This is less discernible in Fig. 9b due to the presence of the consequents of the “bridge”, which are projected on it and have mixed colours. For this reason, in order to better understand the geometry of Fig. 9b, we observe it from two more viewing angles in Fig. 9c,d. This allows us to better understand the toroidal character of the two dense regions (Fig. 9c) and the bridge between them (Fig. 9d). In Fig. 9c, the left toroidal structure of Fig. 9a is projected as a dense ring, being actually in front of another, less dense, toroidal structure. We observe in the centre of the figure an empty region, which corresponds to the hole of the two toroidal objects. In Fig. 9d the dense regions (left and right) are the “tori” while the consequents between them belong to the bridge. We conclude that the phase space in the neighbourhood of the two simple unstable 2D p.o. x1 and x3 is qualitatively very different. It is determined by the presence of other families of p.o. coexisting in the same energy. In the case of x3 we have a “type1” (tangent) bifurcation (Contopoulos, 2004, section 2.4.3) with the stable counterpart being the 2D family x2, while in the case of x1 we have a “type 3”, direct bifurcation of the 3D families x1v1 and x1v1.
Despite the small size of the orbits investigated in this section, compared with the dimensions of the b/p bulge, the orbits in E give an illustrative example of how the dynamical phenomenon of stickiness becomes essential for building the peanut as E increases approaching corotation. Similar phenomena take place at larger E as will be described below. Stickiness is in many cases a significant dynamical phenomenon in the neighbourhood of simple unstable p.o. A detailed study of the sticky chaotic orbits just beyond the first transition is given by Katsanikas & Patsis (2011) and by Katsanikas et al. (2013) in a double MiyamotoNagai rotating potential. Here we find similar results and this indicates a general dynamical behaviour in the phase space at the vertical ILR region of rotating galactic potentials.
5 x1v1 peanuts
Up to now we examined the starting point and the central structure of the peanut. However, our main goal is to study the overall structure of b/p bulges and the X feature. For this, most important are the outermost orbits that determine to a large degree the observed morphology. As mentioned in the introduction, the orbits that have been proposed by most authors of relevant studies as building blocks belong to the stable family of p.o. that is bifurcated at the vertical 2:1 resonance, i.e. the x1v1 family. We want to investigate in what degree x1v1 is sufficient in building peanut and Xshaped bulges with dimensions similar to those observed in edgeon disc galaxies. An obvious obstacle with this family is its complex instability character over large energy intervals. This is a typical property of this family in 3D rotating galactic potentials (Pfenniger, 1985b; Patsis & Zachilas, 1990; Skokos et al., 2002b; Katsanikas et al., 2011a).
In our model the longest projection of x1tree p.o on the major axis of the bar
(yaxis) reaches a radius . This can be roughly considered
as the length of the bar that can be built by the orbits of our model. The bar
extends to 1.5 horizontal scale lengths of the Miyamoto disc. This is a value
close to the mean bar length as fraction of the exponential disc scale length
proposed by Erwin (2005) for early type barred galaxies
The maximum length of the b/p bulge in the sideon projection is expected to be according to Lütticke et al. (2000) (see Sect. 1). Then, by measuring the dimensions of the X features in Aronica (2006) we estimate a height compatible with observations about . This height is reached by x1v1 p.o. at E =. Taking this into account we have built indicative profiles with nonperiodic orbits in the neighbourhood of the x1v1 p.o. Such sideon, , profiles are given in Fig. 10. They have been constructed by taking eight orbits at equally spaced energy intervals between E and integrating them for time equal to 12 periods of the x1v1 p.o. at the energy E =. The profiles have been converted to images by means of the ESOMidas software. The colour bar in the middle of the figure indicates increasing density from left to right. The criterion used for selecting these orbits in each panel differs. In Fig. 10a we have taken perturbed x1v1 orbits by . The perturbations were such as to secure a clearly by eye discernible “” feature formed when we plot the sideon views of the nonperiodic orbits and their symmetric with respect to the equatorial plane of the model. This perturbation is less than 0.1 pc and is added to the initial condition of each p.o. x1v1. The rest of the initial conditions of the used orbits are 0, as for the corresponding x1v1 p.o. The orbital profile in Fig. 10b is constructed by considering perturbations by . This time we obtain “” morphologies for perturbations of the individual orbits in the range . In Figs. 10c,d we start our orbits with =0.05 and =0.1 respectively, while the initial component for orbits of the x1v1 family is 0. In all cases the individual orbits build in their (sideon) projections a “”type profile. However, we have to stress that, like in the case of the profiles with periodic orbits presented in Patsis et al. (2002) and Patsis & Xilouris (2006), the wings of the X in the profiles in Fig. 10 are not along the wings of the individual orbits. In other words by considering all orbits together we do not form a single “” structure. The wings of the X features in the profiles are formed along the zmaxima of the orbits, where (cf with figure 19a in Patsis et al. (2002)). This will become apparent also in the following section, Sect. 5.1. For the time being we conclude that x1v1 p.o. perturbed by or by , with initial form a sharper X profile than the orbits starting away from the equatorial plane . We underline the fact that by constructing these profiles we have applied the perturbations to the x1v1 p.o. independently of their stability.
5.1 Contribution by orbits close to complex instability
Qualitatively, the orbits that participate in the building of “x1v1based” X features reaching heights up to are of two kinds. They are either quasiperiodic orbits trapped by stable x1v1 p.o. or chaotic orbits associated with complex instability. The former are encountered at energies E and E , while the latter at E . The quasiperiodic orbits shape the innermost, as well as the outermost, and thus most important for the overall morphology, part of the peanut. On the other hand, over the larger part of the E interval, the x1v1 family is complex unstable. Nevertheless, our investigation shows that even if we integrate the chaotic orbits participating in Fig. 10 for times equal to 40 periods of the outermost x1v1 p.o. participating in the X structure the boxy morphology is not destroyed.
In the case of complex instability the quantity (Eq. 9) is . We have observed that in our model, for a constant deviation from the initial conditions of the periodic orbit, smaller values of correspond to orbits more extended in configuration space. Fig. 11 shows the individual orbits used to construct the profile of Fig. 10c. In the central panel of Fig. 11 we plot as E varies. From “1” to “8” we mark in this diagram with heavy dots the locations of the orbits we used. The projections of the corresponding p.o. are given in the embedded frame. With black we have drawn the stable ones, while the orange coloured are complex unstable (). The arrows help us understand that the wings of the X feature are formed along the maxima of the successive orbits. In the frames around the central figure, the orbits numbered “1”, “7” and “8” are quasiperiodic and have morphologies resembling those of the x1v1 p.o. Counterintuitively, we observe at some energies that also chaotic orbits in the neighbourhood of complex unstable x1v1 p.o. resemble “”shaped structures (orbits numbered “2”, “4”, “5”, “6”). However, they are much less confined in the configuration space. The orbit in position “3” for E is a particular case, since it is more affected by the presence of the family x1v2, which for E is simple unstable. We will refer to it in detail below (Sec. 6.2). In Fig. 11 we observe that the most extended orbit in the configuration space is orbit “5”, which has the smallest . The same happens with orbits perturbed by and from x1v1 (not shown in Fig. 11).
The explanation for the orbital behaviour of the “” orbits in Fig. 10 is again related to the phenomenon of stickiness. Stickiness associated with complex instability has been recently studied extensively by Katsanikas et al. (2011a). We find in our barred model perfect agreement with the results of that study. The consequents in the 3D projections of the 4D space of section in the neighbourhood of the complex unstable p.o. form spirals (Contopoulos et al., 1994; Papadaki et al., 1995). In the 4D representation by means of the colourrotation method we realise that there is a smooth colour variation along the arms of these spirals. For small perturbations of the initial conditions objects that are called confined tori are formed in the neighbourhood of the “” p.o. in the 3D projections of the spaces of section (Pfenniger, 1985a, b; Jorba & Olle, 2004; Olle et al., 2000). Like in Katsanikas et al. (2011a) we find also in the model of the present paper that these are 4D objects with an internal spiral structure and with a smooth colour variation on their surfaces. For explaining the peanut and the X structure, we are interested in orbits that give a few tens of consequents on the 4D surfaces of section. These are consequents that describe a first set of spirals on the confined tori. Especially the first 1020 consequents are located in parts of the spirals very close to the p.o. During this time the orbital behaviour in the neighbourhood of the “” p.o. is characterised by stickiness (Katsanikas et al., 2011a). This results to morphologies having a certain similarity to quasiperiodic orbits around a stable orbit. This explains the morphologies of the orbits “2”,“4”,“5” and “6” in Fig. 11. Diffusion in phase space happens after completing thousands of consequents. However, for larger perturbations or for larger energies of the orbits, diffusion happens earlier.
6 “OX” vs “CX” profiles
The sideon profiles given in Fig. 10 resemble the kind of X features encountered in b/p bulges, which have wings that do not cross the centre of the galaxy. A typical, example of this kind of b/p bulges is the galaxy NGC 4710 (for a highresolution image of this object, see the Space Telescope Science Institute’s news release STScI200930). This is the kind of X features labelled by Aronica (2006) and Bureau et al. (2006) as “OX” (offcentred). Indeed, in many images of galaxies the direction of the wings of the X, as revealed after applying an unsharpmasking filter on them, clearly shows that their extensions towards the centre of the galaxy will not cross it (see also figure 2 in Patsis & Xilouris, 2006). However, in some other cases we have a genuine X feature, as is for example the edgeon galaxy HCG87A in the Hickson 87 group (Hickson, 1993) (for a high resolution image, see the STScI199931a news release). Aronica (2006) and Bureau et al. (2006) call this type of X, “CX” (centred). The sideon “OX” profiles in Fig. 10 are based on orbits in the neighbourhood of x1v1 p.o. Possible projection effects, that could bring the wings of the “x1v1X” closer and form “CX” instead of ’OX” profiles will be discussed in section 7. Below we examine the possibility of building true “CX” profiles using orbits.
6.1 Possible x1v2 contribution
In nearly sideon views of galaxies, true X’s can be formed only with the help of orbits associated with other families than x1v1. Orbits belonging to the x1tree, bifurcated closer to corotation, with the appropriate shape may have in some models stable parts. Such an example is given by Patsis et al. (2002) for their standard model, based on the x1v4 family. Families bifurcated closer to corotation, if populated, will build b/p bulges in which the ratio of the length of the bar to the length of the b/p distortion is slightly larger than 1. However, in the present study we are looking for b/p structures in which this ratio is , as we mentioned in the introduction.
Under those considerations, a family candidate to play this role would be x1v2, bifurcated from x1 at a US transition, in our model at E (see Fig. 1). The sideon projections of the p.o. of this family have an type morphology (Skokos et al., 2002a). In Fig 12a we plot a typical sideon view of a x1v2 p.o. orbit (red) together with the sideon projections of the p.o. x1v1 and x1v1. Obviously the x1v2 orbits have qualitatively the appropriate shape to support a central type feature in the plane. However, the x1v2 family is initially simple unstable, inheriting the simple unstable character of x1, while for larger E it evolves to double unstable (Fig. 1a). Thus, the mechanism of trapping quasiperiodic orbits around x1v2 p.o. cannot be applied.
Nevertheless we may look for nonperiodic orbits that will fill the empty central parts in the profiles of Fig. 10 converting them from “OX” to “CX”type. The range of energies that has to be covered is E , i.e. roughly the interval within which x1v2 is simple unstable. When this family becomes “DU” we do not expect any contribution to an observed structure, since double instability leads to strong chaoticity (Katsanikas et al., 2013).
We know (Katsanikas et al., 2013) that in a range of E away from the immediate
neighbourhood of the bifurcating point
For energies beyond the bifurcation point, we have integrated x1v2 orbits perturbed by or for equal time as for the orbits composing the profiles of Fig. 10. We first realised that such orbits, with E , have the appropriate sizes to fill the central empty regions of the profiles in Fig. 10. We also found that a type morphology is supported for the longest time by orbits perturbed by . Perturbed x1v2 orbits by are able to reinforce to some degree an type morphology in their sideon projections. For example, by substituting the orbits used for E and in Fig. 10b by x1v2 orbits perturbed by and 0.05 respectively, we obtain the profile given in Fig. 12b. Despite a small asymmetry with respect to the y=0 axis we can observe a “CX” central feature.
By varying E , we realise that the ability of perturbed x1v2 orbits to support an type sideon profile depends in general mostly on the overall structure of the phase space in the region close to the p.o. Critical for this purpose is the relative location of the x1v1 and x1v1 rotational tori with respect to the x1v2 p.o. For energies close beyond the US transition of x1 (E 0.425) the x1v1 tori are still in the neighbourhood of the p.o. x1v2. At E = 0.42398, the stability index of x1v2, which is associated with vertical perturbations, is . Its initial conditions are , while the corresponding initial conditions for x1v1 at the same energy are . Despite the fact that the orbits explore a 4D space, the proximity of the x1v1 tori to the location of the x1v2 p.o. guides the asymptotic curves of the x1v2 manifolds. As long as the unstable manifold of the U p.o. winds around the x1v1 tori, the chaotic orbits follow a fuzzy x1v2 morphology. By this we mean that during the time of integration the orbit has conspicuous x1v2 morphological elements for considerable time intervals. For example, the orbit of Fig. 13a has a hybrid morphology with features mainly from x1v2 and x1v1. We consider as hybrid orbits supporting a “CX” profile, orbits that create a clearly discernible intersection point close to the equatorial plane, i.e. at , as is the orbit in Fig. 13a. We find that about the first 20 consequents of such orbits in the plane are projected very close to the unstable manifold of x1v2. For the orbit of Fig. 13a this can be seen in Fig. 13b, where we give the projection of the unstable manifold (red curve) of x1v2 and the corresponding projections of the first 18 consequents of the orbit (black points). The location of the x1v2 p.o. is marked with an symbol. The red curve corresponds to the part of the manifold that winds around the x1v1 and x1v1 tori.
As the energy increases the distance between the p.o. x1v2 and that of the x1v1, x1v1 tori increases and the diffusion of the chaotic orbits from the neighbourhood of x1v2 into the available phase space is fast (Katsanikas et al., 2013).
However, we find anew an enhancement of x1v2 type morphologies at energies E . The new element that has to be taken into account in this energy interval is the introduction in the system of the planar stable family o1. This family is introduced in the system together with its symmetric with respect to the yaxis family, o2 (Skokos et al., 2002b). In Fig. 1 we can see how o1 is introduced in the system at E , when the stability index of x1 becomes . 3D chaotic orbits, sticky to o1, o2 tori have morphologies resembling that of Fig. 13a.
Our analysis shows that there are orbits starting in the neighbourhood of x1v2 that can contribute to the formation of a “CX”type profile in the sideon views of the models. We find them being sticky to tori of stable families coexisting at the same energy. Then these orbits can be reached also by perturbing the initial conditions of the stable orbits themselves, like the one labelled with “3” in Fig. 11, which we found by perturbing the x1v1 initial conditions. However, the orbits reinforcing a “CX”type profile do not share common features in the two other projections of the configuration space. This underlines the chaotic character of the orbits.
6.2 The relative importance of nonperiodic orbits and the x1mul2 family
Trying to assess the significance of the presence in phasespace of each of the families of p.o. for enhancing a “CX” sideon profile, we investigated perturbations of their common “mother family”, x1, away of the equatorial plane . For this purpose we considered first x1 orbits on the equatorial plane and we explored the phase space in the vicinity of x1 p.o. perturbed by . This gave us insight about the motion of particles, which are following initially a x1 flow and are kicked out of the galactic plane. A typical energy, where we can follow the influence of all these families to the orbital dynamics is E = . It is typical because it is away from bifurcating points and the families of p.o. retain their main stability character. At E = coexist the families x1 (S), x1v1 (S) and x1v2 (U) (Fig. 1). For describing the properties of these orbits which are relevant to the sideon morphology of the peanut structure, we base our description on the projection of the 4D phase space. This is given in Fig. 14. In this figure the x1 p.o. is located at (0,0) marked with a heavy black dot.
At E small perturbations of x1 by bring in the system quasiperiodic orbits of cylindrical shape. They are characterised by a local minimum for in the projections, similar to that of quasiperiodic orbits around x1 (Fig. 2a). We notice though, that at this E we have a second, local, minimum at , best observed in the projection. This means that we have two characteristic sets of dents along the boundaries of the quasiperiodic orbits. Such an orbit is depicted in Fig. 15a and the local minima under discussion for , are marked with arrows. The innermost two elliptical curves around in Fig. 14 correspond to this kind of orbits. The symbols along this and other orbits in Fig. 14 correspond to the consequents of parts of the orbits that are used for presenting their morphologies (see below).
As we increase the perturbation, always for constant E , both local minima deepen, approaching . The increase of the perturbation brings us closer to the initial conditions of x1v2, which is marked in Fig. 14 with a red “” symbol close to (0,0.2). With a green “” symbol close to we mark the x1v2 orbit. Essentially these two orbits are two branches of the same family (Skokos et al., 2002a). We find that starting from x1 and approaching the initial conditions of x1v2 by increasing , the boundaries of the nonperiodic orbits become more and more similar to the shape of the x1v2 p.o. This is similar to what was happening as we were approaching the bifurcation point of x1v1 by increasing the energy in Section 3.
By reaching the perturbed by x1 orbit we encounter a new kind of orbital behaviour. The first 1500 consequents form an elliptical, lemonshaped, ring in the projection, plotted with dark blue colour in Fig. 14. Then the projected consequents drift to the inner part of the ring and the ring thickens. After 2249 intersections with the surface of section, the consequents depart from the area of the ring and fill a region around two “empty” lobes in the left and the right side of Fig. 14. In the interior of these two lobes are projected the initial conditions of x1v1 (right) and x1v1 (left), marked with black and grey “*” symbols respectively. The first 50 consequents of the perturbed orbit after their departure from the ring area are plotted with intermediate sized blue dots in Fig. 14 and can be seen mainly around the left lobe. The morphology of this orbit within 10 periods of x1 (corresponding to the consequents marked with yellow “” symbols on the blue ring) is given in Figs. 15b,c and d for the , and projections respectively. There is a conspicuous enhancement of an type morphology in the view. If we consider the orbit for the whole time it remains on the blue ring, this morphology persists. This can be seen in its projection given in Fig. 15e. We underline also the narrowness of the orbit along the axis. This may be crucial for the appearance of the “CX” morphology in the edgeon views of disc galaxies, since it secures a small distance between the front and the back branches of the  shaped orbit.
By increasing further the perturbation in the direction, we come very close to the initial conditions of x1v2. However, by varying only the coordinate we will never reach the x1v2 p.o. since we keep constant the initial coordinate, equal with the value of x1. This is smaller than the x1v2 initial condition by about 0.011. Orbits in the neighbourhood of the p.o. x1v2 recede away from it following trajectories in the direction specified by its unstable maniflods (Fig. 13b). The “empty” lobes of Fig. 14 are like those given in Fig. 13b. For example, small displacements from x1v2 (red ) in the direction lead to clouds of points around the lobes in Fig. 14. The light blue and magenta coloured points belong to orbits perturbed from x1 by and respectively. The empty in Fig. 14 regions are occupied by tori around x1v1 and x1v1 and are needed thousands of consequents to fill them in the projection. In Fig. 15f we give the projection of the first of these two orbits integrated for 100 x1 periods. On top of it we colour with red the part that corresponds to 10 x1 periods. The consequents of this part of the orbit are given with light blue “” symbols around the right “empty” lobe in Fig. 14. This additional information of the shape of the orbits within realistic times is necessary in order to understand the morphologies they support. Nevertheless, for understanding the dynamical phenomena behind the observed orbital dynamics we frequently need to extend our calculations to times of the order of a Hubble time or more.
Orbits of x1 perturbed by larger than the x1v2 initial condition fill roughly a ring surrounding the area defined by the initial conditions of x1v2 and x1v2 and the “empty” lobes to the right of x1v1 and to the left of x1v1 in Fig. 14. In this area exist four tiny stability islands. Starting close to we find them in azimuthal distances of about (coloured red). These islands are the projections of rotational tori around a stable p.o. of multiplicity four. This orbit enhances the boxiness of the peanut mainly by forming a zone of stickiness connecting the four tori. Orbits with an ab initio strong chaotic behaviour are found for . The morphology of the perturbed by x1 orbit is already quite chaotic in the projection within 10 x1 periods, as can be observed in Fig. 15g. The consequents corresponding to the part of the orbit depicted in Fig. 15g are given with heavy, this time black, dots in Fig. 14. They are distributed in a zone that surrounds the orbits sticky to x1v1 and x1v2. The same holds also for the next 1000 consequents (drawn with lighter black dots). We observe however, that there is no substantial mixing of the outer chaotic zone, which includes the black dots, with the inner structure containing the lemonshaped ring, the x1v2 and x1v2 p.o. and the “empty” lobes around them. Stickiness is ubiquitous in Fig. 14. We practically reach the same topology of phase space by perturbing the x1 orbits in the direction and by perturbing the x1v2 orbit itself. One can compare Fig. 14 with Fig. 13b but also with figure 5 in Katsanikas et al. (2013). The basic feature that structures the phase space is the presence of the area occupied by the x1v1 and x1v1 tori. This is a typical property of the 4D spaces of section of our system.
We have left for the end the description of the two tori that appear centred around in Fig. 14, because they bring in the system an important new family, which has to be discussed separately. Both tori belong to the same family, i.e. it is a p.o. of multiplicity 2. Just the size of the projected tori indicates the importance of this family. Its shape can be seen in the panels of the upper row of Fig. 16 (a,b,c). It is quite similar to the family x2mul2 found by Skokos et al. (2002a) to be bifurcated from x2 (panels Fig. 16d,e,f). Its view has a shape that contributes to the enhancement of a “CX” profile. In accordance with x2mul2 we call this family x1mul2. Both p.o. in Fig. 16 are given at E shortly after the bifurcation of x2mul2 from x2. We note that the vertical extent of x1mul2 is more than 60 times larger than that of x2mul2, hence it is a much more significant family for the overall morphology of the boxy bulges. The x1mul2 family is a vertical bifurcation of x1 at E . In order to describe the connection between x1 and x1mul2 we give in Fig. 17 the evolution of the stability indices of the x1 family when it is considered being of multiplicity 2 (red curves). The stability indices of x1mul2 emerge at E , where the index becomes tangent with the axis. We observe that x1mul2 is introduced in the system before the x1v1 and x1v2 families and extends to a large energy interval. In the interval E it is stable, while for E it becomes complex unstable. We note that the nonzero coordinates of the initial conditions of the p.o. of the x1mul2 family are and like the initial conditions of x1v2.
The structure of the 4D space of section in the neighbourhood of x1mul2 is as foreseen by Katsanikas et al. (2011b). Namely we have a couple of rotational tori with smooth colour variation along them. A typical example is given in Fig. 18. It depicts two 4D rotational tori around x1mul2 at E =. The chosen 3D projection is , while the consequents are coloured according to their values.
Sideon profiles of quasiperiodic and sticky orbits around x1mul2 support boxy structures. On top of this the orbits are selfintersecting close to in their sideon views. Most of the quasiperiodic orbits can be thought as thicker versions of the orbit depicted in the upper row of Fig. 16. Nevertheless, in other projections these orbits may have a complicated or chaotic looking morphology. As an example we give in Fig. 19a,b the quasiperiodic orbit that deviates from x1mul2 by , integrated for 12 periods of x1v2 at E = . In Fig. 19a the projection particular structure. However, in Fig. 19b, right panel, we have a conspicuous support of a boxy feature with an “CX” morphology. The faceon view , left panel of Fig. 19b, does not have a clear x1, character as in Fig. 16a, or Fig. 15b at the same energy. It has a morphology resembling quasiperiodic orbits around 3:1 stable orbits on the equatorial plane. Even so, it retains its narrowness in the direction, especially close to . Similar conclusions we draw by observing orbits in the immediate neighbourhood of x1mul2 close after the S transition. An example is given in Fig. 19c for the orbit that deviates by from the p.o at E . The orbit in Fig. 19c shares the same properties with the one in Fig. 19b in supporting a “CX’ sideon profile.
6.3 The phasespace close to complex unstable p.o. of multiplicity 2
In previous papers the structure of phase space has been studied in the neighbourhood of complex unstable simple periodic orbits, but not close to “” p.o. of multiplicity 2 or higher. Questions that arise concern the shape of the ’confined tori” (Pfenniger, 1985a, b) of multiplicity 2, their internal spiral structure (Contopoulos et al., 1994; Papadaki et al., 1995) and their 4D morphology (Katsanikas et al., 2011a). For this purpose we will briefly describe here the phase space structure in the neighbourhood of “” x1mul2 orbits, which are of multiplicity 2.
In the case of the simple periodic x1v1 family we know (Katsanikas et al., 2011a) that close beyond a S transition the 4D representations of the confined tori are almost 2D, slightly warped, disky objects with a smooth colour variation along their surfaces. Further away from the transition point the confined tori obtain a toroidal structure and finally, beyond a critical energy, we have scattered consequents in the 4D space. Katsanikas et al. (2011a) describe also phenomena of stickiness to confined tori.
This kind of behaviour is reproduced in the case of the x1mul2 family, however with some peculiarities due to the double multiplicity of the p.o. Let us consider an orbit at E that differs from the p.o. only by . We chose here such a small perturbation, because we want at this point to describe in detail the underlying dynamics. We use for our description the 3D projection and we colour the consequents according to their values. After 400 intersections two disky objects tend to form centred at ), as we can see in Fig. 20a. There is no colour mixing as we move across these objects. We have rather a smooth colour variation, meaning that all requirements exist for the formation of a confined torus in 4D. The way the consequents tend to fill these two disky structures follows a rule that can be better understood in the projection, i.e. if we observe the orbit of Fig. 20a from above. We give the first 100 consequents of the orbit in the projection in Fig. 20b. They form a set of 3 spirals with the p.o. in the centre of the configuration, like the spirals in the complex unstable p.o. of single multiplicity described in Katsanikas et al. (2011a). Three successive consequents belong to three different spiral arms. However, along each one of the spiral arms the successive consequents belong alternately to one of the two disky objects depicted in Fig. 20a. Red dots in Fig. 20b belong to the upper one in Fig. 20a , while black to the lower one . As we continue integrating the orbit, after the first 470 consequents, three tips are formed in each disky object. The subsequent consequents initially elongate these three tips and then they diffuse in the 3D phase space with their colours mixed (Fig. 20c). The diffusion starts before the two confined tori are filled by the consequents. The tips of the three extensions of each uncompleted confined torus are indicated with arrows in Fig. 20c. The whole areas of the confined tori as depicted in Fig. 20a are included now within the triangular light blue areas out of which the consequents diffuse in the 4D space. After 800 consequents we clearly see that we observe a unique cloud of points. The final stage of this evolution, reached after 2000 consequents is given in Fig. 20d. The situation remains like this, as we checked, even after 8000 consequents and is characterised by the formation of a volume around the p.o in which the distribution of consequents is sparse. The majority of the consequents is trapped in a cloud of points surrounding this rather “empty” volume (Fig. 20d). For even larger integration times we finally have a uniform distribution of consequents in all 3D projections. These results are of interest for the study of autonomous Hamiltonian systems in general. However, for the particular application we study here, we have to remember that we are interested for the orbital dynamics within tens and not thousands of dynamical times. Thus, we have to keep in mind that also close to multiplicity 2 complex unstable p.o. the spirals in phase space keep the initial tens of consequents in a small area of the phase space. Thus, they prevent the fast diffusion of orbits starting close to the p.o. and hence they reinforce, within this time, structures associated with them.
We note that for E =, besides x1mul2 the x1v1 family is also complex unstable. Perturbations of x1v1 by keep the consequents on a disky confined torus for thousands of consequences. So, during long times we do not have mixing in phase space due to the presence of two “” p.o.
6.4 A new sideon profile
Family x1mul2 is proven to be a very important family of p.o. for three reasons: (i) it exists over a large energy range having extended stable parts. (ii) it starts existing before the E of the vertical 2:1 resonance, i.e. it may trap material around it before x1v1. (iii) it builds a sticky zone connecting its two rotational tori, which is located around the x1v1 and x1v1 tori region (Fig. 14). Then on one hand we have orbits sticky to the x1mul2 tori and on the other hand the sticky zone blocks the easy communication between the x1v1  x1v1 tori region and the outer parts of the phase space.
Under these considerations our model indicates the possibility of b/p bulges consisting partly, or totally, by orbits associated with the x1mul2 family. A sideon profile composed by quasiperiodic and sticky orbits associated with the orbits of x1mul2 is given in Fig. 21a. We have taken five perturbed by x1mul2 orbits for E and . The x1mul2 p.o. at these energies are stable except for E , which is complex unstable. They form a clear butterfly structure with dimensions , i.e. it reaches a height . Including further orbits from the neighbourhood of x1mul2 orbits at higher energies results to a larger boxy bulge, however without retaining the conspicuous butterfly morphology. Combined with orbits building an Xshaped bulge, the “butterfly” will be located in the central area of the profile. How much pronounced it will appear in a model depends on the coexistence of other families at the same regions and on how much each of them will contribute to the building of the boxy bulge. For example the profile in Fig. 21b consists of the orbits close to x1v1 p.o., which have built the profile in Fig. 10c adding the orbits that give us the profile in Fig. 21a. A central butterfly feature is in this case discernible. In other models with with less extended “” parts in their stability diagrams the role of x1mul2 is more pronounced.
7 Geometrical considerations
The projection angle under which we view the b/p bulges is expected to play a role in their observed morphology. At any rate we study the morphology of disc galaxies observed nearly edgeon, thus practically we want to investigate just the changes introduced by the rotation of the galaxy around its rotational axis. By rotating the nonperiodic orbits that compose the profiles we have presented in Fig. 10, Fig. 12 and Fig. 21 we first realise that nonperiodic orbits in the neighbourhood of x1v2 and x1mul2 p.o. contribute to the formation of a b/p or “” feature only close to their sideon views. Away from this projection angle they only contribute to the formation of a boxy bulge, remaining confined in the direction.
As regards the orbits associated with the x1v1 family, besides the sideon view that gives the “OX” profiles (Fig. 10), there is another projection angle for which nonperiodic orbits in the neighbourhood of x1v1 p.o. could reinforce this time a “CX” type profile. For this, one has to consider both representatives of x1v1 and x1v1 at each energy. These projection angles are close to the endon view of the bar, i.e. we have for the box a ratio . In addition the resulting morphology is closer to an “”type or peanutshaped structure, than to an Xshaped feature with straight line segments as wings, as is for the sideon views. In Fig. 22a we give a set of x1v1 and x1v1 p.o. in the energy range E , from the point of view =(90,15) (for the endon view we have =(90,0)). Nonperiodic orbits close to the x1v1 and x1v1 p.o. however, do not build b/p profiles with features as sharp as in their sideon projections when combined and projected in other projection angles. The “”like structures they build are formed when viewed in a range of angles centred around an angle for which the x1v1 p.o. can reproduce the same feature much sharper. This angle is for quasiperiodic orbits after the complex unstable part of x1v1 (Fig. 1). For quasiperiodic orbits before the complex unstable part of the family (E ) however, the “” feature is best viewed with . In Fig. 22b we give the orbit “8” of Fig. 11 viewed with angles =(90,15) and in Fig. 22c the orbit “1” of Fig. 11 viewed with angles =(90,8). They are given together with their symmetric orbits with respect to the equatorial plane.
An angle is proposed by Gardner et al. (2014) for the viewing angle of the “Xshaped” bulge of the Milky Way. A b/p structure made by orbits associated with the x1v1 family projected at this angle is compatible with our findings.
The projections away from the sideon view of the nonperiodic orbits in Fig. 11, which are in the neighbourhood of “” x1v1 p.o. (“3” to “6”), give evidently less pronounced b/p features, although they have always a boxy character. In Fig. 23 we compare the morphology of the orbits “7” and “6” of Fig. 11 when viewed with angles =(90,13). We plot them always together with their symmetric with respect to the equatorial plane counterparts. The orbits in Fig. 23a with E are quasiperiodic, while those in Fig. 23b with E are chaotic. We realise the degree at which they enhance the “”type morphology close to their endon views. We note that if we consider all the orbits of the sample of Fig. 11 at a nearly endon view the orbits with the larger energies are the most important for characterising the overall morphology of the b/p bulge. Orbits at low energies remain in low heights above the plane and thus they do not shape the outer parts of the bulge.
We note that one can in many cases find combinations of projection angles for individual nonperiodic orbits in the neighbourhood of x1v1 that projected resemble structures with a “CX” type morphology. As an example we give the perturbed by orbit corresponding to the fourth invariant curve around x1 in Fig. 4, starting with , from a point of view defined by the angles (60,82) (Fig. 24). It is integrated for about 50 dynamical times of x1 at the same energy (E ). However, away from the edgeon orientation (in our example in Fig. 24 ) the combination of many orbits projected with the same angles in the cases we tried lead in the best case to a boxy, but not to a b/p structure.
8 Discussion and Conclusions
The main result of the present study is that there is a significant contribution from sticky chaotic orbits to the building of the body of b/p structures in the sideon profiles of disc galaxy models. The families of the x1tree can be considered as the backbone of the peanut structure in two ways. On one hand they trap around them quasiperiodic orbits and on the other hand they create obstacles for chaotic orbits that hinder their diffusion in phase space. These obstacles are the tori around orbits belonging to families of the x1tree (mainly of x1v1 and x1v1). The tori act as attractors for chaotic orbits, which become sticky to them. The important region for the peanut structure is the E interval of the radialILR and verticalILR. We find that in the frame of the current paradigm a new family plays a key role, namely x1mul2, which is of multiplicity two. This family of p.o. can either reinforce b/p sideon profiles based mainly on a x1v1x1v1 backbone, or even support a b/p morphology just by itself. In more detail our conclusions can be summarised as follows:

We find nonperiodic orbits supporting a b/p profile in the sideon view before the critical energy in the vILR region at which x1v1 is bifurcated from x1. These orbits are quasiperiodic orbits on tori around stable x1 orbits. Their outline tends to a combined x1v1x1v1 morphology as E increases towards the value at which x1v1 is introduced in the system (Fig. 2). This tells us firstly that the peanut starts before the vILR resonance. In our model, peanutsupporting orbits are found extending along the major axis of the bar already at y=0.13 bar lengths. Secondly, that the structure of the phase space in the neighbourhood of the mother family changes in a way that the coming of the bifurcated family can be foreseen by the gradual change in the shape of its quasiperiodic orbits. We do not have an abrupt change of the phase space structure immediately upon the introduction of the bifurcated family in the system.

Nonperiodic orbits in the neighbourhood of x1 reinforce boxy profiles even in the interval in which x1 is simple unstable. These are orbits, which either become sticky to the nearby x1v1 and x1v1 rotational tori, or in general orbits that can be observed being trapped in cylindrical structures in 3D projections of the 4D spaces of section that include the plane. This is a significant dynamical mechanism for building the body of the peanut.

The neighbourhood of x1 p.o. is a substantial source of nonperiodic orbits that contribute to boxy sideon profiles. They provide such orbits to the system before the vILR and in the E region where it is U. Of particular significance is its contribution in regions where x1 is the only stable family of p.o. in the system (in our model E ). These orbits contribute to a b/p morphology mainly by becoming sticky to x1v1 and x1v1 rotational tori. Else, they contribute to an overall edgeon boxiness.

In the E region, apart from the U x1 p.o. we have also the U x3 p.o. The structure of the phase space around them is different, since they are unstable to vertical and radial perturbations respectively. The structure of the phase space in the neighbourhood of U p.o. is determined by the type of simple instability of the p.o.

Edgeon profiles of nonperiodic orbits in the neighbourhood of x2 and x3 have a certain similarity (Fig. 8). Thus they will have similar contribution to b/p profiles of galactic bars. They reach heights about 100 pc, considerably smaller than the heights of the boxy bulges in images of edgeon galaxies. This indicates that if x2 stellar flows exist in 3D stellar bars, they are embedded inside the edgeon profiles of the b/p structures.

In the case of Xshaped features built by orbits associated with the x1v1 family, the “X” is formed by the tips of orbits with rather “”like or “”like morphologies in their sideon views. Along the wings of “X” we have (Fig. 11). The morphology of a b/p structure is always characterised by the outermost orbits. Models in which we observe sharp wings emerging out of boxy structures can be constructed by populating orbits in the x1v1 neighbourhood at E beyond the part of the x1v1 family. The part of the b/p bulge supported by chaotic orbits at energies where x1v1 is , is expected to have less conspicuous “X  wings”. Bars with moderate amplitudes can support more efficiently a Xshaped b/p bulge based on x1v1 dynamics, because in this case this family is complex unstable over smaller E ranges.

“CXtype” profiles can be supported:

By quasiperiodic orbits around x1 in energies, where the x1v2 family coexists. They become more important as we increase their initial condition from 0 towards the value of x1v2. The strongest reinforcement is given by orbits corresponding to the closest to x1v2 tori around x1.

By sticky chaotic orbits around the x1v1 and x1v1 tori. We find such orbits by perturbing the x1v2 initial conditions, since the unstable asymptotic curves around x1v2 wind around the tori of the x1v1 and x1v1 families.


A very important role for the dynamics of the b/p bulge is played by the family of multiplicity 2, x1mul2. This family exists already at energies smaller than the one at which x1v1 is introduced in the system, being stable over large E energy intervals. Quasiperiodic orbits trapped around the orbits of this family and, mainly, sticky orbits to them reinforce a b/p structure.

We find chaotic orbits in the neighbourhood of complex unstable periodic orbits that can contribute to the reinforcement of a b/p structure for several tens of dynamical times. In our model, in the case of chaotic orbits close to complex unstable x1v1 p.o. there is a correlation between the degree of confinement of the chaotic orbit in the configuration space and the value of the quantity (Eq. 9). As we have seen in Fig. 11, the smaller the quantity is, the largest the deviation of its morphology from that of a quasiperiodic orbit. Interesting for the general theory of Hamiltonian systems is the correspondence found between the structure of the phase space in the neighbourhood of simple and double complex unstable p.o. In the case of p.o. of double multiplicity the 4D spiral structure is developed in two confined tori simultaneously. We see also that tiny perturbations of complex unstable periodic orbits do not lead to an abrupt filling of the available phase space by the consequents of the two orbits. In a case with two “” p.o. (one of simple and one of double multiplicity), which coexist at the same E , we have found orbits in their immediate neighbourhood that stay in different regions of the phase space for thousands of consequents.

Considering nonperiodic orbits at successive energies introduced in the system through one of the dynamical mechanisms described in our paper, one can build conspicuous b/p profiles mainly in their sideon views. Away from the sideon view we have in general just a boxy morphology in the outer envelopes of the edgeon projections. However, in the case of Xshaped boxy bulges based on x1v1 dynamics, we find nearly endon views, where the same orbits can give an “”type morphology.
Acknowledgements
We thank Prof. G. Contopoulos for fruitful discussions and very useful comments. This work has been partially supported by the Research Committee of the Academy of Athens through the project 200/815.
Footnotes
 This ratio is usually larger in body simulations than the corresponding quantity estimated for real galaxies (Erwin, 2005)  see e.g. models that develop peanuts in Athanassoula & Misiriotis (2002) and in Saha & Naab (2013)  and may change during the life of a galaxy (Cheung et al., 2013)
 Tiny perturbations close to the bifurcating point, at which x1v2 is introduced in the system, secure almost regular behaviour for several Hubble times (Katsanikas et al., 2013). Nonetheless, this is not of practical use because of the needed proximity of the initial conditions of the perturbed orbits to the periodic orbit.
References
 Aronica G., 2006, PhD Thesis, RuhrUniversität Bochum
 Athanassoula E., Misiriotis A., 2002, MNRAS 330, 35
 Broucke R., 1969, “Periodic orbits in the elliptic restricted threebody problem”, NASA Techn. Rep. 32, 1360
 Bountis T., Manos T., Antonopoulos C., 2012, Cel.Mech.Dyn.Ast. 113, 63
 Bureau M., Athanassoula E., 2005, ApJ 626, 159
 Bureau M., Aronica G., Athanassoula E., Dettmar R.J., Bosma A., Freeman, K. C., 2006, MNRAS, 370, 753
 Cheung E., Athanassoula E., Masters K.L. et al. 2013, ApJ 779, 162
 Combes F., Debbasch F., Friedli D., Pfenniger D., 1990, A&A 233, 95
 Contopoulos G., 1981, A&A 102, 265
 Contopoulos G., 2004, “Order and Chaos in Dynamical Astronomy”, Springer Verlag, BerlinHeidelberg
 Contopoulos G., Barbanis B., 1985, A&A 153, 44
 Contopoulos G., Grosbøl P., 1989, A&AR, 1,261
 Contopoulos G., Harsoula M., 2010, Cel.Mech.Dyn.Ast. 107, 77
 Contopoulos G., Harsoula M., 2013, MNRAS 436, 1201
 Contopoulos G., Magnenat P., 1985, Celest. Mech. 37, 387
 Contopoulos G., Farantos S.C., Papadaki H., Polymilis C., 1994, Cel. Mech. 37, 387
 Erwin P., 2005, MNRAS 364,283
 Gardner E., Debattista V., Robin A., et al. 2014, MNRAS 438, 3275
 Hadjidemetriou J., 1975, Celest. Mech. 12, 255
 Heggie D.C., 1985, Celest. Mech. 35, 357
 Hickson P., 1993, Astr. Lett. Comm. 29, 1
 Jorba A., Olle M., 2004, Nonlinearity 17, 691
 Katsanikas M., Patsis P.A., 2011, Int. J. Bif. Ch. 2102, 467
 Katsanikas M., Patsis P.A., Contopoulos G., 2011a, Int. J. Bif. Ch. 2108, 2321
 Katsanikas M., Patsis P.A., Pinotsis A.D., 2011b, Int. J. Bif. Ch. 2108, 2331
 Katsanikas M., Patsis P.A., Contopoulos G., 2013, Int. J. Bif. Ch. 2302, 1330005
 Kregel M., van der Kruit P.C., Freeman K.C., 2004, MNRAS 351, 1247
 Lange S., Richter M., Onken F. et al., 2014, arXiv:1311.7632 [nlin.CD]
 Li ZY., Shen J. 2012, ApJL 757L, L7
 Lütticke R., Dettmar R.J., Pohlen M., 2000, A&A 362, 435
 MartinezValpuesta I., Shlosman I., Heller C., 2006, ApJ 637, 214
 Miyamoto M., Nagai R., 1975, PASJ, 27,533.
 Olle M., Pacha J.R., Villanueva J., 2004, Cel. Mech. Dyn. Astr. 90, 89
 Papadaki H., Contopoulos G., Polymilis C., 1995, in “From Newton to Chaos’’, Roy A.E., Steves B.A. (eds), pp 485494, Plenum Press, New York
 Patsis P.A., 2005, MNRAS 358, 305
 Patsis P.A., Grosbøl P., 1996, A&A 315, 371
 Patsis P.A., Katsanikas M., 2014, MNRAS (in press).
 Patsis P.A., Xilouris E., 2006, MNRAS 366, 1121
 Patsis P.A., Zachilas L., 1990, A&A 227, 37
 Patsis P.A., Zachilas L., 1994, Int. J. Bif. Ch. 4, 1399
 Patsis P.A., Skokos Ch., Athanassoula E., 2002, MNRAS 337, 578
 Pfenniger D., 1984, A&A 134, 373
 Pfenniger D., 1985a, A&A 150, 97
 Pfenniger D., 1985b, A&A 150, 112
 Plummer H.C., 1911, MNRAS 71, 460
 Quillen A.C., Minchev I., Sharma S. et al., 2014, MNRAS 437, 1284
 Richter M., Lange S., Bäcker A., Ketzmerick R., 2014, PhRvE, 89,022902
 Saha K., Gerhard O., 2013, MNRAS 430, 2039
 Saha K., Naab T., 2013, MNRAS 434, 1287
 Saito R.K., Zoccali M., Mc William A. et al., 2011, AJ 142, 76
 Skokos Ch., 2001, Physica D 159, 155
 Skokos Ch., Patsis P.A., Athanassoula E., 2002a, MNRAS 333, 847
 Skokos Ch., Patsis P.A., Athanassoula E., 2002b, MNRAS 333, 861
 Vásquez S., Zoccali M., Hill V. et al. 2014, A&A 555, 91
 Vrahatis M., Isliker H., Bountis T.C., 1997, Int. J. Bif. Ch. 7, 2707
 Wakamatsu K., Hamabe M., 1984, ApJS 56, 283
 Wegg C., Gerhard O., 2013, MNRAS 435, 1874
 Williams M. J, Zamojski M. A., Bureau M. et al., 2011, MNRAS 414, 2163
 Zachilas L., A&AS 97, 549
 Zachilas L., Katsanikas M., Patsis P. A. 2013, Int. J. Bif. Ch. 23, 1330023