Chaos and dynamical trends in barred galaxies: bridging the gap between N-body simulations and time-dependent analytical models
aa r X i v : . [ a s t r o - ph . GA ] D ec Mon. Not. R. Astron. Soc. , 000–000 (2013) Printed 17 September 2018 (MN LaTEX style file v2.2)
Chaos and dynamical trends in barred galaxies: bridging the gapbetween N -body simulations and time-dependent analytical models T. Manos , ⋆ and Rubens E. G. Machado CAMTP - Center for Applied Mathematics and Theoretical Physics, University of Maribor, Krekova 2, SI-2000 Maribor, Slovenia. School of Applied Sciences, University of Nova Gorica, Vipavska 11c, SI-5270 Ajdovˇsˇcina, Slovenia. Instituto de Astronomia, Geof´ısica e Ciˆencias Atmosf´ericas, Universidade de S˜ao Paulo, R. do Mat˜ao 1226, 05508-090 S˜ao Paulo, Brazil.
Received 2013 November 17
ABSTRACT
Self-consistent N -body simulations are efficient tools to study galactic dynamics. However,using them to study individual trajectories (or ensembles) in detail can be challenging. Suchorbital studies are important to shed light on global phase space properties, which are theunderlying cause of observed structures. The potentials needed to describe self-consistentmodels are time-dependent. Here, we aim to investigate dynamical properties (regular/chaoticmotion) of a non-autonomous galactic system, whose time-dependent potential adequatelymimics certain realistic trends arising from N -body barred galaxy simulations. We constructa fully time-dependent analytical potential, modeling the gravitational potentials of disc, barand dark matter halo, whose time-dependent parameters are derived from a simulation. Westudy the dynamical stability of its reduced time-independent 2-degrees of freedom model,charting the different islands of stability associated with certain orbital morphologies and de-tecting the chaotic and regular regions. In the full 3-degrees of freedom time-dependent case,we show representative trajectories experiencing typical dynamical behaviours, i.e., interplaybetween regular and chaotic motion for different epochs. Finally, we study its underlyingglobal dynamical transitions, estimating fractions of (un)stable motion of an ensemble of ini-tial conditions taken from the simulation. For such an ensemble, the fraction of regular motionincreases with time. Key words: galaxies: kinematics and dynamics - galaxies: structure – galaxies: evolution –galaxies: haloes – methods: numerical
Orbits are the fundamental building blocks of any galactic struc-ture and their properties give important insight for understandingthe formation and evolution of such structures (Binney & Tremaine1987). Our understanding relies significantly on the adequacy andefficiency of the models used either in the time-dependent (TD)self-consistent models or in the rather ‘simpler’ analytical time-independent (TI) ones. The presence of chaos, manifested as un-stable orbital motion, and expressed by exponential divergence ofnearby trajectories, is broadly studied over the past years (for agood review on chaos in galaxies, see the book by Contopoulos2002). It is by now well accepted that the chaotic or regular na-ture of orbits influences the general stability of the N -body sim-ulations, which is straightforwardly related to the underlying dy-namics. Therefore, studying the general stability and the detailedstructure of the phase (but also of the configuration) space of ana-lytical models can be proven to be very useful, provided that these ⋆ E-mail: [email protected] potentials are realistic in terms of representing density distributionprofiles close to those derived by simulations.The general nature of an orbit, in conservative (TI) sys-tems, can only be one of the following: periodic (stable or unsta-ble), quasi-periodic or chaotic (Lichtenberg & Lieberman 1992).Nevertheless, there are cases where chaos can be characterizedas weak , suggesting that orbits spend a significant fraction oftheir time in confined regimes and do not fill up phase spaceas ‘homogeneously’ as the strongly chaotic ones. In these case,the different rate of diffusion in the phase space plays an im-portant role, associated for example the weak chaotic motionwith barred or spiral galaxy features, giving rise to a numberof interesting results (see e.g. Athanassoula et al. 2009, 2010,2009; Harsoula & Kalapotharakos 2009; Harsoula et al. 2011a,b;Contopoulos & Harsoula 2013; Kaufmann & Contopoulos 1996;Patsis et al. 1997; Patsis 2006; Romero-G´omez et al. 2006, 2007;Tsoutsis et al. 2009; Brunetti et al. 2011; Bountis et al. 2012).There are also several results in the recent literature show-ing that strong local instability does not necessarily implywidespread diffusion in phase space (Cachucho et al. 2010;Giordano & Cincotta 2004). In Contopoulos & Harsoula (2008, c (cid:13) Manos & Machado N -body simulations due to the exchange of angular momentum (seee.g. Athanassoula & Misiriotis 2002; Athanassoula 2003) and (b)a case where the bar gets weaker by losing mass (see e.g. Combes2008; Combes 2011). There, a new reliable way of using the Gener-alized Alignment Index (GALI) chaos detection method was usedfor estimating the relative fraction of chaotic vs. regular orbits insuch TD potentials. We stress here that in the TD models, indi-vidual trajectories may display sudden transitions from regular tochaotic behavior and vice versa during their time evolution and ingeneral the ‘sticky’ behaviour, as discussed in the literature, is lesspronounced. This is also the typical case in the N -body simula-tions where, generally speaking, the motion may also be either: (i)regular throughout the whole evolution, (ii) chaotic throughout thewhole evolution, (iii) alternate between chaotic and regular motionwith simultaneously orbital shape change (but not necessarily), e.g.from disc to bar like, etc...Completing furthermore the picture on previous studies on thestudy of (ir)regular motion in TD systems proposed for cosmologi-cal and galactic models, we could refer to a number of publicationsstarting with the definition of the orbital complexity n ( k ) of anorbital segment which corresponds to the number of frequenciesin its discrete Fourier spectrum that contain a k -fraction of its to-tal power, used by Kandrup et al. (1997); Siopis et al. (1998). Thequantity n ( k ) was later associated with the short-time evolutionof the Lyapunov Exponents (LEs) for TD models in (Siopis et al.1998). A study in a cosmological system where trajectories werechanging their dynamical nature (from regular to chaotic and/orvice-versa) was performed by Kandrup & Drury (1998). The roleof friction, noise, periodic driving, black holes was studied bySiopis & Kandrup (2000) while later on this was extended (inKandrup et al. 2003; Terzi´c & Kandrup 2004) following the tran-sient chaos due to damped oscillations. In Sideris (2009) an ex-ponential function of time was added to the H´enon-Heiles poten- tial, using the so-called ‘pattern method’ as chaos detection tooland more recently the dynamics of some simple TD galactic mod-els was investigated (in Caranicolas & Papadopoulos 2003; Zotos2012).Regarding the galaxies’ evolution and formation of their sev-eral features, it is generally accepted that the most appropriateway to study them is by analyzing N -body simulations. The self-consistency of the models in this approach captures much betterseveral details of the general dynamics. The direct application ofchaos detection methods to individual orbits is still a rather diffi-cult task while, for a large ensemble of particles, it is even harderif not unfeasible. To overcome this obstacle, mean field poten-tials have been used in the literature in order to study in moredetail the dynamical properties of a specific N -body simulation.These potentials are referred to as ‘frozen’ and they are TI andare derived at specific snapshots of the simulations. Hence, onecan apply chaos detection tools to the mean field potential in-stead of the N -body simulation. For example, Muzzio et al. (2005)used an elliptical galaxy simulation (no bar or halo) without dis-sipation which collapses and eventually reaches an equilibriumstate. Then, by taking a quadrupolar expansion of the frozen snap-shot, they derive a stationary smooth potential. In Voglis et al.(2006) and references therein, the authors deal with disc galax-ies, focusing mainly on the spiral structures rather than bars (nohalo) while the extraction of the mean field potential is again per-formed in a similar manner. Following this approach, the role ofchaotic motion and diffusion rate in barred spiral galaxies has alsobeen studied (Harsoula & Kalapotharakos 2009; Harsoula et al.2011a,b; Maffione et al. 2013; Contopoulos & Harsoula 2013)while some applications to the Milky Way bar can be found inWang et al. (2012) and recently a new code for orbit analysis andSchwarzschild modelling of triaxial stellar systems was given inVasiliev (2013). Nevertheless, following this approach one onlyderives a stationary mean field model for an equilibrium state ofa simulation under study. Furthermore, it does not incorporate anappropriate type set of parameters that would be able to describeand reproduce the time-dependencies in axis ratios, masses, patternspeed, etc. of the several components of a model, like for exam-ple the growth of the bar component or the evolution of the disc intime. Let us point out here, the fact that in all these approaches theorbits under study can be only either regular or chaotic. The latterones may be further distinguished to strongly or weakly chaotic,depending on their diffusion properties, sticky effects etc... duringthe whole evolution.In this paper, we consider an N -body simulation of a discgalaxy embedded in a live halo (i.e. both the stellar disc and thedark matter halo are represented by responsive particles). Disc-halointeraction leads to the formation of a strong bar. We then measurehow the galaxy components vary in time during the simulation. Thetime evolution of the structural parameters is provided as input tothe analytical TD model we build. This ‘candidate’ analytical meanfield potential is meant to mimic the N -body simulation evolutionand more importantly to generate orbits with more similar (and insome sense ‘richer’) morphological behaviour to those of the N -body simulation, i.e., permitting for individual orbits the interplaybetween regular and chaotic epochs as time evolves and providinga stable structure at the same time. Note that, in TI frozen modelsan orbit cannot convert from chaotic to regular. Our TD model iscomposed of three components (bar, disc and halo) whose param-eters were fitted with the N -body measurements, chiefly via therotation curves. Note that many simplifying assumptions are made.For example, our TD model considers an (ellipsoidal) analytical bar c (cid:13) , 000–000 haos and dynamical trends in barred galaxies Figure 1. (Colour online) Snapshots of the N -body simulation (upper panels) at four different times, displaying stellar density, on the same range, projectedon the xy plane. Each frame is 40 by 40 kpc. To illustrate bar lengths and shapes, we overlay ellipses (which are not isophotal fits). Rather, their semi-majoraxes are obtained from the radial m = 2 Fourier component of the mass distribution, and from shape measurements via the inertia tensor (see text). The lowerpanels display the result of evolving an ensemble of initial conditions in the presence of the constructed TD analytical potential. component which is not always an excellent approximation of theshape of the actual N -body bar. Likewise, the analytical descriptionof the halo and the disc cannot be expected to behave identically,either. However, our goal is to study the general dynamical impactin stability caused by the bar’s growth in time (as it happens in the N -body simulation). Thus, by using a realistic TD model, withoutaiming to describe of the exact detailed dynamics yielding from thesimulation, we can use chaos detection tools and quantify generaltrends of relative regular and chaos in the phase and configurationspace. Keeping this in mind, we draw (disc) initial conditions di-rectly from the simulation and we evolve in time them with themean field TD potential.The paper is organized as follows: In Section 2 we describethe N -body simulation we used, and we already start presentingthe first part of our results, on the construction of a novel time-dependent analytical model and we present its dynamical proper-ties with respect to this simulation. In Section 3 we give brieflythe definitions of the chaos detection methods employed and theirbehaviors for chaotic and regular motion. In Section 4 we explorethe global phase space dynamics of the derived analytical model,first for the stationary case (frozen potential), e.g. for different (butfixed) sets of parameters which correspond in a sense to differentsnapshots of the simulation. In Section 5 we study the dynamicaltrends and stability of the 3-d.o.f. time-dependent analytical modelfor both a sample of single orbits but also for a large ensembleof initial conditions from the simulation. Finally, in Section 6 wesummarize our findings. N -BODY SIMULATION AND CONSTRUCTION OFTHE ANALYTICAL MODEL Our aim is to construct a time-dependent analytical model thatrepresents the evolution of a barred galaxy. In order to have anastrophysically well motivated model, we supply it with (time-dependent) structural parameters measured from the output of an N -body simulation. N -body simulation To serve as the base reference for the analytical model, we use oneof the simulations described in Machado & Athanassoula (2010).For simplicity, we select initial conditions with a spherical halo.The mass of the stellar disc is M d = 5 × M ⊙ , with an ex-ponential density profile of radial scale length R d = 3 . kpc, andvertical scale height z = 0 . kpc: ρ d ( R, z ) = M d πz R d exp (cid:18) − RR d (cid:19) sech (cid:18) zz (cid:19) , (1)The spherical dark matter halo has a Hernquist (1993) density pro-file and it is five times more massive than the disc: ρ h ( r ) = M h π / αr c exp ( − r /r c ) r + γ ′ , (2)where M h = 2 . × M ⊙ is the mass of the halo, γ ′ = 1 . kpcis a core radius and r c = 35 kpc is a cutoff radius. The normalisa- c (cid:13) , 000–000 Manos & Machado tion constant α is defined by α = { − √ πq exp ( q )[1 − erf ( q )] } − (3)where q = γ ′ /r c . For additional details on the initial conditions,see Machado & Athanassoula (2010).This is a fairly representative collisionless simulation of astrongly barred galaxy. Four snapshots of the disc particles are dis-played in the upper row of Fig. 1. It was performed with the N -body code GYRFALCON (Dehnen 2000, 2002) using a total of 1.2million equal-mass particles, with a gravitational softening lengthof 0.175 kpc, resulting in 0.1 per cent energy conservation. Thesimulation was carried out for approximately one Hubble time. N -body simulations have been employed to study chaotic mo-tion in simplified models of disc galaxies. For example, Voglis et al.(2006) study chaos and spiral structure in rotating disc galaxies, butthose galaxies are not embedded in dark matter haloes.The connection between chaos and bars was also analysed byEl-Zant & Shlosman (2002), with models were set up by the addi-tion of disc, bar and halo components. They found that in centrallyconcentrated models, even a mildly triaxial halo lead to the onsetof chaos and the dissolution of the bar in a timescale shorter thanthe Hubble time. We construct an analytical model that is described by its total grav-itational potential V = V B ( t ) + V D ( t ) + V H ( t ) , where the threecomponents correspond to the potentials of the bar, disc, and halo,respectively. These components will evolve in time, in accordancewith the behaviour we measure from the simulation. Each of thesethree components is represented in the following way: (i) A triaxial Ferrers bar (Ferrers 1877), whose density is givenby: ρ ( x, y, z ) = (cid:26) ρ c (1 − m ) if m < , if m > , (4)where ρ c = π GM B ( t ) abc is the central density, M B ( t ) is the massof the bar, which changes in time, and m = x a + y b + z c , a > b > c > , with a, b and c being the semi-axes of the ellip-soidal bar. The corresponding bar potential is: V B ( t ) = − πGabc ρ c Z ∞ λ du ∆( u ) (1 − m ( u )) , (5)where G is the gravitational constant (set to unity), m ( u ) = x a + u + y b + u + z c + u , ∆ ( u ) = ( a + u )( b + u )( c + u ) , and λ is the unique positive solution of m ( λ ) = 1 , outside of the bar( m > ), while λ = 0 inside the bar. The analytical expressionof the corresponding forces are given in Pfenniger (1984). In ourmodel, the shape parameters (i.e. the lengths of the ellipsoid axes a , b and c are) are also functions of time. (ii) A disc, represented by the Miyamoto–Nagai potential(Miyamoto & Nagai 1975): V D ( t ) = − GM D ( t ) q x + y + ( A + √ z + B ) , (6)where A and B are its horizontal and vertical scale-lengths, and M D ( t ) is the mass of the disc. Here, ‘disc mass’, refers to the stellarmass excluding the bar. As the bar grows, its mass increases at theexpense of the remainder of the disc mass, such that the total stellarmass is constant: M B ( t )+ M D ( t ) = 5 × M ⊙ . The parameters A and B are also functions of time. h a l o a H ( k p c ) . . . . h a l o γ . . . . . . d i s c A , B ( k p c ) AB M D , M B ( M ⊙ ) M D M B a , b , c ( k p c ) abc Ω b ( k m s − k p c − ) Figure 2. (Colour online) Time evolution of the halo, disc and bar param-eters, measured from the N -body simulation (points) and supplied to theanalytical model (fitted polynomials). First and second panel: parametersof the Dehnen halo profile. Third: parameters of the Miyamoto–Nagai disc.Fourth: bar mass and disc mass. Fifth: semi-major axes of the bar. Sixth:bar pattern speed. (iii) A spherical dark matter halo, represented by a Dehnen po-tential (Dehnen 1993): V H ( t ) = GM H a H × − − γ (cid:20) − (cid:16) rr + a H (cid:17) − γ (cid:21) , γ = 2 , ln rr + a H , γ = 2 .(7) c (cid:13) , 000–000 haos and dynamical trends in barred galaxies M H is the halo mass, a H is a scale radius and the dimensionlessparameter γ (within γ < ) governs the inner slope. Thehalo mass is constant throughout, but the parameters a H and γ arefunctions of time. For γ < its finite central value is equal to (2 − γ ) − GM H /a H .Instead of attempting to use the (disc and halo) profiles fromthe N -body simulations, we opted to represent the bar, disc andhalo using respectively the Ferrers, Miyamoto–Nagai and Dehnenprofiles. There are two reasons for such a choice. First, our ap-proach requires analytical simplicity that could not be afforded bythe profiles used in the initial conditions of the numerical simula-tion. Secondly, due to bar formation and evolution, the initial discprofile in the simulation soon becomes a poor representation in theinner part of the galaxy, where the bar resides. In this sense, it isnot advantageous to continue using the initial profiles to model latertimes. A Miyamoto–Nagai disc provides a sufficient approximationfor our purposes. Likewise, even though the Hernquist (1993) haloprofile is well suited for numerical purposes, it is inconvenient fromthe analytical point of view. We experimented with simple logarith-mic halo profiles (because their rotation curves are also appropri-ate), but the Dehnen (1993) profile was preferable, as it providedequally acceptable rotation curves, and a more satisfactory globalapproximation of the mass distribution. Similarly, fitting the bar bya Ferrers ellipsoid is a justifiable approximation. Surely, it fails tocapture the N -body bar in all its complexity, particularly after thebuckling instability, when the bar is substantially strong and devel-ops the peanut-shaped feature. In general, the N -body bar will bemore boxy than the Ferrers shape would allow. Nevertheless, fit-ting an ellipsoid of the same extent allows us to obtain plausibleshapes, to determine the bar orientation and to estimate its massadequately. Ultimately, regardless of small deviations in the den-sity profiles, our goal is to obtain an analytical total potential thatis approximately comparable to the overall potential of the simula-tion. From the simulation, we are able to measure several quantitiesas a function of time, which are then used to inform the analyticalmodel.For the bar, the required parameters are the bar mass, the barshape and the bar pattern speed Ω b . First, we estimate the bar lengthas a function of time. This is done by measuring the relative con-tribution of the m = 2 Fourier component of the mass distributionas a function of radius, for each time step, and finding the radiusat which the m = 2 has its most intense drop after the peak. Thisradius a is associated with the bar length. Then we estimate the barshape by calculating the axis ratios b/a and c/a from the eigenval-ues of the inertia tensor in this region. The bar mass is measuredby simply adding up the mass enclosed within an ellipsoid of axes a, b, c . Finally, the successive orientations of the major axis as afunction of time are used to compute the pattern speed. The re-sulting time evolution of all these quantities are displayed in thefourth, fifth and sixth panels of Fig. 2. Each of these parameters ismeasured at several time steps, a sample of which is shown, alongwith the resulting polynomial fits.The disc mass M D ( t ) is known once the bar mass has beenmeasured, and the halo mass M H is constant. One still requires thetime evolution of two disc parameters ( A , B ) and two halo param-eters ( a H , γ ). This is achieved by measuring the rotation curvesdirectly from the simulation (at each time step), and then fitting theanalytical v c ( R ) to these data. Since the disc and halo potentialsare known from equations (6) and (7), we obtain their respectiveanalytical circular velocities from v c = R dVdR : v c,D ( R ) = R GM D [ R + ( A + B ) ] / (8) v c,H ( R ) = GM H r − γ ( r + a H ) − γ (9)Fitting Eq. (9) to the measured halo rotation curve, we obtain a H and γ . In the case of the disc, it is not enough to fit Eq. (8). Onemust simultaneously fit the Miyamoto–Nagai density profile to dis-ambiguate the A + B (ignoring the inner part of the disc). Whenfitting the disc rotation curve, we assume the total stellar mass (i.e.we take both disc and bar mass into account). Since the circularvelocities rely on azimuthally averaged quantities, the presence ofthe bar does not greatly interfere with the quality of the fits, whileits removal would lead to spurious results. The measured rotationcurves (disc, halo and total), as well as the resulting fitted circularvelocities, are displayed in the fourth column of Fig. 3 (at four illus-trative instants in time). Errors in the fitted parameters of rotationcurves were typically of about 5 per cent or less. Also shown in thefirst, second and third panels of Fig. 3, are the disc (radial and verti-cal) and halo density profiles. The points correspond to simulationmeasurements and the lines give the resulting fits.One of the main arguments in favor of the adequacy of our an-alytical model is evidenced by the fact that its total rotation curvesare in good agreement with those measured from the simulation.This indicates that the choices of profiles were not unreasonable,as they result in a globally similar gravitational potential. Even ifindividually the densities of the components are idealized simplifi-cations, the similarity of the total potential ensures that the overalldynamical evolution should be sufficiently well approximated.Finally, the resulting time evolution of the halo and disc struc-tural parameters, measured in the manner described above, are dis-played in the first to fourth panels of Fig. 2. With these, the time-dependence of the analytical model is fully specified. In Table 1 wesummarize the analytical model by showing a sample of param-eters for the Ferrers bar, Miyamoto–Nagai disc and Dehnen halopotential as fitted by the N -body simulation, at four times. In order to measure the bar strengths in analytical models,Manos & Athanassoula (2011) had employed the Q b parameter(Buta et al. 2003, 2004), which is a measure of the relative strengthof the non-axisymmetric forces. Here, we opt instead to use amethod more familiar to N -body simulations, namely measure-ments of the m = 2 Fourier component of the mass distribution.For the N -body simulation, we measure this component straight-forwardly as a function of radius and then take the maximum am-plitude to be the A (see e.g. Athanassoula, Machado & Rodionov2013). We refer to this quantity as the bar strength.For the analytical model, we proceed in a way that allows usto treat it as if it could be represented by particles. We extract fromthe simulation a random sample of 100 000 initial conditions (i.e.positions and velocities of disc particles) at a time t = 1 . Gyr(see Sect. 5, where we use this ensemble of orbits to further studydynamical trends). The orbits of each of these ‘test particles’ arethen evolved forward in time in the presence of our time-dependentanalytical potential. Their successive positions can be treated as ifthey were simulation particles. By stacking them at each time step,we produce the snapshots in the bottom row of Fig. 1. These mocksnapshots display a striking resemblance to the N -body snapshots, c (cid:13) , 000–000 Manos & Machado
Figure 3. (Colour online) Parameters of the Dehnen halo profile and of the Miyamoto–Nagai disc are obtained by fitting the circular velocity curves at eachtime. Here we display four different times. The fourth column exhibits total, halo and disc circular velocity curves. The first and second columns show the discdensity profiles (radial and vertical, respectively). The third column has the halo density profile. Points come from measurements of the N -body simulation,while lines are fitted profiles. Table 1.
The parameters for the Ferrers bar, Miyamoto–Nagai disc and Dehnen halo potential as fitted by the N -body simulation, at four times.Bar Disc Halotime a b c Ω b M B A B M D a H γ M H (Gyr) (kpc) (kpc) (kpc) (km s − kpc − ) ( M ⊙ ) (kpc) (kpc) ( M ⊙ ) (kpc) ( M ⊙ )1.4 2.24 0.71 0.44 52 1.40 1.92 0.22 3.96 3.90 0.23 254.2 5.40 1.76 1.13 24 2.36 0.95 0.53 2.64 5.21 0.71 257.0 7.15 2.38 1.58 14 3.02 0.78 0.56 1.98 5.77 0.85 2511.2 7.98 2.76 1.93 9 3.30 0.71 0.59 1.70 5.95 0.89 25 specially bearing in mind that they were obtained by very indirectmeans. While this comparison cannot be expected to yield a perfectmorphological equivalence, one notices that the bar lengths are inquite good agreement, and that in both cases rings are present (al-though not of the same extent). The point is that the dynamics thatarises from the analytical model will give rise to very similar discand bar morphologies. In fact, the relative importance of the baris also quite comparable, as indicated by the A parameter. Anal-ogously to the N -body case, we compute the A of these mocksnapshots and compare them in Fig. 4.We must stress here that this comparison is an a posteriori verification, i.e. the bar strength of the N -body simulation was infact not used as an input to the analytical model. The fact the A do agree well counts as a further sign of the consistency of theconstructed analytical model.It is clear, of course, that the variation of the bar strength mod-ifies the values of several parameters and yields richer information about the dynamics of a self-consistent model. N -body simulationsshow that in general, variations of the bar mass also change themass ratios of the model’s components, the bar shape and the pat-tern speed of the galaxy. Hence, if one wishes to use a mean fieldpotential to ‘mimic’ a self-consistent model as accurately as pos-sible, one should allow for all the parameters that describe the bar(together with all other axisymmetric components) to depend ontime, assuming that the laws of such dependence were explicitlyknown. In our case, however, we adopt a simpler approach and varyonly the masses of the bar and the disc, as a first step towards in-vestigating such models when time-dependent parameters are takeninto account. Thus, we do not pretend to be able to reproduce theexact dynamical evolution of a realistic galactic simulation. Rather,we wish to understand the effects of time dependence on the gen-eral features of barred galaxy models and compare the efficiencyof chaos indicators like the Generalized ALignment Index (GALI) c (cid:13) , 000–000 haos and dynamical trends in barred galaxies . . . . A N -body simulationanalytical model Figure 4. (Colour online) Bar strength, measured from the simulation (line)and from the analytical model (points). A is the maximum relative contri-bution of the m = 2 Fourier component of the mass distribution in thedisc. method and the Maximal Lyapunov Exponent (MLE) in helping usunravel the secrets of the dynamics in such problems.We stress that the method we introduced to construct the ana-lytical model does not rely – at all – on frozen potentials. Instead,it is grounded on the detailed features of a fully time-dependent,self-consistent N -body simulation.Unless otherwise stated, the units of the analytical model aregiven as: 1 kpc (length), 1000 km · sec − (velocity), 1 Myr (time), × M J (mass) and km · sec − · kpc − ( Ω b ) while the param-eter G = 1 . The total mass M tot = M B ( t ) + M D ( t ) + M H isset equal to × M J and since the halo’s mass M H is keptconstant, the disc’s mass M D ( t ) is varied as M D ( t ) = M tot − ( M H + M B ( t )) . Let us here, for the sake of completeness, briefly recall how the twomain chaos detection methods used throughout the paper, namelythe GALI and the MLE, are defined and calculated. Consideringthe following TD 3-d.o.f. Hamiltonian function which determinesthe motion of a star in a 3 dimensional rotating barred galaxy: H = 12 ( p x + p y + p z ) + V ( x, y, z, t ) − Ω b ( t )( xp y − yp x ) . (10)The bar rotates around its z –axis (short axis), while the x directionis along the major axis and the y along the intermediate axis of thebar. The p x , p y and p z are the canonically conjugate momenta, V is the potential, Ω b ( t ) represents the pattern speed of the bar and H is the total energy of the orbit in the rotating frame of reference(equal to the Jacobi constant in the TI case).The corresponding equations of motion are: ˙ x = p x + Ω b ( t ) y, ˙ y = p y − Ω b ( t ) x, ˙ z = p z , ˙ p x = − ∂V∂x + Ω b ( t ) p y , ˙ p y = − ∂V∂y − Ω b ( t ) p x , ˙ p z = − ∂V∂z , (11)while the equations governing the evolution of a deviation vector w = ( δx, δy, δz, δp x , δp y , δp z ) needed for the calculation of the MLE and the GALI, are given by the variational equations: ˙ δx = δp x + Ω b ( t ) δy, ˙ δy = δp y + Ω b ( t ) δx, ˙ δz = δp z , ˙ δp x = − ∂ V∂x∂x δx − ∂ V∂x∂y δy − ∂ V∂x∂z δz + Ω b ( t ) δp y , ˙ δp y = − ∂ V∂y∂x δx − ∂ V∂y∂y δy − ∂ V∂y∂z δz − Ω b ( t ) δp x , ˙ δp z = − ∂ V∂z∂x δx − ∂ V∂z∂y δy − ∂ V∂z∂z δz. (12)Regarding the estimation of the value of the MLE, λ , of anorbit under study we follow numerically its evolution in time to-gether with its deviation vectors w , by solving the set of equa-tions (11) and (12) respectively. For this task we use a Runge-Kuttamethod of order 4 with a sufficiently small time step, which guaran-tees the accuracy of our computations, ensuring the relative errorsof the Hamiltonian function (in the TI case) are typically smallerthan − . Furthermore, we need to have a fixed time step in orderto ensure that in the TD case the orbits vary simultaneously withthe potential.In general, the derivatives of the potential V depend explicitlyon time and and the ordinary differential equations (ODEs) (11) arenon-autonomous. Hence, one has to solve together the equations forthe deviation vectors (12) with the equations of motion (11). Trans-forming the Eqs. (11) [and consequently (12)] to an equivalent au-tonomous system of ODEs by considering time t as an additionalcoordinate (see e.g. Lichtenberg & Lieberman 1992, section 1.2b),is not particularly helpful, and is better to be avoided as shown inGrygiel & Szlachetka (1995).So, in order to compute the MLE and the GALI we numeri-cally solve the time-dependent set of ODEs (11) and (12). Then,according to Benettin et al. (1976); Contopoulos et al. (1978);Benettin et al. (1980) the MLE λ is defined as: λ = lim t →∞ σ ( t ) , (13)where: σ ( t ) = 1 t ln k w ( t ) kk w (0) k , (14)is the so-called ‘finite time MLE’, with k w (0) k and k w ( t ) k beingthe Euclidean norm of the deviation vector at times t = 0 and t > respectively. A detailed description of the numerical algorithm usedfor the evaluation of the MLE can be found in Skokos (2010).This computation can be used to distinguish between regularand chaotic orbits, since σ ( t ) tends to zero (following a power law ∝ t − ) in the former case, and converges to a positive value in thelatter. But Hamiltonian (Eq. 10) is TD, which means that its orbitscould change their dynamical behavior from regular to chaotic andvice versa, over different time intervals of their evolution. In suchcases, the MLE (13) does not behave exactly as in TI model (pre-senting in general stronger fluctuations) and its computation mightnot be able to identify the various dynamical phases of the orbits,since by definition it characterizes the asymptotic behavior of anorbit (see e.g. Manos et al. 2013). Nevertheless, we will show theMLE for a number of orbits throughout the paper, for a more globaldiscussion of the several dynamical properties observed.Thus, in order to avoid such problems in our study, we use theGALI method of chaos detection (Skokos et al. 2007). The GALIindex of order k (GALI k ) is determined through the evolution of k N initially linearly independent deviation vectors w i (0) , c (cid:13) , 000–000 Manos & Machado i = 1 , , . . . , k , with N denoting the dimensionality of the phasespace of our system. Thus, apart from solving Eqs. (11), which de-termines the evolution of an orbit, we have to simultaneously solveEqs. (12) for each one of the k deviation vectors. Then, accord-ing to Skokos et al. (2007), GALI k is defined as the volume ofthe k -parallelogram having as edges the k unit deviation vectors ˆ w i ( t ) = w i ( t ) / k w i ( t ) k , i = 1 , , ..., k . It can be shown, that thisvolume is equal to the norm of the wedge product (denoted by ∧ )of these vectors:GALI k ( t ) = k ˆ w ( t ) ∧ ˆ w ( t ) ∧ . . . ∧ ˆ w p ( t ) k . (15)We note that in the above equation the k deviation vectors are nor-malized but their directions are kept intact.In principal and for TI systems, the GALI k ( t ) for regular or-bits remains practically constant and positive if k is smaller or equalto the dimensionality of the torus on which the motion occurs, oth-erwise, it decreases to zero following a power law decay. For thechaotic ones, all GALI k ( t ) tend exponentially to zero with expo-nents that depend on the first k LEs of the orbit (Skokos et al. 2007,2008). Nevertheless, in the TD case studied in Manos et al. (2013)and also here, the way the theoretical estimation of the GALI’s ex-ponential rates are strongly related to the LEs, being more compli-cated and still open to further inquiry.The procedure used for the detection of the several differentdynamical epochs of the TD system is the following: We evolve theGALI k with k = 2 or k = 3 (i.e., using 2 or 3 deviation vectors)for the 2-d.o.f or 3-d.o.f. cases respectively and whenever GALI k reaches very small values (i.e. GALI k − ) we re-initialize itscomputation by taking again k new random orthonormal deviationvectors, which resets the GALI k = 1 . We allow then these vectorsto evolve under the current dynamics. For time intervals where theindex decays exponentially corresponds to chaotic epochs while inthe other cases to non-chaotic. The reason in doing this is that weneed to follow the current dynamical stability of an orbit under-study which in principle can interplay between chaotic and regu-lar for different epochs during the total time evolution. Thus, letus assume that a trajectory experiences chaotic dynamics and lateron drifts to a regular regime. The volume formed by the deviationvectors will first shrink exponentially to very small values and re-main small throughout the whole evolution unless one re-initializesthe deviation vectors and their volume, in order to allow them to‘feel’ the new current dynamics. However, when we are interestedin more general dynamical trends in time (less details), we will fixtime intervals and we will re-initialize the deviation vectors in thebeginning of each one. Shedding some light on the underlying TI dynamics is an importantstep for understanding the more complicated case of the fully 3-d.o.f. TD model where all parameters vary simultaneously in time.By setting z, p z equal to zero at t = 0 (remaining zero at alltimes) in the Hamiltonian Eq. (10), the orbits’ motion is then re-stricted in the 2 dimensional ( x, y ) space. Note that here t = 0 refers to the t = 1 . Gyr of the N -body simulation. We canthen study, in a frozen potential, individual orbits and the stabil-ity of the phase space, in terms of detecting and locating chaoticand regular motion, for the several sets of the potential parametersat deferent times, as derived from the N -body simulation. We shallbegin by choosing fixed parameters from four time snapshots, i.e.,at t = 1 . Gyr, t = 4 . Gyr, t = 7 . Gyr and t = 11 . Gyr and we integrate orbits for 10 Gyr. In Fig. 5, we present the Poincar´e Sur-face of Section (PSS) defined by x = 0 , p x > with H = − . ,for three typical orbits being integrated for 10 Gyr. The set of pa-rameters for the bar, disc and halo components are chosen from thefits with the 3-d.o.f. TD Hamiltonian at t = 7 . Gyr of the N -bodysimulation (see Table 1 for more details). The blue ( ∗ ) points on thePSS correspond to a disc regular orbit, forming a curve by the suc-cessive intersections with the plane x = 0 , having initial condition ( y, p y ) = ( − . , . with x = 0 and p x = H ( x, y, p y ) (called‘DR’ from now on). Its projection on the ( x, y ) -plane is shown inthe top left inset panels of Fig. 5 and coloured in blue. The GALI for this orbit confirms that its motion is regular by oscillating to apositive value during its evolution in time as well as the MLE σ following a power law decay [see second and third top inset panelsof Fig. 5 respectively (blue), note that the axis here are in lin-logscale]. The three small black curves in the central part of Fig. 5,marked with ( × ), are formed by the successive intersections of aninitial condition with ( y, p y ) = ( − . , . (we will call it ‘BR’).From its projection on the ( x, y ) -plane [first bottom inset panel ofFig. 5 (black)], it is evident that it is a bar-like orbit elongated alongthe long x -axis. It surrounds a stable periodic orbit of period 3 in thecenter of these islands and its regular dynamics is clearly revealedfrom the evolution of the GALI and the MLE σ evolution [sec-ond and third top (black) inset panels of Fig. 5 respectively, notethat the axis are again in lin-log scale]. The central scattered redpoints on the PSS, marked with (+) , correspond to a chaotic orbitwith initial condition ( y, p y ) = (2 . , . (called ‘DC’ from nowon). In the three inset panels positioned vertically in the right partof the Fig. 5 (red), we depict its projection on the ( x, y ) -plane (toppanel). Its GALI successive and exponential decrease to zero intime (middle panel) indicates its chaotic nature. Notice that we re-initialize the deviation vectors each time the index becomes small ( − ) . Its MLE σ , as expected, converges to a positive value(bottom inset).Since in this case, the potential has no TD parameters, the gen-eral asymptotic nature can be either regular or chaotic (weakly orstrongly). This of course will change when the parameters start tovary in time and the motion can convert from one kind to the other,though this will be driven by the momentary underlying dynamics.Hence, it is important to have a good idea of how the phase space(through the potential parameters) evolves in time and for differentvalues of the Hamiltonian function. Regarding the former we willpick four sets of parameters given at different times as extractedfrom the N -body simulation, i.e., at t = 1 . Gyr when the bar isstill relatively small, at the end of the simulation ( t = 11 . Gyr)and at two intermediate states ( t = 4 . Gyr and t = 7 . Gyr). Thenext step is to pick some ‘representative’ or ‘most significant’ en-ergy (Hamiltonian function) values from the total available energyspectrum at each time and set of parameters that would be relevantfor the N -body simulation as well. For this reason, we first calcu-lated the energy dispersion of the 100 000 disc particles from thesimulation. The energies are calculated from the TD Hamiltonianfunction Eq. (10) at t = 1 . Gyr. In the inset panel we show thecumulative distribution of this set of orbits. From the energy dis-tribution histogram in Fig. 6, we can see that the majority of thetrajectories of the N -body simulation are concentrated in relativelylarge values.Exploiting this information, we first choose a sample of 2-d.o.f. Hamiltonian function values (for the four times mentionedabove and the respective sets of parameters), from the interval ofenergies where the majority of the N -body simulation’s particles ismore probable to be found. From this sample, we select a subset of c (cid:13) , 000–000 haos and dynamical trends in barred galaxies −1.5−1−0.5 0 0.5 1 1.5−4 −2 0 2 4 6 p y y −4−2 0 2 4 −4 −2 0 2 4 xy (DR) −8−6−4−2 0 0 2 4 6 8 10 time (Gyr)Log (GALI ) (DR) −4−3−2−1 0 2 4 6 8 10 time (Gyr)Log ( s ) (DR) −4−2 0 2 4 −4 −2 0 2 4 xy (BR) −8−6−4−2 0 0 2 4 6 8 10 time (Gyr)Log (GALI ) (BR) −4−3−2−1 0 2 4 6 8 10 time (Gyr)Log ( s ) (BR) −4−2 0 2 4 −4 −2 0 2 4 xy (DC) −10−8−6−4−2 0 2 0 2 4 6 8 10 time (Gyr)Log (GALI ) (DC) −4−3−2−1 0 2 4 6 8 10 time (Gyr)Log ( s ) (DC) Figure 5. (Colour online) The Poincar´e Surface of Section defined by x = 0 , p x > with H = − . , for three typical orbits (two regular and onechaotic) being integrated for 10 Gyr. The set of parameters for the bar, disc and halo components are chosen from the fits with the 3-d.o.f. TD Hamiltonian at t = 7 . Gyr of the N -body simulation. In the insets we depict their projection on the ( x, y ) -plane together with the GALI and MLE σ evolution in time(see Table 1 for the exact parameters and text for more details on these trajectories). % o f o r b it s ( e n s e m b l e fr o m t h e N - body ) energy Cumulative distribution
Figure 6. (Colour online) Dispersion energy histogram for an ensemble of100 000 initial conditions from the N -body simulation. The energies arecalculated from the TD Hamiltonian function Eq. (10) at t = 1 . Gyr. Inthe inset panel we show the cumulative distribution of this set of orbits. representative energies to illustrate typical phase space structures,focusing at this point on the underlying dynamics. Then, we chartthe regular and chaotic regimes of the phase space with GALI .In Fig. 7, we have used a grid of 100 000 initial conditions on the ( y, p y ) -plane of the corresponding PSS and we have constructed a(colour online) chart of the chaotic and regular regions similar tothe PSS, but with more accuracy and higher resolution, in a sim-ilar manner just like in Manos & Athanassoula (2011), using theGALI method. The different colour corresponds to the different final value of the GALI after 10 Gyr (10 000 time units) for orbitsrepresenting each cell of the grid. The yellow (light-grey in b/w)colour corresponds to regular orbits (and areas) where the GALI oscillates around to relatively large positive values, the black colorrepresents the chaotic orbits where GALI tends exponentially tozero (10 − ), while the intermediate colors in the colour-bars be-tween the two represent ‘weakly chaotic or sticky’ orbits, i.e. or-bits that ‘stick’ onto quasi-periodic tori for long times but theirnature is eventually revealed to be chaotic. Note that the modelin this case is TI and hence there is no need for re-initializationof the deviation vectors, since the asymptotic dynamical natureof the orbits does not change in time. The first row refers to aset of potential parameter given at t = 1 . Gyr for Hamilto-nian values H = − . , − . , − . , − . , the second rowat t = 4 . Gyr for H = − . , − . , − . , − . , the thirdrow at t = 7 . Gyr for H = − . , − . , − . , − . and thefourth row at t = 11 . Gyr for H = − . , − . , − . , − . .Using the above approach, we can measure and quantify thevariation of the percentage of regular orbits in the phase space asthe total energy increases for a specific choice of potential param-eters at same fixed times. The chosen values of the Hamiltonianfunctions cover the range of the available energy interval up tothe value of the escape energy which in general is different. Al-though the main general trend is that this percentage decreases asthe energy grows, its behavior changes at high energy values whereit is no more monotonic. Note that this happens for energy val-ues H > − . out of the range of N -body simulation orbits.In Fig. 8 we show the variation of percentages of regular motionas a function of the energy H for the different sets of parameters c (cid:13) , 000–000 Manos & Machado
Figure 7. (Colour online) The (in)stability map for the 2-d.o.f. frozen potential case using a grid of ≈
100 000 initial conditions on the PSS and integrating themfor 10 Gyr. The colour-bar represents the final GALI values of each initial condition in the end of the iteration. The first row refers to a set of potential parame-ter given at t = 1 . Gyr for Hamiltonian values H = − . , − . , − . , − . , the second row at t = 4 . Gyr for H = − . , − . , − . , − . ,the third row at t = 7 . Gyr for H = − . , − . , − . , − . and the fourth row at t = 11 . Gyr for H = − . , − . , − . , − . . The yellow(light-grey in b/w) colour corresponds to regular orbits where the GALI oscillates around to relatively large positive values, the black color represents thechaotic orbits where GALI tends exponentially small values, while the intermediate colors in the colour-bars between the two represent ‘weakly chaotic orsticky’ orbits. The exact set of parameters used at each t = 1 . , . , . , . Gyr are given in table 1. at t = 1 . Gyr, t = 4 . Gyr, t = 7 . Gyr and t = 11 . Gyr.The threshold GALI > − was used to characterize an orbitas regular and GALI < − as chaotic which will also be thechaos criterion/threshold throughout this paper. We should empha-size that these percentages refer to a set of initial conditions thatcover uniformly the whole 2-d.o.f. phase space. On the other hand,an ensemble of trajectories extracted from the N -body simulationdoes not necessarily populate ‘democratically’ the phase space. Bysimply inspecting these percentages one can claim that the fractionof regular motion is systematically larger for later times, and for allenergies. When looking and comparing the phase space for earlytimes, i.e., t = 1 . , . Gyr (first and second row) and late times,i.e., t = 7 . , . Gyr (third and fourth row) in Fig. 7, we maysee that the central island of stability, originating bar-like orbits, isbecoming larger as the time grows and this is even more evident forthe relatively larger energies (see the third and fourth row from topto bottom). This indicates that the bar component becomes gradu-ally more important and dominant.Considering the N -body simulation’s bar growth evolution(Fig. 1), this fact seems to be already in a quite good agreement,since one would expect the barred morphological features to be-come more evident as time increases and for the relevant energies,in both the simulation and the derived TD analytical model. Wemay here note that the energies of the N -body simulation parti-cles (being 3 dimensional objects) are concentrated in the interval − . . H . − . (as shown in Fig. 6). A more ‘realistic’ en-semble of initial conditions will be used for the 3-d.o.f. case in thenext section. Despite the fact that these percentages refer to the 2-d.o.f. case, one may get a brief idea of how the relative stability ofthe phase space changes as the energy varies in the ‘full’ 3-d.o.f.model and expect a relatively large fraction of regular motion foran ensemble chosen within this energy interval. For the global study of the dynamics of the phase and config-uration space of a 3-d.o.f. model (TI and/or TD), the choice ofinitial conditions plays a crucial role. When one seeks to ex-plore the whole available phase space the orbits might occupy, afew useful approaches have been proposed and used in the lit-erature (e.g. Schwarzschild 1993; Papaphilippou & Laskar 1998;El-Zant & Shlosman 2002), where some adequate (but different)ways in populating several different families of orbits were pro-posed by giving initial conditions on several appropriate chosenplanes in positions and momenta with the total energy restrictiontaken into account (also used recently in Maffione et al. 2013). InManos & Athanassoula (2011) and Manos et al. (2013) the authorsused two similar distributions as in the latter references and also athird one, derived by a random set of orbits related to the densitymass distribution of the model. Their momenta were set zero along c (cid:13) , 000–000 haos and dynamical trends in barred galaxies % o f r e gu l a r m o ti on energy t = t = t = t = Figure 8. (Colour online) Percentages of regular motion for the 2-d.o.f.frozen case as a function of the energy H for the different sets of parametersat t = 1 . Gyr, t = 4 . Gyr, t = 7 . Gyr and t = 11 . Gyr. the y, z directions while the p x was estimated by the total Hamil-tonian function value, randomly chosen from the available energylevel. The approaches described above may supplementary explorewell the phase space of a model. However, they do not ensure thatthe particles of a realistic N -Body simulation necessarily populatethe phase space in the same uniform way. This becomes even morecomplicated when the potential under consideration ceases to be TIand the parameters’ time-dependencies in time cause consecutivealternations in the phase space.For these reasons, here, and in order to study global stabilityand dynamical trends in our 3-d.o.f. TD analytical model, we ex-tract an ensemble of 100 000 initial conditions directly from the N -body simulation at the time t = 1 . Gyr where the bar has al-ready started to be formed and starts growing from that point on.Note that these orbits are chosen only from the disc density dis-tribution of the simulation and do not include halo particles. Ourgoal is twofold: (i) to estimate quantitatively the fraction of chaoticand regular motion as time increases and (ii) to check whether andhow this variation is associated now to the bar strength for this TDmodel. From this point and on, we attribute the time t = 1 . Gyrto t = 0 . Then, in the TD analytical model, the orbits are inte-grated for Gyr (a bit less than a Hubble time), a realistic upperlimit related to the N -body simulation’s set up.Let us point out here that the straightforward comparisonof the two approaches is rather difficult, since several assump-tions have been made for the construction of the TD analyticalmodel. For example, the total energy of the particles is not ex-pected to be conserved, like in the N -body simulation where theself-consistency ensures this with a good accuracy. Moreover, haloparticles, present in the N -body simulation, are not taken into ac-count as point mass particles at this point. The reason, as explainedearlier, is that we are mainly interested in the bar’s growth and evo-lution in time. Hence, when one allows all the TD model param-eters of the several components of model to vary in time for a setof initial conditions, the total energy is not conserved. Neverthe-less, and for the sake of a more general study, we tried two alterna-tive cases; let us call them dynamical scenarios whose parametersare given in such a way that the total energy can be conserved bymaking reasonable assumptions. We then checked how sensitivethe general dynamical trends are to these choices of potentials. Scenario A : In this case, all the parameters of the three compo-nents (bar-disc-halo) of the potential vary simultaneously in time, following the fitting functions derived from the simulation. In thiscase the value of the Hamiltonian in general is not conserved intime for a single trajectory. Here, the TD model incorporates thetime-dependencies of the N -body simulation for each particle,however in the latter the sum of the energy of all the particlesremains approximately constant. This is the evolutionary scenariopresented in the panels of the lower row in Fig. 1. We cansee that this TD model can capture quite successfully the basictrends and the bar formation of the N -body simulation (upper row). Scenario B : In this case, the parameters of the bar and disc varyas in case ‘Scenario A’ but the halo ones are adjusted in such a waythat the value of the total Hamiltonian function, for each initialcondition, remains constant in time. In order to achieve that, wefirst allow the bar and disc potential to vary in time as before,keeping constant the halo parameter γ ( ≈ . at t = 1 . Gyr)and then we allow the a H to vary at each time step. From theinitially known Hamiltonian value, we calculate the value of thetotal potential at every time step and then find the value of the V H ( t ) such as the total energy remains constant. Having thisvalue, we can numerically estimate the appropriate value of a H ( t ) that fulfills the imposed energy conservation condition. This alsoimplies that the scale radius of the halo might vary in a similar wayfor all the orbits but not identically so. Scenario C : In this case, the bar parameters are completely time-dependent, the disc parameters are fixed in an approximately aver-aged value, with respect to the whole evolution, i.e., A = 0 . and B = 0 . , while the halo parameters are estimated just like in case‘Scenario B’, obtaining again the energy conservation.Before focusing on the general dynamical trends of the abovedifferent scenarios , let us first study the evolution of a few typi-cal examples of orbits. We picked a random but representative (interms of radial distance) number of orbits from our 100 000 ini-tial conditions and calculated their orbital evolution together withtheir GALI and MLE. Let us mention here that for all three Sce-narios the number of escaping orbits is negligible. From this sub-set of orbits, we chose to show here, in Fig. 9 and Fig. 10, twotypical examples of trajectories which present also a rather ‘rich’and interesting morphological behaviour. Each figure is divided inthree main blocks, each one corresponding to the different evolu-tion Scenarios A, B and C respectively. Then, each block is furthercomposed of three parts where the orbital evolution is depicted onthe ( x, y ) -plane in the course of time (first rows of each subpart),together with the GALI and MLE σ in the next two rows respec-tively.Thus, in Fig. 9 we show the evolution of an orbit fromthe ensemble of the N -body simulation with initial condition ( x, y, z, p x , p y , p z ) ≈ ( − . , . , − . , . , − . , . (we will refer to this orbit fromnow on as ‘ B ’), which is iterated for 10 Gyr (10 000 time units) forthe TD potentials mentioned above. The B orbit belongs to a setof initial conditions which is representative of the imposed scenarioof the bar growth by our TD potential. The starting and completeset of parameters for the TD model is taken at t = 1 . Gyr bythe fits with the N -body simulation with the procedure describedin Section 2. Note that in all figures’ panels, we set everywhere the t equal to zero instead of t = 1 . Gyr.In the first row of the first block in Fig. 9(a,b,c,d), we showits projection on the ( x, y ) -plane for four successive time intervalsof ∆ t = 2 . Gyr when evolved by the ‘Scenario A’, i.e., all the c (cid:13) , 000–000 Manos & Machado
Scenario A (bar-like orbit B ) -8-6-4-2 0 L og ( GA L I ) (e) -3-2-1 0 2 4 6 8 10 L og ( s ) time (Gyr) (f) Scenario B (bar-like orbit B ) -8-6-4-2 0 L og ( GA L I ) (k) -3-2-1 0 2 4 6 8 10 L og ( s ) time (Gyr) (l) Scenario C (bar-like orbit B ) -8-6-4-2 0 L og ( GA L I ) (q) -3-2-1 0 2 4 6 8 10 L og ( s ) time (Gyr) (r) Figure 9. (Colour online) The 3-d.o.f. orbit B evolved with the Scenario A, B and C . Its different projection on the ( x, y ) -plane in different time windowsis depicted in the first (four-panel) row (from top block-part to the bottom respectively) and the colour bar corresponds to the time (in Gyr). Their GALI andMLE σ evolution in time is shown for each case just below them. Note that the orbit (in all three cases) starts as a disc-like and gradually its shape turns tobarred, displaying the bar’s growth through the parameters of the Ferrers’ potential (see text for more details). c (cid:13) , 000–000 haos and dynamical trends in barred galaxies Scenario A (disc-like orbit D ) -8-6-4-2 0 L og ( GA L I ) (e) -3-2-1 0 2 4 6 8 10 L og ( s ) time (Gyr) (f) Scenario B (disc-like orbit D ) -8-6-4-2 0 L og ( GA L I ) (k) -3-2-1 0 2 4 6 8 10 L og ( s ) time (Gyr) (l) Scenario C (disc-like orbit D ) -8-6-4-2 0 L og ( GA L I ) (q) -3-2-1 0 2 4 6 8 10 L og ( s ) time (Gry) (r) Figure 10. (Colour online) Same as in Fig. 9 but for the 3-d.o.f. disc-like orbit D evolved again with the ‘Scenario A, B and C’. Note, that here the disc-likepattern slightly varies from cases to case. The different degree of chaoticity can be accurately captured by the frequency and fast decay to zero of the GALI [Fig. 10(e,k,q)], indicating that the orbits is relatively ‘strongly chaotic’ under the ‘Scenario A, C’ while under the ‘Scenario B’ is relatively ‘weakly chaotic’.This information can not be revealed in such a way by the MLE σ shown in Figs. 10(f,l,r) (see text for more information).c (cid:13) , 000–000 Manos & Machado three potential components V B , V D , V H are time-dependent andthe total energy is not in general conserved. The colour bar nextto each panel corresponds to the time (in Gyr), hence the most re-cent epochs of the orbit are coloured with yellow (light-grey in b/w)while those in the earlier ones with dark blue or black (dark-greyor black in b/w). In Fig. 9(e) we show its GALI , capturing accu-rately the chaotic nature of the orbit during the first [Fig. 9(a)], third[Fig. 9(c)] and fourth [Fig. 9(d)] time windows by decaying expo-nentially to zero. On the other hand, in the second [Fig. 9(b)] timewindow its regular (even by just looking its projected morphologi-cally on the ( x, y ) -plane) behaviour is successfully revealed by thefluctuates to a non-zero value of the index. Note that the plot is inlin-log scale and the deviation vectors are re-initialized, by takingagain k new random orthonormal deviation vectors, each time theGALI becomes very small (i.e. GALI − ). It turns out thatthe B begins as a regular disc-like orbit during the first 2.5 Gyrand, as the bar starts forming and growing, it gradually evolves toa chaotic bar-like orbit until the end of the integration. We may no-tice how hard it is for the finite time MLE σ [Fig. 9(f)] to capturethese different dynamical different transitions and epochs due to itstime-averaged definition (see also Manos et al. 2013). Furthermore,its power law decay for regular time intervals and its tendency topositive values for chaotic ones are of the same order of magni-tude making it rather hard to use the temporary value σ as a safecriterion of regular and weak or strong chaotic motion.In order to see how morphologically sensitive our ‘full’ TDparameter model is (as described in ‘Scenario A’) to the severalcomponent parameters and the lack of energy conservation, let usevolve the same initial condition ( B ) with the ‘Scenario B’. Wecan see that, the general shape of the orbits is more or less sim-ilar, starting again as an disc-like orbit and drifting to a bar-likeone in time [first row of the second block in Fig. 9(g,h,i,j)]. Ofcourse, the dynamical epochs are a bit different. In the first rowof the third block in Fig. 9(m,n,o,p), we show its evolution whenintegrated with the time-dependencies described in the ‘ScenarioC’ and recovering again similar morphological behaviour. In theelongated panels, beneath the projections on the ( x, y ) -place of the‘Scenario B, C’, we depict their corresponding GALI [Fig. 9(k,q)respectively] and MLE σ [Fig. 9(l,r) respectively] just like for the‘Scenario A’.Of course, we should not expect such a good agreement in theseveral morphological behaviours in general and for all orbits ofthe initial ensemble. More specifically, and in cases where chaosis strong from the very early moments, the orbits are expected todiffer significantly depending on the evolutionary scenario. As alsodiscussed in Carpintero & Wachlin (2006), the chaotic behaviouris sensitive to the choice of the potential, even if it is frozen like intheir case.In Fig. 10, and in a similar manner as in Fig. 9,we show another characteristic disc-like orbit formost of the total of integration with initial condition ( x, y, z, p x , p y , p z ) ≈ ( − . , − . , . , . , − . , . (we will refer to this orbit fromnow on as ‘ D ’). The evolutionary scenarios are again the sameas before, i.e., in the first row of the first block in Fig. 10(a,b,c,d),we present its projection on the ( x, y ) -plane for different timewindows. The D orbit experiences a regular epoch during its first2.5 Gyr, then gradually becomes chaotic switching to a bar-likeshape and finally becomes a chaotic but disc-like now orbit. Itsregular and chaotic epochs are accurately captured by the GALI [Fig. 10(e)], fluctuating to constant value for the first 2.5 Gyr andthen successively decaying exponentially to zero for the rest of the integration. The MLE σ [Fig. 10(f)] also reveals this dynamicalevolution, by decaying with a power law for the regular part andconverging to non-zero value for the three last time windows.However, here the motion does not present any further transitionand/or interplay between regular and chaotic motion and again (asfor the B orbit) the order of magnitude for the σ is not varyingsufficiently enough to lead to a safe conclusion at certain timeswithout seeing its whole time evolution.Its evolution with the ‘Scenario B’, is shown in the first rowof second middle block in Fig. 10(g,h,i,j)] where it turns out thatpresents different morphological shape. In more detail, it still be-gins similarly but for later times it forms a disc-like orbit (lesschaotic also) which does not visit the central region compared towhat happens in ‘Scenario A’. Moreover, when looking at its shapewhen integrated with the ‘Scenario C’ [Fig. 10(m,n,o,p)], we canagain conclude that it is a chaotic disc-like orbit but not identicalto the other cases. Notice that, the different degree of chaoticity canbe accurately captured by the frequency and fast decay to zero ofthe GALI [Fig. 10(e,k,q)], indicating that the orbits is relatively strongly chaotic under the ‘Scenario A, C’ while under the ‘Sce-nario B’ is relatively weakly chaotic . This information can not berevealed in such a way by the MLE σ [Fig. 10(f,l,r)].Furthermore, let us clarify the following for such TD systems.When speaking more strictly and using more rigorous notions fromnon-linear dynamics theory, those orbits experiencing such an in-terplay between regular and chaotic epochs in different time inter-vals are in general (asymptotically) chaotic. Their MLE is positive(like in our examples in Fig. 9 and Fig. 10). However, when oneis interested in astronomical time scales, one may split the time insmaller windows and see how an ensemble of orbits evolve. There,one may detect the different orbital trends (general at first), in termsof chaotic or regular by checking the local exponential divergence .The GALI method is one good way to capture these phenomenaquite successfully, after giving enough time to the deviation vec-tors to detect the stability of the local dynamics. Since one cannotknow or predict in advance, when exactly such transitions betweenregular and chaotic motion will occur for each individual trajectory,one is limited to the study of rather average trends.One common approach broadly performed in the recent litera-ture is to measure the relative fraction of regular and chaotic motionin a sequence of frozen snapshots (e.g. Valluri et al. 2010). How-ever, this would not allow us to capture an important and abundantorbital behaviour in the N -body simulations, namely orbits that al-ternate their nature from regular to chaotic and vice-versa as well astheir morphological shape, e.g. from disc-like shape to barred-likeone. Turning now to the aspect of the global (in)stability of themodel and this set of initial conditions in particular, one wouldbe interested to monitor the variation of the total fraction of reg-ular vs. chaotic motion in the course of time. Since our model isTD the percentage of chaotic orbits following the potentials (de-scribed as Scenarios A, B and C) is expected to change in time. Inorder to estimate this we adopt the following strategy: We dividethe total integration time of 10 Gyr in four successive time win-dows of length ∆ t = 2 . Gyr time units. At the beginning of eachtime window, we re-initialize GALI to unity and follow a newset of three orthonormal deviation vectors for each orbit. Then foreach time window we calculate the current percentage of regular(non-chaotic) orbits as the fraction of orbits whose GALI remains > − during that time interval, i.e., GALI does not decay ex-ponentially (fast) to small values. In this way we allow the GALI to have enough evolutionary time to reveal the chaotic or regular c (cid:13) , 000–000 haos and dynamical trends in barred galaxies
75 80 85 90 95 100 0 2.5 5 7.5 10 % o f r e gu l a r o r b it s time (Gyr) Scenario AScenario BScenario C
Figure 11. (Colour online) Percentages of regular motion for different timewindows. nature of the orbit within this interval. The results of this procedureare shown in Fig. 11, where we depict the percentages of regularmotion for different time windows for the ensembles of initial con-ditions evolved under the three different evolutionary processes.Though, it turns out that all three dynamical evolutionary casesshow quite similar trends, in terms of fraction of chaotic and regularmotion, implying that the variations of the halo and disc parame-ters effect is relatively weak with respect to the bar’s. The amountof regular motion for the three different cases systematically growsin time, starting from ≈ to ≈ after . Gyr reaching ≈ to ≈ after Gyr. The ‘Scenario A’ turns out to beslightly more chaotic compared to other two. Moreover, let us stressthe fact that the presence of regular motion is abundant in this TDmodel when comparing to similar studies (Manos & Athanassoula2011; Manos et al. 2013) where chaos turned out to be dominant.The next step is to see how the bar strength is relatedto the general (average) stability. As found and discussed inManos & Athanassoula (2011), and for TI modes, one should ex-pect to find a strong correlation between the strong bars and thetotal amount of chaos. This result was confirmed for a simple TDmodel, where only the mass of bar was considered to grow in timein Manos et al. (2013), in a potential also composed of a bar and adisc but of a bulge instead a halo component. By simply looking atthe percentages of regular motion in Fig. 11, one would conjecturethat the relative non-axisymmetric forcings are decreasing slightlyin time, implying that the bar gets weaker.However, this is not what happens here. As shown in Fig. 4,the A parameter (maximum relative contribution of the m = 2 Fourier component of the mass distribution in the disc) increasesin time during the first part of the evolution and stays more or lessconstant till the end, for both the N -body simulation and the TDanalytical model. In order to interpret appropriately the increase ofregularity with the simultaneous increase of the non-axisymmetricforcings, we should look back to the 2-d.o.f. case and the PSSsin Fig. 7. There, one may see that for the relative energy intervalvalues of the ensemble of 3-d.o.f. initial conditions (chosen fromthe N -body simulation) the stable island , associated to the bar-likeorbits, gets larger in the course of time. Combining this with thefact that both the self and non-self consistent models enhance thebarred morphological feature as time grows, we may conclude thatthe fraction of regular motion increase in time is linked to the un-derlying growth of the central ‘barred’ island of stability. Moreover, it implies that these initial conditions populate with higher proba-bility this part of the (6 dimensional) phase space of our 3-d.o.f.model which tends in time to encompass larger regular area withbar-like orbits within.Thus, it becomes evident that for a TD model the relative frac-tion of regular (or chaotic) motion is not straightforwardly corre-lated to the bar’s strength. In order to understand such dynamicaltrends, one has to shed some light to the underlying time-dependentdynamics. This can be done by considering how the ensemble isdistributed in the energy interval, and the specific dynamical trendsof the model. These trends are manifested by specific morpholog-ical properties, such as growth of the ‘barred’ stable island in ourcase. They may also affect the global (in)stability in different man-ners depending on the model. In order to carry out all the analyses summarized below, we em-ployed an analytical model that was specially tailored for this pur-pose. It was based on the results of a self-consistent N -body simu-lation of an isolated barred galaxy. We measured several structuralparameters of this simulation, as a function of time, and then usedthem to set up the analytical gravitational potential of a galacticmodel. This model was composed of three components, represent-ing the disc, the bar and the dark matter halo. It implements ana-lytical potentials that are meant to be conveniently simple, whileat the same time being able to mimic the simulation to a reason-able degree of accuracy. It should de pointed out that our analyticalmodel is not based on frozen potentials, in any way. Instead, it isa fully time-dependent model, that relies upon the detailed featuresof the N -body simulation. This entails that, contrary to the frozen(and more accurate) ones, where the variety of orbital motion is re-stricted and destined to be either regular or chaotic, our TD modelis equipped with some extra orbital behavior, namely trajectorieswhich may alternate nature, behaving regularly for some epochsand chaotically for others, and in a not necessarily in a monotonicway. This is also what broadly happens in the N -body simulations.The adequacy of the model we constructed is verified in atleast three essential ways. First, the similarity of the rotation curvesensures the global dynamics should be well approximated. Second,if we compute the orbits of an ensemble of test particles subject tothis potential, they give rise to morphological disc and bar featuresremarkably similar to those of the N -body simulation. Third, thelength and strength of the bar in the resulting mock snapshots (fromthe analytical model) are in very good quantitative agreement withthe N -body bar. Such comparisons indicate that our model is ableto adequately capture both the dynamics and the morphology of thebarred galaxy model in question.Starting with the reduced 2-d.o.f. and time-independent caseof the model, we used the GALI method to successfully survey theunderlying dynamics of the phase space by mapping the chaoticand regular regimes (and motion). By measuring these, we foundthat the fraction of regular orbits increases in time, i.e., for modelparameters at later times, for almost all energies. Moreover, theisland of stability associated to the bar-like trajectories also getslarger in time, implying that the bar features are enhanced as timeevolves. This tasks allowed us to get a brief idea also for the full3-d.o.f. time-dependent as well as the N -body simulation, exhibita bar growth evolution.Regarding this TD analytical model, we similarly estimated stability trends in terms of estimating the amount of regular and c (cid:13) , 000–000 Manos & Machado chaotic motion in different time-windows. In this case, we used amore realistic set of initial conditions coming directly from the sim-ulation itself and iterated them under the constructed TD potential.In order to do that we firstly monitored the GALI’s detection ef-ficiency to a representative sample of orbits and presenting in thispaper two of them which show some typical generic behaviour. Wealso discussed its advantages with the traditional MLE in capturingsuch sudden dynamical variations. Of course we did not manage tospan the whole orbital richness of our ensemble of initial conditionsbut we have rather achieved to give a flavor of the possible evolu-tions for individual trajectories. It turns out that the complete set oforbits tends to become relatively more ‘regular’ in time. To furtherexamine this, and since now the total energy is not conserved, wetried two alternative but similar models whose Hamiltonian func-tion value can be conserved by adjusting the disc’s and/or the haloparameters appropriately. In all cases the trend was found to be thesame.Moreover for this 3-d.o.f. TD model presented here (‘ScenarioA’), we envisage that the detailed study of orbital dynamics can beextended further by, for example, classifying the disc and bar orbitsfrom a morphological point of view as the time varies. Addition-ally we can do the same for the other two ‘Scenarios’ described inSection 5. Furthermore, we have already a quite large ensemble oforbits to study diffusion properties in the phase (also in configura-tion) space, just like in other papers in the literature, and compareit with the one observed in the N -body simulation.Approaches such as this are potentially suited to a broad classof applications in galactic dynamics (e.g. double bars, central massconcentrations and even bar dissolution). Once the output of an N -body simulation has been modeled into a time-dependent an-alytical potential, a variety of analyses could then be undertaken,particularly regarding orbital studies. Work that generally relied onhighly simplified (and usually frozen) analytical potentials couldtake advantage of more astrophysically realistic galaxy models.This bridging would afford an approximation to the richness of de-tail of an N -body simulation, at a lower computational cost andwith the versatility of a simple analytical formulation. ACKNOWLEDGEMENTS
We would like to thank Ch. Skokos, T. Bountis and M. Romero-G´omez, for their comments and fruitful discussions. REGM ac-knowledges support from FAPESP (2010/12277-9). TM was sup-ported by the Slovenian Research Agency (ARRS) and partiallyby a grant from the Greek national funds through the Opera-tional Program ‘Education and Lifelong Learning’ of the Na-tional Strategic Reference Framework (NSRF) - Research Fund-ing Program: THALES, Investing in knowledge society throughthe European Social Fund. This work has made use of the com-puting facilities of the Laboratory of Astroinformatics (IAG/USP,NAT/Unicsul), whose purchase was made possible by the Brazilianagency FAPESP (grant 2009/54006-4) and the INCT-A. We alsoacknowledge support from FAPESP via a Visiting Researcher Pro-gram (2013/11219-3), and we would like to thank the Institut HenriPoincar´e for its support through the ‘Research in Paris’ programand hospitality in the period during which part of this work tookplace.
REFERENCES
Athanassoula E., 2003, MNRAS, 341, 1179Athanassoula E., Machado R. E. G., Rodionov S. A., 2013, MN-RAS, 429, 1949Athanassoula E., Misiriotis A., 2002, MNRAS, 330, 35Athanassoula E., Romero-G´omez M., Bosma A., Masdemont J. J.,2009, MNRAS, 400, 1706Athanassoula E., Romero-G´omez M., Bosma A., Masdemont J. J.,2010, MNRAS, 407, 1433Athanassoula E., Romero-G´omez M., Masdemont J. J., 2009,MNRAS, 394, 67Benettin G., Galgani L., Giorgilli A., Strelcyn J.-M., 1980, Mec-canica, 15, 9Benettin G., Galgani L., Strelcyn J.-M., 1976, Phys. Rev. A, 14,2338Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton Uni-versity PressBountis T., Manos T., Antonopoulos C., 2012, Celest. Mech. Dyn.Astron., 113, 63Brunetti M., Chiappini C., Pfenniger D., 2011, A&A, 534, A75Buta R., Block D. L., Knapen J. H., 2003, AJ, 126, 1148Buta R., Laurikainen E., Salo H., 2004, AJ, 127, 279Cachucho F., Cincotta P. M., Ferraz-Mello S., 2010, Celestial Me-chanics and Dynamical Astronomy, 108, 35Caranicolas N. D., Papadopoulos N. J., 2003, A&A, 399, 957Carpintero D. D., Wachlin F. C., 2006, Celestial Mechanics andDynamical Astronomy, 96, 129Combes F., 2008, in J. Funes & E. Corsini ed., Formation andEvolution of Galaxy Discs Vol. 396 of ASP conference series.p. 325Combes F., 2011, Mem. S.A.It.Suppl., 18, 53Contopoulos G., 2002, Order and Chaos in Dynamical Astron-omy. Springer-Verlag, BerlinContopoulos G., Galgani L., Giorgilli A., 1978, Phys. Rev. A, 18,1183Contopoulos G., Harsoula M., 2008, Int. J. Bif. Chaos, 18, 2929Contopoulos G., Harsoula M., 2010, Celestial Mechanics and Dy-namical Astronomy, 107, 77Contopoulos G., Harsoula M., 2013, MNRASDehnen W., 1993, MNRAS, 265, 250Dehnen W., 2000, ApJL, 536, L39Dehnen W., 2002, J. Comput. Phys., 179, 27El-Zant A., Shlosman I., 2002, ApJ, 577, 626Ferrers N. M., 1877, Quart. J. Pure Appl. Math., 14, 1Giordano C. M., Cincotta P. M., 2004, A&A, 423, 745Grygiel K., Szlachetka P., 1995, Acta Phys. Pol. B, 26, 1321Harsoula M., Kalapotharakos C., 2009, MNRAS, 394, 1605Harsoula M., Kalapotharakos C., Contopoulos G., 2011a, MN-RAS, 411, 1111Harsoula M., Kalapotharakos C., Contopoulos G., 2011b, Int. J.Bif. Chaos, 21, 2221Hernquist L., 1993, ApJS, 86, 389Kandrup H. E., Drury J., 1998, Annals of the New York Academyof Sciences, 867, 306Kandrup H. E., Eckstein B. L., Bradley B. O., 1997, A&A, 320,65Kandrup H. E., Vass I. M., Sideris I. V., 2003, MNRAS, 341, 927Katsanikas M., Patsis P. A., 2011, Int. J. Bif. Chaos, 21, 467Katsanikas M., Patsis P. A., Contopoulos G., 2011, Int. J. Bif.Chaos, 21, 2321Katsanikas M., Patsis P. A., Pinotsis A. D., 2011, Int. J. Bif. c (cid:13) , 000–000 haos and dynamical trends in barred galaxies Chaos, 21, 2331Kaufmann D. E., Contopoulos G., 1996, A&A, 309, 381Lichtenberg A. J., Lieberman M. A., 1992, Regular and ChaoticDynamics, 2nd edn. No. 38 in Applied Mathematical Sciences,Springer-Verlag, New York, NYMachado R. E. G., Athanassoula E., 2010, MNRAS, 406, 2386Maffione N. P., Darriba L. A., Cincotta P. M., Giordano C. M.,2011, Celestial Mechanics and Dynamical Astronomy, 111, 285Maffione N. P., Darriba L. A., Cincotta P. M., Giordano C. M.,2013, MNRAS, 429, 2700Manos T., Athanassoula E., 2011, MNRAS, 415, 629Manos T., Bountis T., Skokos C., 2013, J. Phys. A: Math. Theor.,46, 254017Manos T., Skokos C., Antonopoulos C., 2012, Int. J. Bif. Chaos,22, 1250218Miyamoto M., Nagai R., 1975, PASJ, 27, 533Muzzio J. C., Carpintero D. D., Wachlin F. C., 2005, CelestialMechanics and Dynamical Astronomy, 91, 173Papaphilippou Y., Laskar J., 1998, A&A, 329, 451Patsis P. A., 2006, MNRAS, 369, L56Patsis P. A., Athanassoula E., Quillen A. C., 1997, ApJ, 483, 731Pfenniger D., 1984, A&A, 134, 373Romero-G´omez M., Athanassoula E., Masdemont J. J., Garc´ıa-G´omez C., 2007, A&A, 472, 63Romero-G´omez M., Masdemont J. J., Athanassoula E., Garc´ıa-G´omez C., 2006, A&A, 453, 39Schwarzschild M., 1993, ApJ, 409, 563Sideris I., 2009, Chaos Analysis Using the Patterns Method..p. 347Siopis C., Eckstein B. L., Kandrup H. E., 1998, Annals of the NewYork Academy of Sciences, 867, 41Siopis C., Kandrup H. E., 2000, MNRAS, 319, 43Skokos C., 2010, Lect. Notes Phys., 790, 63Skokos C., Bountis T., Antonopoulos C., 2008, Eur. Phys. J. Spec.Top., 165, 5Skokos C., Bountis T. C., Antonopoulos C., 2007, Physica D, 231,30Terzi´c B., Kandrup H. E., 2004, MNRAS, 347, 957Tsoutsis P., Kalapotharakos C., Efthymiopoulos C., ContopoulosG., 2009, A&A, 495, 743Valluri M., Debattista V. P., Quinn T., Moore B., 2010, MNRAS,403, 525Vasiliev E., 2013, MNRAS, 434, 3174Voglis N., Stavropoulos I., Kalapotharakos C., 2006, MNRAS,372, 901Wang Y., Zhao H., Mao S., Rich R. M., 2012, MNRAS, 427, 1429Zotos E. E., 2012, New Astron., 17, 576This paper has been typeset from a TEX/ LaTEX file prepared by theauthor. c (cid:13)000