The phase-space of boxy-peanut and X-shaped bulges in galaxies I. Properties of non-periodic orbits

The investigation of the phase-space 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 boxy-peanut (b/p) and"{\sf X}-shaped"structures in rotating potentials of galactic type. We underline the significance of the rotational tori around the 3D families x1v1 and x1v1$^{\prime}$ that have been bifurcated from the planar x1 family. These tori play a multiple role: (i) They belong to quasi-periodic 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$^{\prime}$). Their unstable asymptotic curves wind around the x1v1 and x1v1$^{\prime}$ 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. Our work shows also that there are peanut-supporting orbits before the vertical ILR. Non-periodic orbits associated with the x1 family secure this contribution as well as the support of b/p structures at several other energy intervals. Non-linear 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$^{\prime}$ periodic orbits are pronounced only in side-on or nearly end-on views of the bar.


INTRODUCTION
The skeleton of the boxy and peanut-like (b/p) structures observed in edge-on disc galaxies is considered to be the families of periodic orbits (hereafter p.o.), which constitute the x1-tree. These are the families bifurcating from the planar x1 family (Contopoulos 1981) which lies on the equatorial plane, at the n : 1 vertical resonances (Skokos et al. 2002a,b;Patsis et al. 2002). Despite the fact that the periodic orbits determine their environment in the phase-space, the main body of the structures are non-periodic orbits. It is known that even in 2D Hamiltonian systems quasi-periodic 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 quasi-periodic orbits on invariant curves close to it. Contrarily, quasi-periodic orbits on invariant curves at ⋆ patsis@academyofathens.gr (PAP) † mkatsan@academyofathens.gr (MK) the outer parts of stability islands may support more complicated structures. For example quasi-periodic 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. 2011aKatsanikas et al. ,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 colour-rotation 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 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 colour-rotation technique in order to investigate the phase space of the b/p and X-shaped structures in our model. This method has been proven to be especially efficient in estimating the role of sticky orbits .
Referring to b/p-shaped 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 edge-on 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 off-centred 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 R max and a maximum height vertical to the equatorial plane z max . 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 edge-on 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 2.7±0.3. Aronica (2006) (his appendix C) presented the unsharp-masked profiles of a sample of 25 edge-on 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 unsharp-masked images of this sample, one can empirically estimate roughly a variation of the r b/p = R max /z max ratio, to be 0.8 r b/p 1.5. We estimated a mode value close to r b/p = 1.4. 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 end-on and large for bars viewed nearly side-on, since z max 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). Non-periodic 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)  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 edge-on 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 side-on 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.

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: where M D is the total mass of the disc, A and B are the horizontal and vertical scale lengths, and G is the gravitational constant. For the bulge we used a Plummer sphere (Plummer 1911) with potential: where ǫ s is the scale length of the bulge and M S is its total mass. Finally, the density ρ(x) of a triaxial Ferrers bar is: where a, b, c are the semi-axes and M B is the mass of the bar component. The corresponding potential Φ B 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 a : b : c = 6 : 1.5 : 0.6, as in Pfenniger (1984) and Skokos et al. (2002b). The masses of the three components satisfy G(M D + M S + M B ) = 1. The length unit is taken as 1 kpc, the time unit as 1 Myr and the mass unit as 2 × 10 11 M ⊙ . 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 The 3D bar is rotating around its short z axis. The x axis is the intermediate and the y axis the long one. The system is rotating with an angular speed Ω b and the Hamiltonian governing the motion of a test-particle can be written in the form: where p x , p y , and p z are the canonically conjugate momenta. We will hereafter denote the numerical value of the Hamiltonian by E J and refer to it as the Jacobi constant or, more loosely, as the 'energy'. The corresponding equations of motion are: In Table 1 we summarise the parameters of the specific model we will use in our paper.

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 (x 0 , p x 0 , z 0 , p z 0 ) in the plane y=0, for p y > 0. 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 −11 . For the integration of the orbits we used a fourth order Runge-Kutta scheme. The relative error in the energy remains in our calculations less than 10 −15 .
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 T : R 4 → R 4 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: ξ = M ξ 0 . Here ξ is the final deviation, ξ 0 is the initial deviation and M is a 4× 4 matrix, called the monodromy matrix. It can be shown that the characteristic equation is written in the form λ 4 +αλ 3 +βλ 2 +αλ+1 = 0. Its solutions (λ i , i = 1, 2, 3, 4) obey the relations λ 1 λ 2 = 1 and λ 3 λ 4 = 1 and for each pair we can write: where b i = 1/2 (α ± ∆ 1/2 ) and The quantities b 1 and b 2 are called the stability indices. One of them is associated with radial and the other one with vertical perturbations. If ∆ > 0, |b 1 | < 2 and |b 2 | < 2, the 4 eigenvalues are on the unit circle and the periodic orbit is called 'stable'. If ∆ > 0, and |b 1 | > 2, |b 2 | < 2, or |b 2 | > 2, |b 1 | < 2, two eigenvalues are on the real axis and two on the unit circle, and the periodic orbit is called 'simple unstable'. If ∆ > 0, |b 1 | > 2, and |b 2 | > 2, all four eigenvalues are on the real axis, and the periodic orbit is called 'double unstable'. Finally, ∆ < 0 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 stable, simple unstable, double unstable and complex unstable 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 J ), while all other parameters remain constant (Contopoulos & Barbanis 1985;Contopoulos & Magnenat 1985). In other words it gives the variation of the stability indices b 1 and b 2 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 b = −2 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 b = 2 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 "x1-tree" is composed by the x1 family itself together with the 3D families, which are generated at specific S→U or U→S transitions. At these transitions the stability index b 2 associated with the vertical perturbations has two intersections, or a tangency, with the b = −2 axis. They happen at the regions of the vertical resonances (see e.g Patsis & Grosbøl 1996). The "x1-tree" is the backbone of b/p bulges and X features (Patsis et al. 2002). In standard models the first S→U and U→S transitions due to an intersection of b 2 with the b = −2 axis are encountered near the vertical 2:1 resonance and the bifurcated simple periodic 3D families of p.o. of the x1-tree 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 x1-tree as E J varies is given in Fig. 1a. The two stability curves, that exist already for E J ≦ −0.45 (left side of the figure), belong to x1. The stability index b 1 is associated with the radial, while b 2 with the vertical perturbation. The families x1v1, x1v2 etc. are introduced at the intersections of b 2 with the b = −2 axis. In Fig. 1b we zoom into the ILR region, which is of special importance for the b/p dynamics. We will refer Table 1. The parameters used in our model. G is the gravitational constant, M D , M B , M S are the masses of the disk, the bar and the bulge respectively, ǫ s is the scale length of the bulge, Ø b is the pattern speed of the bar, E J (r-ILR) and E J (v-ILR) are the values of the Jacobi constant for the radial and vertical 2:1 resonances, given with higher accuracy than in Skokos et al. (2002b), and R c is the corotation radius. In the grey shaded region ∆E vILR (see text) x1 is vertically unstable. (c) The four eigenvalues of the family x1v1 (λ i , 1/λ i , i = 1, 2) in its complex unstable part.
to it in detail in section 4. The straight line segments in Fig. 1a and b at b 1,2 = 0 for E J −0.3965, labelled "∆(x1v1)" is drawn to mark the area over which the x1v1 family is complex unstable. In this region the inidices b 1 and b 2 in Eq. 8 cannot be defined since we have ∆ < 0 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.

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 J (r-ILR) in Table 1 is taken as the left limit of the x2-x3 loop in the characteristic or in the stability diagram. E J (v-ILR) is taken at the energy where we have the S→U 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 −0.438317 <E J < −0.424945 (Fig. 1a). One more 3D family (x1v2) associated with the vertical ILR starts at the right border of this interval where we have the U→S transition (Fig. 1a).
It is generally believed that the peanut starts at the vertical 2:1 (v-ILR) resonance, where we encounter 3D p.o. with elliptical projections on the equatorial plane and "smiles" or "frowns" in their side-on 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;Martinez-Valpuesta et al. 2006). We note that Combes et al. (1990) emphasise that at the end of their N-body models the radial and vertical 2:1 resonances coincide.
The bar is elongated along the y-axis and its backbone is the x1-tree. 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 non-periodic orbits that support the disc in the same region. Among them there are 3D quasi-periodic orbits with E J < −0.438317 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 side-on profiles. For E J < −0.438317, the only existing bar-supporting family of p.o. belonging to the x1-tree 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, quasi-periodic 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 z or p z 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 S→U transition of x1 at the vertical 2:1 resonance (E J ≈ −0.438317).
In Fig. 2a we depict the (y, z) projections of two quasi periodic orbits around x1 at E J = −0.527983. Both of them have the same x 0 initial condition as the x1 p.o. at this energy (x 0 = 0.072730404). The narrow, darker one, is perturbed in the z direction by ∆z = 0.01, while the thicker one by ∆z = 0.1. 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 z has a shorter extent along the y direction. In the 3D configuration space these orbits have a toroidal morphology. In the (y,z) projections they are characterised by a minimum z for y = 0 and a maximum z at maximum |y| as in Fig. 2a. This is a typical, general morphology of the quasi-periodic orbits around x1 away from the ILRs of the system. However, as energy increases we observe that initially for small perturbations, the minimum z for y = 0 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 J = −0.45 having initial conditions that deviate from those of the x1 p.o by ∆z=0.01. In order to clearly observe the formation of the local maximum, the scale of the z-axis is 0.02 times that of the y-axis.
By studying the quasi-periodic orbits around x1 for energies just before the S→U transition of x1 at the vertical 2:1 resonance, we realise that the local z maximum at y = 0 becomes prominent for quasi-periodic orbits in the neighbourhood of x1. The orbit in Fig. 2c is for E J =−0.438666 and occurs by perturbing the x1 p.o. with initial conditions (x 0 , z 0 , p x , p z )=(0.0930694109,0,0,0) by ∆z=0.01. The projection of the orbit has taken a shape like an " " symbol. The p.o. x1 at this energy is still stable (b 1 ≈ −1.98, b 2 ≈ 1.98). Despite the fact that we are before the S→U 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 quasi-periodic 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 J = −0.438225) when they have been introduced in the system as stable p.o. a f ter the S→U 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 colour-rotation method (Patsis & Zachilas 1994), these tori appear in the case of family x1, before the v-ILR, as typical "tube tori". This happens because there is a 3D projection, (x, z, p x ), in which the torus appears as folded and self-intersecting (Vrahatis et al. 1997). In Fig. 3a we depict the perturbed x1 orbit by ∆p z = 0.01 at E J = −0.438666. For constructing this 4D object the orbit is integrated for 500 periods of the x1 p.o. at the same energy. We use the (x, z, p x ) projection and the consequents are coloured according to their p z 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 (θ, φ) = (270 • , 0 • ) 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 pseudo-intersection. The figure gives a typical 4D Poincaré surface of section for all 3D bar-supporting quasi-periodic orbits before the first S→U 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 (θ, φ) = (60 • , 266 • ). The quasi-periodic orbit is drawn with black. The red curve is the x1 p.o. on the z = 0 plane at the same energy. One can easily understand that by viewing this orbit side-on, the shape will be as the one depicted in Fig. 2c.
The study of successive (y, z) projections shows that the shape  of the quasi-periodic 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 y = 0 in the (y, z) projections, we can say that we have peanut-supporting orbits already at a distance of about 0.6 kpc along the major axis of the bar.

Non-periodic orbits in the neighbourhood of x1
Thus far we explored the contribution of the quasi-periodic 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 S→U transition of x1, at E J ≈ −0.438317 marks the beginning of the vertical ILR region, while in our model the radial ILR already exists at this energy. Beyond E J ≈ −0.438317 the x1 family becomes simple unstable, with b 2 < −2, for −0.438317 < E J < −0.424945. Let us call this energy interval ∆E vILR . At the end of the ∆E vILR interval we have the U→S transition of x1. Within ∆E vILR 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 vILR 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 b 1 index of x1v1 appears coinciding with the b = 2 axis. In reality b 1 just remains always very close to the b = 2 axis, having two tangencies at E J ≈ −0.4104 and at E J ≈ −0.3982. At E J ≈ −0.396 b 1 meets b 2 and the x1v1 family becomes complex unstable. The stability indices associated with the x2-x3 loop are plotted with red.
In ∆E vILR x1 is unstable relative to vertical perturbations, because b 2 < −2. However, for perturbations on the equatorial plane x1 remains stable, because −2 < b 1 < 2. A (x, p x ) cross section in ∆E vILR with orbits deviating from x1 by ∆x or ∆p x , and thus remaining on the equatorial plane, will be as the one in the example depicted in Fig. 4 for E J = −0.438225. We note that despite the fact that for this E J 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 J = −0.438225 (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 J = −0.438317. 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 z and in p z directions, for which the integrated orbits will be orbits on the z Figure 5. Perturbed x1 orbits at E J =-0.438225. (a) The red curves numbered "1","2", "3" and "4", correspond to x1 orbits perturbed by ∆z = ±0.0001, ±0.005, ±0.01, ±0.015. The blue curve is for the x1 orbit perturbed by ∆z = 0.02. The locations of x1, x1v1 and x1v1 ′ are indicated with arrows. The dots in the centres of the "stability islands" correspond to tiny tori around x1v1 and x1v1 ′ . (b) The same orbits as in (a) in the 3D (z, p z , x) projection. (c) The (y, z) projection of the orbit numbered "3" in (a). (d) The (y, z) projection of the orbit on the torus surrounding all three p.o. in (a) and (b). It has a peanut-like shape.
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 x 0 as x1 and a small z 0 0. In Fig. 5a we give the (z, p z ) projection of the consequents of the x1 orbits perturbed successively by ∆z = ±0.0001, ±0.005, ±0.01, ±0.015 (red curves) and 0.02 (blue curve), at E J =-0.438225. The location of x1 is at (z, p z ) = (0, 0). The orbits on the red curves, for ∆z = ±0.0001 to 0.015, 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 0.02 surrounds in the (z, p z ) plane all three p.o. In the 3D (z, p z , x) projection we observe (Fig. 5b) that all these orbits are on tori, the thickness of which (x-dimension) increases with increasing ∆z. The side-on morphology of the orbits on the red curves is as the one in Fig. 5c (z 0 = 0.01) 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 z 0 = 0.02, have a clear peanut structure, similar to the one we have seen in Fig. 2c, at an energy before the S→U transition of x1. This is given in Fig. 5d.
If we include in the vertical (by ∆z, ∆p z ) perturbations of x1 also ∆x or ∆p x displacements, i.e. if we essentially perturb vertically quasi-periodic 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 ∆z 0.1 starting in this region belong to x1v1 invariant tori. Differentiations from this behaviour is observed when we start integrating orbits with z or p z 0 starting at and beyond the fourth closed curve around x1 in the x1 "stability island" of Fig. 4. There, for the same range of ∆z or ∆p z 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 phase-space 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 ∆z = 0.01 (a) and by ∆z = 0.1 (b). The side-on 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 side-on 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 quasi-periodic orbit morphology similar to the one in Fig. 2a even if we continue its integration for 10 5 dynamical times. In Patsis & Katsanikas (2014) (paper II), we will return to these orbits discussing orbits that combine boxy side-on 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 side-on 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 (x, p x ) 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 side-on morphologies with the orbit in Fig. 6b for thousands of dynamical times. Their sticky character is revealed very illustrative by means of the colour-rotation method. The 3D projections of the space of section show us that these orbits remain confined for large times (of the order of 10 2 to 10 4 dynamical times) in certain regions of the phase space. We remind that usually in Galactic Dynamics we are interested about time-scales of tens or in extreme cases a few hundreds of dynamical times.
As an example, let us consider the orbit with initial conditions (x 0 , z 0 , p x 0 , p z 0 ) = (0.304, 0.1, 0, 0) and integrate it for 10 4 dynamical times. In Fig. 7a we plot this orbit with blue points in the (x, p x ) 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 (x, p x , z) projection (Fig. 7b) we realise that the two orbits ("red" and "blue") build in fact cylindrical structures that remain confined also in the z-direction. 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 colour-rotation method. They are given in Fig. 7c and Fig. 7d for the "red" and "blue" orbits respectively. In Fig. 7c we use (x, z, p x ) as spatial coordinates and we colour the consequents according to their p z values for the "red" orbit. In Fig. 7d we use (x, z, p z ) as spatial coordinates and we colour the consequents according to their p x 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 (x, p x ) 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 side-on profiles. In 3D projections of the 4D spaces of section, which include the (x, p x ) 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.

The x2, x3 region
Apart from x1, the other simple unstable 2D family in the ∆E vILR region is x3. Unlike x1 it is radially unstable. On the (x, p x ) 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 ∆z the initial conditions of x3 we find in its neighbourhood orbits confined in the z-direction of the configurations space. However, they are not as sharply defined as the quasi-periodic orbit around x2. For comparison with x2, we give the morphology of two perturbed by ∆z = 0.1 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 x2-flow (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 x-axis (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 non-periodic orbits around x3 to the local density away from the equatorial plane is comparable with that of quasi-periodic orbits around x2. However, they have a clearly different morphology in their face-on 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 colour-rotation method to visualise the 4D space. In this particular case we use the (x, z, p z ) projection and we colour the consequents according to their p x variation, normalised in the [0,1] interval. The x1 orbit is perturbed by ∆z = 0.1. 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 Figure 7. A chaotic orbit at E J = −0.438225 trapped for thousands of dynamical times of x1 in the region between the last drawn invariant curve around x1 and a chain of eight stability islands in the (x, p x ) projection of the 4D space of section (Fig. 4). It is given with blue coloured consequents in (a) and (b). In (c) and (d) we give 4D representations for the "red" and "blue" orbits of the panels (a) and (b) respectively.  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 (x, z, p z ) 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 p x 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. co-existing 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 vILR give an illustrative example of how the dynamical phenomenon of stickiness becomes essential for building the peanut as E J increases approaching corotation. Similar phenomena take place at larger E J 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 S → U transition is given by Katsanikas & Patsis (2011) and by Katsanikas et al. (2013) in a double Miyamoto-Nagai 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.

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 X-shaped bulges with dimensions similar to those observed in edge-on 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 x1-tree p.o on the major axis of the bar (y-axis) reaches a radius y ≈ 4.6. 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 1 .
The maximum length of the b/p bulge in the side-on projection is expected to be according to Lütticke et al. (2000) R max ≈ 1.7 (see Sect. 1). Then, by measuring the dimensions of the X features in Aronica (2006)  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 ∆z. The perturbations were such as to secure a clearly by eye discernible " " feature formed when we plot the side-on views of the non-periodic orbits and their symmetric with respect to the equatorial plane of the model. This ∆z perturbation is less than 0.1 pc and is added to the z 0 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 ∆x. This time we obtain " " morphologies for ∆x perturbations of the individual orbits in the range 0.05 < ∆x < 0.25. In Figs. 10c,d we start our orbits with p z =0.05 and p z =0.1 respectively, while the initial p z 0 component for orbits of the x1v1 family is 0. In all cases the individual orbits build in their (y, z) (side-on) 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 z-maxima of the orbits, where p z = 0 (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 ∆x or by ∆p z , with initial z 0 = 0 form a sharper X profile than the orbits starting away from the equatorial plane z 0 0. We underline the fact that by constructing these profiles we have applied the perturbations to the x1v1 p.o. independently of their stability.

Contribution by orbits close to complex instability
Qualitatively, the orbits that participate in the building of "x1v1based" X features reaching heights up to z = 1.2 are of two kinds. They are either quasi-periodic orbits trapped by stable x1v1 p.o. or chaotic orbits associated with complex instability. The former are encountered at energies −0.438 < E J < −0.396 and E J > −0.301, while the latter at −0.396 < E J < −0.301. The quasi-periodic 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 J 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 ∆ < 0. 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 J varies. From "1" to "8" we mark in this diagram with heavy dots the locations of the orbits we used. The (y, z) 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 |z| 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 J = −0.37 is a particular case, since it is more affected by the presence of the family x1v2, which for E J = −0.37 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 ∆ < 0. The same happens with orbits perturbed by ∆z and ∆x 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 colour-rotation 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 10-20 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). y z Figure 11. In the central frame we give the variation of ∆ with the energy, for 8 x1v1 p.o. The embedded frame depicts the side-on projections of these p.o. With black are the stable and with orange colour the complex unstable ones. The dashed arrows indicate that the wings of the X are formed along the z maxima of the orbits. The panels surrounding the central frame give the side-on views of the non-periodic orbits used in the profile of Fig. 10c. The numbers in the lower right corner of these panels corresponds to the energies numbered along the ∆ curve. We observe that smaller ∆ values correlate with more diffused orbits in phase space at the corresponding energy. Orbit "3" is a particular case (see text). 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.

"OX" VS "CX" PROFILES
The side-on 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 high-resolution image of this object, see the Space Telescope Science Institute's news release STScI-2009-30). This is the kind of X features labelled by Aronica (2006) and Bureau et al. (2006) as "OX" (off-centred). Indeed, in many images of galaxies the direction of the wings of the X, as revealed after applying an unsharp-masking 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 edge-on galaxy HCG87A in the Hickson 87 group (Hickson 1993) (for a high resolution image, see the STScI-1999-31a news release). Aronica (2006) and Bureau et al. (2006) call this type of X, "CX" (centred). The side-on "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 "x1v1-X" closer and form "CX" instead of 'OX" profiles will be discussed in section 7. Be-low we examine the possibility of building true "CX" profiles using orbits.

Possible x1v2 contribution
In nearly side-on 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 x1-tree, 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 2.7 ± 0.3, as we mentioned in the introduction.
Under those considerations, a family candidate to play this role would be x1v2, bifurcated from x1 at a U→S transition, in our model at E J ≈ −0.425 (see Fig. 1). The side-on projections of the p.o. of this family have an ∞-type morphology (Skokos et al. 2002a). In Fig 12a we plot a typical side-on view of a x1v2 p.o. orbit (red) together with the side-on 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 (y, z) plane. However, the x1v2 family is initially simple unstable, inheriting the simple unstable character of x1, while for larger E J it evolves to double unstable (Fig. 1a). Thus, the mechanism of trapping quasi-periodic orbits around x1v2 p.o. cannot be applied.
Nevertheless we may look for non-periodic 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 −0.425 < E J < −0.35, 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 .
We know ) that in a range of E J away from the immediate neighbourhood of the bifurcating point 2 the unstable asymptotic curves of the (U) x1v2 p.o. wind around the rotational tori of the stable x1v1 and x1v1 ′ families that exist at the same energy. The construction of the asymptotic curves is extensively described in Katsanikas et al. (2013). Here, we will examine how this applies in the case of our 3D bar model.
For energies beyond the bifurcation point, we have integrated x1v2 orbits perturbed by ∆x, ∆p x , ∆z, or ∆p z for equal time as for the orbits composing the profiles of Fig. 10. We first realised that such orbits, with E J < −0.35, 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 ∆x. Perturbed x1v2 orbits by ∆x 0.1 are able to reinforce to some degree an ∞-type morphology in their side-on projections. For example, by substituting the orbits used for E J = −0.41, −0.39 and −0.37 in Fig. 10b by x1v2 orbits perturbed by ∆x = 0.1, 0.1 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.
2 Tiny perturbations close to the bifurcating point, at which x1v2 is introduced in the system, secure almost regular behaviour for several Hubble times . Nonetheless, this is not of practical use because of the needed proximity of the initial conditions of the perturbed orbits to the periodic orbit. By varying E J , we realise that the ability of perturbed x1v2 orbits to support an ∞-type side-on 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 U→S transition of x1 (E J ≈ -0.425) the x1v1 tori are still in the neighbourhood of the p.o. x1v2. At E J = -0.42398, the stability index of x1v2, which is associated with vertical perturbations, is b 2 = −2.2. Its initial conditions are (x 0 , z 0 , p x 0 , p z 0 ) ≈ (0.10476, 0, 0, 0.055), while the corresponding initial conditions for x1v1 at the same energy are (0.10376, 0.159, 0, 0). 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 |z| 0.1, as is the orbit in Fig. 13a. We find that about the first 20 consequents of such orbits in the (z, p z ) 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 (z, p z ) 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 ). However, we find anew an enhancement of x1v2 type morphologies at energies −0.39 E J −0.37. 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 y-axis family, o2 (Skokos et al. 2002b). In Fig. 1 we can see how o1 is introduced in the system at E J ≈ −0.39, when the b 1 stability index of x1 becomes b 1 < −2. 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 side-on views of the models. We find them being sticky to tori of stable families co-existing 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.

The relative importance of non-periodic orbits and the x1mul2 family
Trying to assess the significance of the presence in phase-space of each of the families of p.o. for enhancing a "CX" side-on 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 ∆p z . 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 J = −0.41. It is typical because it is away from bifurcating points and the families of p.o. retain their main stability character. At E J = −0.41 coexist the families x1 (S), x1v1 (S) and x1v2 (U) (Fig. 1). For describing the properties of these orbits which are relevant to the side-on morphology of the peanut structure, we base our description on the (z, p z ) 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 J = −0.41 small perturbations of x1 by ∆p z bring in the system quasi-periodic orbits of cylindrical shape. They are characterised by a local |z| minimum for y = 0 in the (y, z) projections, similar to that of quasi-periodic orbits around x1 (Fig. 2a). We notice though, that at this E J we have a second, local, |z| minimum at x = 0, best observed in the (x, z) projection. This means that we have two characteristic sets of dents along the boundaries of the quasi-periodic orbits. Such an orbit is depicted in Fig. 15a and the local minima under discussion for x = 0, are marked with arrows. The innermost two elliptical curves around (z, p z ) = (0, 0) 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 ∆p z perturbation, always for constant E J = −0.41, both local minima deepen, approaching z = 0. The increase of the ∆p z 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 (0, −0.2) 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 p z , the boundaries of the non-periodic 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 p z = 0.195 x1 orbit we encounter a new kind of orbital behaviour. The first 1500 consequents form an elliptical, lemon-shaped, ring in the (z, p z ) 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 (x, y), (x, z) and (y, z) projections respectively. There is a conspicuous enhancement of an ∞type morphology in the (y, z) 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 (y, z) projection given in Fig. 15e. We underline also the narrowness of the orbit along the x axis. This may be Figure 14. The (z, p z ) projection of the 4D phase space of x1 orbits perturbed in the p z direction, at E J = −0.41. We indicate with black and grey " * " symbols the location of the x1v1 and x1v1 ′ p.o. respectively. Red and green "×" symbols indicate the two branches of x1v2. The " " at (z, p z ) = (0.1, 0.42) marks the location of the p.o. x1mul2. Consequents corresponding to parts of the orbits that are plotted in other figures are marked with ⊙ symbols. crucial for the appearance of the "CX" morphology in the edge-on 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 p z direction, we come very close to the initial conditions of x1v2. However, by varying only the p z coordinate we will never reach the x1v2 p.o. since we keep constant the x 0 initial coordinate, equal with the x 0 value of x1. This is smaller than the x1v2 x 0 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 z 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 (∆z, ∆p z ) = (0.15, 0.2) and (0.2, 0.2) 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 (z, p z ) projection. In Fig. 15f we give the (y, z) 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 p z larger than the p z 0 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 (z, p z ) ≈ (0, 0.3) we find them in azimuthal distances of about 90 • (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 p z 0.5. The morphology of the perturbed by p z = 0.5 x1 orbit is already quite chaotic in the (y, z) 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 lemon-shaped 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 p z direction and by perturbing the x1v2 orbit Figure 15. The evolution of the morphology of quasi-periodic orbits around x1, as we approach the initial conditions of x1v2 at E J = −0.41. (a) A typical quasi-periodic orbit corresponding to the green invariant curve around x1 in Fig. 14. (b,c,d) The sticky orbit that forms the dark blue ring in Fig. 14 after 10 x1 periods. (e) The (y,z) projection of the same orbit after 2249 intersections with the surface of section. It still has an "∞"-type morphology. (f) The side-on view of the orbit coloured light blue in Fig. 14, which is sticky to the x1v1 and x1v1 ′ tori. The overplotted red part corresponds to 10 x1 periods. (g) The (y,z) projection of the x1 orbit perturbed by p z =0.5. It has a chaotic character. Note the different scales of the axes on the panels. 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 (z, p z ) = (0.1, ±0.42) 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 (y, z) 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 J = −0.3976 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 J ≈ −0.479. 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 J ≈ −0.479, where the b 2 index becomes tangent with the b = −2 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 −0.48 E J −0.392 it is stable, while for E J −0.392 it becomes complex unstable. We note that the nonzero coordinates of the initial conditions of the p.o. of the x1mul2 family are x 0 and p z 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 J =−0.41. The chosen 3D projection is (x, z, p z ), while the consequents are coloured according to their p x values.
Side-on profiles of quasi-periodic and sticky orbits around x1mul2 support boxy structures. On top of this the orbits are selfintersecting close to (y, z) = (0, 0) in their side-on views. Most of the quasi-periodic 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 ∆x = −0.111, integrated for 12 periods of x1v2 at E J = −0.4. In Fig. 19a the projec- tion particular structure. However, in Fig. 19b, right panel, we have a conspicuous support of a boxy feature with an "CX" morphology. The face-on view (x, y), 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 quasi-periodic orbits around 3:1 stable orbits on the equatorial plane. Even so, it retains its narrowness in the x direction, especially close to y = 0. 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 ∆p z = 0.2 from the p.o at E J = −0.39. The orbit in Fig. 19c shares the same properties with the one in Fig. 19b in supporting a "CX' side-on profile.

The phase-space 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 J = −0.392 that differs from the p.o. only by ∆x = −10 −5 . 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 (x, p x , p z ) 3D projection and we colour the consequents according to their z values. After 400 intersections two disky objects tend to form centred at (x, p x , p z ) ≈ (0.118, 0, ±0.47), 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 (x, p x ) projection, i.e. if we observe the orbit of Fig. 20a from above. We give the first 100 consequents of the orbit in the (x, p x ) 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 (p z > 0), while black to the lower one (p z < 0). 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 J =−0.392, besides x1mul2 the x1v1 family is also complex unstable. Perturbations of x1v1 by ∆x = 10 −4 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.

A new side-on 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 J 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 side-on profile composed by quasi-periodic and sticky orbits associated with the orbits of x1mul2 is given in Fig. 21a. We have taken five perturbed by p z ≈ 0.025 x1mul2 orbits for E J = −0.46, −0.44, −0.42, −0.40 and −0.39. The x1mul2 p.o. at these energies are stable except for E J = −0.39, which is complex unstable. They form a clear butterfly structure with dimensions (∆y × ∆z) ≈ (1.6 × 1.2), i.e. it reaches a height |z| ≈ 0.6. 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 X-shaped bulge, the "butterfly" will be located in the central area of the profile. How much pronounced it will ap-  Fig. 10c. Both profiles support a "∞"type feature. pear 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.  Fig. 11 projected with θ = 8 • . These are cases in which we have an "∞" or "X" feature away from the side-on views. The orbits in all three frames are given in the same size.

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 edge-on, thus practically we want to investigate just the changes introduced by the rotation of the galaxy around its rotational axis. By rotating the non-periodic orbits that compose the profiles we have presented in Fig. 10, Fig. 12 and Fig. 21 we first realise that non-periodic orbits in the neighbourhood of x1v2 and x1mul2 p.o. contribute to the formation of a b/p or "∞" feature only close to their side-on views. Away from this projection angle they only contribute to the formation of a boxy bulge, remaining confined in the z direction.
As regards the orbits associated with the x1v1 family, besides the side-on view that gives the "OX" profiles ( Fig. 10), there is another projection angle for which non-periodic 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 end-on view of the bar, i.e. we have for the box a ratio r b/p = R max /z max 1. In addition the resulting morphology is closer to an "∞"-type or peanut-shaped structure, than to an X-shaped feature with straight line segments as wings, as is for the side-on views. In Fig. 22a we give a set of x1v1 and x1v1 ′ p.o. in the energy range −0.41 <E J < −0.27, from the point of view (φ, θ)=(90 • ,15 • ) (for the endon view we have (φ, θ)=(90 • ,0 • )). Non-periodic orbits close to the x1v1 and x1v1 ′ p.o. however, do not build b/p profiles with features as sharp as in their side-on 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 θ ≈ 15 • for quasi-periodic orbits after the complex unstable part of x1v1 (Fig. 1). For quasiperiodic orbits before the complex unstable part of the family (E J < −0.3965) however, the "∞" feature is best viewed with θ ≈ 8 • . 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 θ ≈ 15 • is proposed by Gardner et al. (2014) for the viewing angle of the "X-shaped" 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 side-on view of the non- Figure 23. Comparison of a quasiperiodic and a chaotic orbit supporting an "∞"-type morphology close to their end-on views. In (a) we give orbit "7" and in (b) figure "6" of Fig. 11. In both cases we have θ = 13 • periodic 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 J = −0.29 are quasi-periodic, while those in Fig. 23b with E J = −0.31 are chaotic. We realise the degree at which they enhance the "∞"-type morphology close to their end-on views. We note that if we consider all the orbits of the sample of Fig. 11 at a nearly end-on 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 z = 0 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 non-periodic orbits in the neighbourhood of x1v1 that projected resemble structures with a "CX" type morphology. As an example we give the perturbed by ∆z = 0.01 orbit corresponding to the fourth invariant curve around x1 in Fig. 4, starting with p x = 0, 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 J = −0.438225). However, away from the edge-on orientation (in our example in Fig. 24 φ = 60 • ) 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.

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 side-on profiles of disc galaxy models. The families of the x1-tree can be considered as the backbone of the peanut structure in two ways. On one hand they trap around them quasi-periodic 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 x1-tree (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 J interval of the radial-ILR and vertical-ILR. 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 side-on profiles based mainly on a x1v1-x1v1 ′ backbone, or even support a b/p morphology just by itself. In more detail our conclusions can be summarised as follows: (i) We find non-periodic orbits supporting a b/p profile in the side-on view before the critical energy in the v-ILR region at which x1v1 is bifurcated from x1. These orbits are quasi-periodic orbits on tori around stable x1 orbits. Their outline tends to a combined x1v1-x1v1 ′ morphology as E J increases towards the value at which x1v1 is introduced in the system (Fig. 2). This tells us firstly that the peanut starts before the v-ILR 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 quasi-periodic orbits. We do not have an abrupt change of the phase space structure immediately upon the introduction of the bifurcated family in the system.
(ii) Non-periodic 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 (x, p x ) plane. This is a significant dynamical mechanism for building the body of the peanut.
(iii) The neighbourhood of x1 p.o. is a substantial source of nonperiodic orbits that contribute to boxy side-on profiles. They provide such orbits to the system before the v-ILR and in the ∆E vILR 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 −0.332 <E J < −0.302). These orbits contribute to a b/p morphology mainly by becoming sticky to x1v1 and x1v1 ′ rotational tori. Else, they contribute to an overall edge-on boxiness.
(iv) In the ∆E vILR 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.
(v) Edge-on profiles of non-periodic 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 edge-on galaxies. This indicates that if x2 stellar flows exist in 3D stellar bars, they are embedded inside the edge-on profiles of the b/p structures.
(vi) In the case of X-shaped 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 side-on views. Along the wings of "X" we have p z = 0 (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 J 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 X-shaped b/p bulge based on x1v1 dynamics, because in this case this family is complex unstable over smaller ∆E J ranges.
(vii) "CX-type" profiles can be supported: (a) By quasi-periodic orbits around x1 in energies, where the x1v2 family co-exists. They become more important as we increase their p z 0 initial condition from 0 towards the p z 0 value of x1v2. The strongest reinforcement is given by orbits corresponding to the closest to x1v2 tori around x1.
(b) 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.
(viii) 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 J energy intervals. Quasiperiodic orbits trapped around the orbits of this family and, mainly, sticky orbits to them reinforce a b/p structure.
(ix) 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 quasi-periodic 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 co-exist at the same E J , we have found orbits in their immediate neighbourhood that stay in different regions of the phase space for thousands of consequents.
(x) Considering non-periodic 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 side-on views. Away from the side-on view we have in general just a boxy morphology in the outer envelopes of the edge-on projections. However, in the case of X-shaped boxy bulges based on x1v1 dynamics, we find nearly end-on views, where the same orbits can give an "∞"-type morphology.