Chaotic mixing and the secular evolution of triaxial cuspy galaxy models built with Schwarzschild's method

1 MB·English

by  E. Vasiliev

Mon.Not.R.Astron.Soc.000,1–14(2011) Printed4January2012 (MNLATEXstylefilev2.2)
Chaotic mixing and the secular evolution of triaxial cuspy
galaxy models built with Schwarzschild’s method
E. Vasiliev1,2,3⋆, E. Athanassoula3
2 1Lebedev Physical Institute, Leninsky prospekt 53, Moscow, Russia
1 2Rochester Institute of Technology, 76 Lomb Memorial drive, Rochester, NY 14623, USA
0 3Laboratoire d’Astrophysique de Marseille (LAM), UMR6110, CNRS/Universit´e de Provence, 38 rue Joliot Curie, 13388
2
Marseille C´edex 13, France

n
a
J
August8,2011

3

] ABSTRACT
A
WeusebothN-bodysimulationsandintegrationinfixedpotentialstoexplorethesta-
G bilityandthelong-termsecularevolutionofself-consistent,equilibrium,non-rotating,
. triaxialspheroidalgalacticmodels.Morespecifically,weconsiderDehnenmodelsbuilt
h
with the Schwarzschild method. We show that short-term stability depends on the
p
degree of velocity anisotropy (radially anisotropic models are subject to rapid devel-
-
o opment of radial-orbit instability). Long-term stability, on the other hand, depends
r mainly on the properties of the potential, and in particular, on whether it admits a
t
s substantial fraction of strongly chaotic orbits. We show that in the case of a weak
a density cusp (γ = 1 Dehnen model) the N-body model is remarkably stable, while
[
the strong-cusp (γ = 2) model exhibits substantial evolution of shape away from tri-

1 axiality, which we attribute to the effect of chaotic diffusion of orbits. The different
v behaviour of these two cases originates from the different phase space structure of
7 the potential; in the weak-cuspcase there exist numerous resonantorbitfamilies that
6 impede chaotic diffusion. We also find that it is hardly possible to affect the rate of
6
this evolution by altering the fraction of chaotic orbits in the Schwarzschild model,
0
whichisexplainedbythefactthatthechaoticpropertiesofanorbitarenotpreserved
1. by the N-body evolution.There are, however,parameters in Schwarzschildmodelling
0 thatdoaffectthestabilityofanN-bodymodel,sowediscusstherecipeshowtobuild
2 a ‘good’ Schwarzschild model.
1
: Key words: stellar dynamics – galaxies: structure – galaxies: kinematics and dy-
v
namics – galaxies: elliptical – methods: numerical – methods: N-body simulations
i
X
r
a
1 INTRODUCTION namicalinstabilities,suchastheradial-orbitinstability(e.g.
Polyachenko& Schukhman 1981), which manifest them-
The Schwarzschild method is an important tool for con-
selves on a rather short time, typically of the order of a
structing self-consistent equilibrium models of galaxies,
crossingtime.Theotherreason concernsmoregradualevo-
whenthesystemdoesnothaveananalyticdistributionfunc-
lutionandisrelatedtotheexistenceofchaoticorbitsinmost
tion (Schwarzschild1979).Models oftriaxial galaxies areof
non-integrablepotentials,especiallythosehavingahighcen-
particularinterest,sincetheyoftenadmitalargefractionof
tralmassconcentration (densitycusporblackhole).Orbits
chaotic orbits.
intheSchwarzschildmodel(SM)areregardedasstationary
The aim of the Schwarzschild method is to construct
‘building blocks’, which ensures that the distribution func-
self-consistentequilibriummodels,inwhichorbitsinagiven
tion of the whole model satisfies the time-independent col-
potentialarearrangedsothattheresultingdensitymatches
lisionless Boltzmann equation (Binney & Tremaine 2007),
the potential via the Poisson equation. On the other hand,
df/dt=0. By construction, orbits are evolved for a certain
dynamicalstabilityofthesemodelsisoutofthescopeofthe
time, typically of order 102 dynamical times, which is ex-
Schwarzschild method.
pected to be sufficient to sample the available phase space.
There are two reasons why such a model might turn
Thisisareasonableassumptionforregularorbits,whichfill
out to be non-stationary. First this could be due to dy-
theirinvariant tori more or less uniformly duringthis time.
But chaotic orbits have a much larger region of available
phase space (of higher dimensionality) and may sample it
⋆ E-mail:[email protected](EV),[email protected](EA)
(cid:13)c 2011RAS
2 E. Vasiliev, E. Athanassoula
in a very non-uniform way, remaining in a confined portion times, despite having large (& 50%) fraction of chaotic or-
of this region for many hundreds of dynamical times. This bits. All these results are based on N-body modelling, and
may lead to the effect of chaotic diffusion, when these or- the methods used for creating their initial conditions can-
bits eventually escape to a different region of phase space not make a system with predefined properties, unlike the
and,accordingly,changetheirshapeinconfigurationspace, iterativemethodofRodionov et al.(2009),discussedinSec-
which leads to the breakdown of the self-consistency of the tion 7.
model. Poon & Merritt (2004) constructed models for triaxial
The stability of models built with Schwarzschild’s scale-free cusps with γ = 1 and γ = 2 around a black hole
method has been tested by N-body simulations as early as using the Schwarzschild method. Their solutions contained
inSmith & Miller(1982),albeitwitharatherroughN-body about half of the mass in chaotic orbits, but nevertheless
code and on a short timescale. Zhao (1996) confirmed the were found to be reasonably stable when evolved by an N-
stabilityofagalacticbarmodelusingtheself-consistentfield body tree-code during 6T (measured at few times the
dyn
(SCF)method.Morerecently,Antoniniet al.(2009)showed black hole influence radius), except for prolate (T = 0.75)
that the radial-orbit instability (ROI) exists also in triax- models. They conclude that their chaotic orbits were suf-
ial systems, if the velocity anisotropy coefficient β & 0.3. ficiently mixed during the T = 100 orbital times used for
Wu et al. (2009) explored the existence and stability of tri- integration in SM,so that theyrepresented reasonably sta-
axial Dehnen models in MOND gravity and found it to be tionarybuildingblocks.However,thetimeintervalwasquite
acceptable. Nevertheless, their N-bodymodels initially dis- short (only afew dynamicaltimes for orbitsoutside thera-
played some evolution towards a more equilibrium configu- diusof influence,where most chaoticorbits arefound),and
ration, presumably due to difficulties associated with con- the change of axial ratios was small yet non-negligible. In
structing self-consistent models in MOND gravity. addition, they followed the evolution of models in a fixed
On the other hand, the effect of chaotic diffusion on smooth potential for 100 dynamical times, and found no
∼
the evolution of a SM is usually assumed to lead towards changeinshape,asisreasonabletoexpect,giventhatorbits
amoresphericalmass distribution.Supportforthisconjec- in the SM were evolved for a similar time. Another reason
ture comes from thefact that the ‘average’ shape of a fully forthelessapparentshapeevolutionintheirmodelsisthat
chaotic orbit (which fills almost ergodically the equipoten- forscale-freepotentialstheequipotentialsurfaceisnotmuch
tial surface, excludingtheparts of phase space occupied by rounderthantheequidensitysurface,soeven‘fully chaotic’
regularorbits)isgenerallyrounderthantheequidensitysur- orbitsdosupportthenecessarymodelshapetosomedegree.
face. Schwarzschild (1993) confirmed this effect for a scale- AsomewhatdifferentissueisaddressedinValluri et al.
freelogarithmic potential,corresponding toaγ =2density (2010). They investigated the change of shape of triaxial
cusp, by creating a self-consistent SM and then following dark matter halos in response to the growth of a compact
the ensemble of orbits for an interval of time three times massinthecentre,thatis,duetoanadiabaticchangeofthe
longer than was used in creating the model. He recorded potential. They compared theorbit population of themod-
the shape obtained from the superposition of these orbits els before the growth of the central mass, after the growth
for this longer interval and concluded that it was evolving (intermediatestage),andafteranadiabatic‘evaporation’of
towards sphericity, but that this change was quitesmall. thismass(finalstage).TheyusedanN-bodytree-codeboth
A number of studies used N-body simulations to ad- to follow the evolution of the ‘live’ system and to analyze
dress the long-term stability of triaxial models created by thepropertiesoftheorbitsinthe‘frozen’N-bodypotential.
various methods, and their evolution caused by chaotic or- Valluri et al. (2010) found that, in the case with no strong
bits.Holley-Bockelman et al.(2001)constructedanequilib- central mass concentration, the evolution of orbit shapes is
rium triaxial model with axial ratio a:b:c=1:0.85:0.7 mostly reversible, despite the fact that many of the orbits
based on the spherical Hernquist model (γ = 1) by apply- are chaotic in the intermediate stage, and attributed this
ing artificial squeezing along two axes, and evolved it with reversibility to resonant trapping of orbits in the course of
a SCF N-body code for 8 half-mass dynamical times to theslow change of the potential.
∼
confirm that thereis nosignificant change in shape. We continue and extend these studies in two intercon-
Later, Holley-Bockelman et al. (2002) extended their nected aspects. Namely, we perform N-body simulations of
study to include a supermassive black hole with mass M• triaxial cuspy Dehnen models built with the Schwarzschild
equal to 0.01 of the total model mass. The growth of this method, and find that, in the cases when there is no rapid
central point mass destabilizes thepopulation of box orbits onset of radial-orbit instability, they are remarkably stable
andconvertsmostofthemintoratherstronglychaoticones, overmany dynamical times. There is, however, a long-term
which,inturn,quicklydrivestheinnerregionsofthemodel evolutionoftheshapeofthesemodelsandwefindthatitis
towards almost spherical shape. Again using SCF N-body caused bythe influenceof chaotic orbits.
code,Kalapotharakos et al.(2004)exploredthedependence We begin by introducing our triaxial Dehnen models
of suchevolution on theblack holemass, and found thatin constructed with the Schwarzschild method, and briefly re-
all cases thereexisted alarge fraction of chaotic orbits, but view their properties in section 2. Then in section 3 we re-
theoverallshapeofthemodelsubstantiallyevolvedonlyfor view the concept of chaos, in particular, thedistinction be-
M• > 0.005, when these orbits were more strongly chaotic tween regular, weakly (sticky) chaotic and strongly chaotic
(asmeasuredbyLyapunovexponents),whilesmallerM•did orbits, and the quantities that are intended to represent
not cause much secular evolution. Muzzio et al. (2009) cre- chaotic properties of an individualorbit.
atedastronglytriaxialmodelwithaγ 1cuspbycoldcol- The long-term evolution of orbits in a fixed potential
≃
lapseandconfirmeditsstabilitybyquadrupolar(arestricted is the subject of section 4. It turns out that in the poten-
variant of SCF) N-bodycode overseveral hundredcrossing tials considered here, most chaotic orbits are in fact quite
(cid:13)c 2011RAS,MNRAS000,1–14
Chaotic mixing and the evolution of triaxial galaxies 3
sticky, and the timescale for chaotic diffusion and the as- (whichalsofixesthetimeunit)1,andchoosetheaxialratios
sociated change of orbit shapes is rather long, of order 103 b/a = p5/8 0.79, c/a = 1/2, which corresponds to a
dynamical times. However, in the strong-cusp model and triaxiality para≈meter T = (a2 b2)/(a2 c2) equal to 1/2
− −
in the outer parts of the weak-cusp model this stickiness is (maximal triaxiality)2. For the cusp slope γ we choose two
less prominent, and the chaotic orbits do exhibit evolution values:γ =1(weak cusp) and γ =2 (strongcusp).It turns
towards a more spherical shape. We also discuss how our out that there are substantial differences between the two
findings can be explained by the amount of ‘complexity’ of models, which will be discussed in thefollowing sections.
the phase space, that is, the presence and importance of All integration times and orbit frequencies in the SM
resonant orbit families. are measured in units of dynamical time T (E), which is
dyn
Next, we describe our N-body experiments which defined as the period of the long-axis orbit with the given
basically confirm the expectations derived from the energy: T = 4 rmax(E)[2(E Φ(r))]−1/2dr. Other pa-
dyn R0 −
fixed-potential evolution (section 5). Unless the velocity pers often use the period of the x y plane closed loop
−
anisotropy inthemodel issufficientlybiased towardsradial orbit,whichissomewhat shorter,butourdefinitionhasthe
velocities to allow for a rapid development of the radial- advantage that all natural frequencies of orbits are greater
orbit instability, the shape of the density distribution in thanT−1.ForourmodelsT (r)maybeapproximatedas
dyn dyn
the N-body model (NM) remains in agreement with that 4.4(r1/2+r3/2) for the weak-cusp case and 4.4(r2+r3)1/2
of the original Schwarzschild model (SM) for a long time. for the strong-cusp case, where time and length units are
The weak-cusp model is essentially stable over the course dimensionless time units.
of the simulation, while for the strong-cusp case there is Themodelswereconstructedusing50radialshells,2400
a gradual evolution towards a more spherical (or, rather, gridcellsand3 104 orbitsevolvedfor500dynamicaltimes.
oblateaxisymmetrical)shape,whichweattributetochaotic Theinitialcond·itionsfortheorbitswereassignedrandomly
diffusion. using the following recipe: sample the elliptical radius uni-
In addition, a very important, and somewhat disap- formly in the value of enclosing mass, randomly choose an-
pointing, finding is that we are hardly able to control the gles to put the orbit onto the equidensity ellipsoid, and as-
degree of evolution due to this chaotic diffusion by alter- signeachcomponentofthevelocityaccordingtoaGaussian
ing the properties of the SM. This is because the orbits distribution with the dispersion of the equivalent spherical
in theNM,while basically resembling theircounterpartsin model. This is different from the traditionally used method
the SM in shape and orbital class, do not inherit the at- of sampling initial conditions at a small number of energy
tributeofchaoticityfrom thesmooth-potentialmodel.This levelsandin fixedpositions on thegrid; wefindthat aran-
argues that the overall evolution of the NM is determined dompositionand,moreimportantly,continuousdistribution
mostly bythe gross properties of the potential, rather than in energy yields bettermodels.
by any specific arrangement of orbits (as long as this sat- In the following three sections we will discuss mainly
isfies self-consistency). Nevertheless, we do present recipes onevariantforeachmodel,namelythe‘unconstrained’vari-
forbuildingbetterSM (inthesensethatthecorresponding ant with no restrictions on the orbit population or on the
NM are more stable) in section 6. fraction of chaotic orbits, and with a velocity anisotropy
In section 7 we compare the evolution of models built β = 1 σt2 (Binney & Tremaine 2007) varying from 0 in
with the Schwarzschild method with that of a model con- − 2σr2
thecentreto 0.6intheouterparts.Weconsidertheeffect
structed with the iterative method (Rodionov et al. 2009), ∼
of variation of these parameters in section 6.
whichiscompletelydifferentfromtheSM asitreliesonthe
‘guided evolution’ of an N-body model towards a specific
equilibrium.Weshow that,despiteconceptually beingvery
different, they perform similarly in terms of stability of the
3 GENERAL REMARKS ON THE CHAOTIC
resulting model, and have similar properties of ensemble of
PROPERTIES OF AN ORBIT
orbits.
Finally, we present our conclusions. There are several methods for quantifying the degree of
chaoticity of agiven orbit,which may beclassified intotwo
broad groups. One group deals with a single orbit and con-
siders its spectrum, that is, the Fourier transform of some
quantity along the trajectory (e.g. the x coordinate, or the
distance from the centre) sampled at equally spaced mo-
ments of time. It is based on the fact that all regular or-
2 SCHWARZSCHILD MODELS
bits are multiply-periodic, with their spectra containing in
Inthispaperwerestrictourattentiontonon-rotatingthree- 3D only linear combinations of no more than three funda-
dimensional systems with mild flattening. Namely, we con- mentalfrequencies,thesefrequenciesbeingconstantsof the
sidertwovariantsofthetriaxialDehnenmodel,withdensity motioninatime-independentpotential.Anydeviationfrom
profile
(3 γ)M 1
ρ(r)= − , (1) 1 Note that for dimensionless radius we use the long-axis scale
4πabc mγ(1+m)4−γ
radius a and not the half-mass radius, the latter being 2.4a for
γ=1andaforγ=2models,measuredalongthexaxis.
wherem=[(x/a)2+(y/b)2+(z/c)2]1/2istheellipticradius. 2 T =0correspondstooblateandT =1toprolateaxisymmetric
WeadoptdimensionlessunitsinwhichM =1,a=1,G=1 models,soT =1/2iscalled‘maximallytriaxial’.
(cid:13)c 2011RAS,MNRAS000,1–14
4 E. Vasiliev, E. Athanassoula
this simple spectrum is an indication of chaos. A quantita- the threshold separating regular from chaotic orbits, which
tiveestimateofchaosmaybeobtainedeitherasthenumber roughly corresponds to the distinction based on the Lya-
of spectral lines containing a specified fraction (say, 0.9) of punov exponent, if both are measured on the interval of
the total power (Kandrup et al. 1997), or as the variation 100T .
dyn
in these fundamental frequencies calculated e.g. from the Unfortunately,FDRitselfisnotastrictlydefinedquan-
firstandsecondhalvesoftheintegrationtime(Laskar1993; tity: if we measure ∆ω for two successive integration inter-
Valluri & Merritt 1998), dubbed ‘frequency diffusion rate’ vals,orevenforasomewhat differentdurationoftime(say,
(FDR). 100 and 110 T ), or change the sampling rate, this quan-
dyn
Thesecondgroupofmethodsconsidersthedeviationof tity may change by a factor of few5. This is not surprising,
nearbyorbits;theorbitinquestionisaccompaniedbyoneor since this quantity by definition measures the difference in
moreadjacentorbits,andtheevolutionofdeviationvectors two ’instantaneous’ values of a fluctuating variable ω, and
istracked.Inthecaseofaregularorbitthesedeviationvec- is itself a random quantity with an uncertainty of about
torsshouldgrownofasterthanlinearly,andifmorethanone 0.3 0.5 orders of magnitude.
−
such vector is considered, they should remain correlated in The shape of an orbit in configuration space may be
direction.Conversely,inthechaoticcasethedeviationstarts described byits inertia tensor components,
to grow exponentially after some time. This class includes
manedthaoldigsnbmaesnedt ionndetxheesc(aelc.gu.laStkioonkoosf2L0y1a0p)u.nTohveerexpisonqeunittseΛa Iij ≡ N1s XNs x(in)x(jn), (3)
n=1
good correspondence between FDR and Λ as chaos indica-
tors in a smooth potential. To construct SM with different wherex(in) isthei-thcomponentofthen-thsamplingpoint
fraction of chaotic orbits (Section 6), we use here the Lya- of the trajectory. If we consider only diagonal components,
punovexponentstodistinguishbetweenregularandchaotic √Iii givesan estimate of theextentof thisorbit in thei-th
orbits. coordinate. Likewise, thequantities
For an N-body system, however, all trajectories have I
S ii , i=x,y,z (4)
large positive Lyapunov exponents. Their timescale for i≡ Ixx+Iyy+Izz
exponential deviation is a fraction of the crossing time,
describe the orbit shape (flattening in the i-th direction).
and furthermore, they do not decrease with increas-
Tracingtheirchangeovertimemaybeusedtoestimatethe
ing N (Kandrup & Sideris 2001; Sideris & Kandrup 2002;
Valluri & Merritt 2000)3. This does not preclude a more shape evolution.
Finally, it is necessary to note that theaccuracy of de-
regular behaviour of the system with more particles, since
terminationoforbitalfrequenciesislimitedbytheaccuracy
Lyapunovexponents,bydefinition,measurethegrowthrate
ofenergyconservation(typically∆ω ∆E/E ).Whilethis
of infinitely small perturbations, but do not tell by what ≃| |
isnotalimitationinthecaseofasmoothtime-independent
distance two nearby orbits will be separated after a fi-
potential (the error in energy conservation in the integra-
nite time. It appears that if the phase space has a com-
tor is . 10−9), it becomes a major obstacle when we come
plex structure of stable resonant islands, then chaotic or-
to analyzing the orbits in a live N-body simulation, where
bits usually tend to be confined to small regions of the en-
particlesexperiencerandomvariationsinenergy,e.g.dueto
tire phase space, demonstrating the so-called phenomenon
two-body relaxation.
of stickiness, and resemble regular orbits for many dynam-
ical times (e.g. Contopoulos 1971, 2002; Valluri& Merritt
2000;Harsoula & Kalapotharakos 2009).
Therefore, for the purpose of comparing the chaotic 4 SECULAR SHAPE EVOLUTION INDUCED
properties of orbits evolved in a smooth potential and in BY CHAOTIC ORBITS
anN-bodymodel,wewillrestrictourselvestothefirstclass
TheconstructionofequilibriummodelsbytheSchwarzschild
of chaos detection methods. Namely, we use the following
methodrelies on findingtime-independent‘building blocks’
definition for thefrequency diffusion rate (FDR):
asrequiredbyJeans’theorem.Regularorbitsobviouslysat-
1 3 ω(1) ω(2) isfy this requirement (as long as we ensure sufficiently uni-
∆ω X | i − i | (2) form coverage of theinvariant torus), butfor chaotic orbits
≡ 3 (ω(1)+ω(2))/2
i=1 i i it’s not the case. It has been suggested that all chaotic or-
bitsofagivenenergymaybeinfactconsideredasrepresen-
Here ω(1) and ω(2) are the leading frequencies in Cartesian
i i tatives of one ‘super-orbit’ (in 3D, the chaotic part of the
coordinates(i=x,y,z)forthefirstandthesecondhalvesof
phase space at a given energy is one interconnected region,
theintegrationtime,respectively4.Weadopt∆ω=10−3 as
theso-calledArnoldweb),andhence,ifweaverageallthese
3 AmethodtoestimatetrueLyapunovexponentforanorbitina forthefirsthalf,theotherforthesecondhalf,andthedifference
frozen-N-body system has been proposed by Kandrup&Sideris willbelarge.(Ofcourse,thisstillmeansthatthespectrumexpe-
(2003), but it is unclear whether it can be easily applied to live riences changes, but notnecessarilythat rapidly).Soweusethe
simulations. averagedvalueforallthreecoordinates,andfurthermore,discard
4 This is different from what was used in Valluri&Merritt thelineswithrelativevariationgreaterthan0.5(whichareindeed
(1998), where only the largest of the three differences ∆ωi was rare).
tracked.Wefindthattheirdefinitionmaysometimesexhibitun- 5 However, if the sampling rate is changed by an integer fac-
wantedfluctuations, whenthespectrumhastwodistinctlinesof tor, the frequencies are likely to remain the same, as tested in
comparable amplitudes: one may become the leading frequency Vallurietal.(2010).
(cid:13)c 2011RAS,MNRAS000,1–14
Chaotic mixing and the evolution of triaxial galaxies 5
orbits, the resulting building block may also be used as a
0.8
time-independent unit. In practice, however, the usefulness
S
x
of this approach is doubtful, for a number of reasons. First 0.7
of all, it is not easy to distinguish unambiguously regular
0.6
fromchaoticorbits,becausesomechaoticorbitsmaybeex-
tremelystickyandconservetheirshapeandfrequenciesover 0.5
many orbital periods. Secondly, such ‘fully mixed’ models, 0.4
in which all chaotic orbits with the same energy are aver-
0.3
aged, cannot be built self-consistently (Merritt & Fridman
1996; Siopis 1999): the resulting super-orbit is featureless 0.2
androunderthanthedensitydistribution,andthereforenot 0.1
very suitable for themodel, butregular orbits alone do not bin
0
have a sufficient variety to support the shape, especially in 1 2 3 4 5 6 7 8 9 10
theouter parts of Dehnen models.
0.8
If, on theother hand,we choose to treat chaotic orbits S
x
in the same way as regular ones in constructing SM, the 0.7
model is no longer guaranteed to be stationary. It means
0.6
that even in a constant potential, the shape of the mass
0.5
distribution of a given chaotic orbit may slowly evolve in
time,intheprocessofchaoticdiffusion.Thisshouldleadto 0.4
a global evolution of the model shape, but it is not easy to
0.3
predict how strong thiswill be and in what sense.
0.2
In general, we may expect that a fully chaotic orbit
fills all available configuration space within its correspond- 0.1
bin
ingequipotentialsurface,whichisrounderthantheequiden-
0
sity surfaces of thetriaxial model. Ifwe somehow divideall 1 2 3 4 5 6 7 8 9 10
orbits into three classes – regular, sticky and fully chaotic
– then in the course of time there will be exchange of or- Figure1. ShapechangeforchaoticorbitsduringT =1000Tdyn.
bits between the last two classes: some sticky orbits may Eachpairoferrorbarsshowsthespreadindistributionofthex-
eventually become unstuck, and vice versa. The outcome axis flattening Sx (Eq. 4) for initial (left, red) and final (right,
of this ‘diffusion’ depends on the initial orbit population of green) 100Tdyn intervals of time, averaged over the ensemble of
the model: if initially there were few strongly chaotic or- chaotic(∆ω>10−3)orbitsinagivenenergybin.Thehorizontal
bits (which is to be expected in a triaxial self-consistent axis corresponds to 10 bins, each of which contains 10% of the
totalmass,withtheinnermostparticlesintheleftbin.
model, since such models are not very suitable to support
Toppanel:γ=1,bottom:γ=2models.Thedecreaseintheav-
the shape of the density distribution), then we may expect
eragevalueineachpairmeansthatorbitsbecomerounder(both
that their fraction increases in time, and the overall shape
y-andz-axiscomponents increaseattheexpense ofx-axiscom-
becomes rounder. However, the timescale for this chaotic
ponent).Itisevidentthatintheweak-cuspcaseonlythechaotic
mixing may be rather long, and the change of the overall orbits in the outer shells do change shape systematically to be-
shape of density distribution does not necessarily happen comerounder,whileinthestrong-cuspcasethistendency exists
on this timescale. Merritt & Valluri (1996) find that an en- formostoftheradialshells.
semble of points initially confined to a small region in the
chaotic part of the phase space evolves to a nearly time-
invariant mixed distribution in 100T . However, such
dyn
∼
mixingdoesnotnecessarily distributeorbitsuniformly over between these two intervals, as described in Sect. 3). The
all the accessible chaotic part of phase space – it may well systematic decrease of the x-flattening (and corresponding
be confined to a region surrounded by resonant tori, which increaseinbothyandz,notshownhere)meansthatorbits,
significantly slow down further diffusion (Valluri& Merritt onaverage,becomerounder.Butitisimportanttonotethat
2000). this effect is not observed for all energies: in the weak-cusp
To study this further we perform the following experi- case, only the chaotic orbits in the outer 25% of energy
∼
ment:we create N-bodyrealizations of SM with 104 equal- binswerefoundtoexhibitasubstantialevolution,whilefor
weight particles representing a self-consistent solution, in thestrong-cusp case most chaotic orbits became rounder.
whichorbitswereintegratedfor100T ,Thenwecontinue The difference between these cases may be attributed
dyn
toevolvetheminthesamepotentialfor10subsequentinter- tothepresenceof arichnetwork ofresonances intheγ =1
valsof100T ,andcomparethepropertiesofeachorbitin model. Fig. 2 shows frequency maps for the two models.
dyn
thefirst and thelast interval. The results show that indeed In the weak-cusp model there is a considerable fraction of
thechaoticorbitsshowthetendencytobecomerounder.We points lying on or near resonant lines, corresponding to
may quantify this by measuring the change in flattening of either regular or sticky chaotic orbits. Consequently, they
an orbit in each direction by using S (eq. 4). Fig. 1 shows have rather low values of FDR (right panel, blue), with
i
the mean value and spread of flattening along the x axis, ∆ω . 10−2. On the other hand, γ = 2 model, as well as
measured for the first and the last interval of 100T , av- the outer parts of γ = 1 model, have no major resonances
dyn
eraged overtheensembleof chaoticorbits (definedasthose (apartfrom1:1x andz axistubesand1:2x zbananaor-
that changed their leading frequencies by more than 10−3 bits).Most non-t−ubeorb−itsinthesecasesare−quitestrongly
(cid:13)c 2011RAS,MNRAS000,1–14
6 E. Vasiliev, E. Athanassoula
1.0 w /w 1.0 w /w
y z 4,-3,0 0,1,-1 y z 0,1,-1 bits
0.9 2,0,-1 12,,--23,,11 0.9 2,0,-1 1,-2,1 on of or
0.8 1,-1,0 0.8 1,-1,0 fracti
00..67 5,-1,-23,-1,-14,-2,-1,5-3,0 3,0,-221,1,,-22,-23,1,-32,2,-3 wx /wz 00..67 3,-1,-1 2,1,-2 2,2,-3 wx /wz
0.5 0.6 0.7 0.8 0.9 1.0 0.5 0.6 0.7 0.8 0.9 1.0 10-8 10-6 10-4 10-2 Dw 1
Figure 2. Frequency maps for the γ = 1 (left) and γ = 2 (middle) model. Blue dots mark regular and red dots chaotic orbits (as
distinguishedbythevalueofLyapunovexponent). TherightpanelshowstheFDRdistribution(∆ω)forthesetwocases(blue–γ=1,
green–γ=2), separately fortube-likeorbits(shortaxistubes, longaxis tubes andchaotic orbits closeto1:1resonances, solidlines)
andforbox-likeorbits(dot-dashed lines).
chaoticwith∆ω&10−2.Wemaythusinferthatindeedthe bins have T .1, so we may expect the chaotic diffusion
dyn
existence of resonances slows down chaotic diffusion. to havean impact on theshape evolution.
Toquantifytheimpactthatchaoticdiffusionmayhave
However, it may also be considered somewhat differ-
ontheshapeofamodel,weusedthefollowingapproach.We
ently, as sketched in Fig. 3. Examination of the shapes of
took the same set of 104 orbits from a self-consistent solu-
low-energy chaotic orbits of the γ = 1 model reveals that
tionasabove,andre-integratedthemintwovariantsoffixed
they are mostly sticky, having a more-or-less well defined
(time-independent)potential:asmoothDehnenanalyticpo-
shapeforanyintervalof100T .Thismeansthatevenifan
dyn tential (same as was used in SM) and a frozen-N-body po-
orbitisscatteredtoadifferentregionofphasespace,itmay
tential obtained from a rigid distribution of 106 particles
stillretainadistinctshapeandsorepresentuseful‘building
representingthedensityoftheDehnenmodel.Wecalculate
blocks’.Itcouldbeexpectedthatthereisroughlyabalance
theorbitsfor10subsequentintervalsoftime,100timeunits
between orbits which became rounder and those which be-
each (not to be confused with 100T in the previous ex-
cameflatter.Ontheotherhand,orbitsfromthehigh-energy dyn
periment!).Oneachoftheseintervals,wesample100points
partoftheγ =1modelareinitiallycloser tohavingawell-
fromeachorbit,therebycreatinganN-bodymodelwith106
defined shape, but become more fuzzy-looking at the end.
particles,andrecordtheshapeofthismodel.Thedifference
They are not balanced by the opposite flux of orbits be-
from the previous experiment is that we evolved the orbits
coming more regular, since there were initially very few of
forfixed‘physical’time,ratherthanforfixednumberofdy-
thosenear-roundchaoticorbitsintheself-consistentequilib-
namical times which depend on energy. This way,theinner
rium.AsimilarargumentisusedinMuzzio et al.(2005)and
parts were much ‘older’ in terms of dynamical time.
in Aquilano et al. (2007) to explain why it is advantageous
to have a density model which is less triaxial in the outer TheresultsofthiscalculationarepresentedlaterinSec-
parts,andthereforecanadoptalargeinitialfractionoffully tion 5, for theγ =2case only.They showthat evolution of
chaotic(almost round)orbitstobalance this‘puffingup’of shape during 1000 time units is substantial, and compara-
initially more elongated chaotic orbits. Anyway, the exis- ble to theevolution exhibited in the self-consistent N-body
tence of shape-supporting resonances with their associated simulation. For the γ = 1 case, no evolution was observed,
families of orbits, even chaotic, seems to be an important as expected from theaboveconsideration.
condition for slowing down chaotic diffusion. There are obvious ways to suppress this chaotic dif-
fusion in SM. First of all, one may reduce the fraction of
On the other hand, orbits with ∆ω < 10−3 did not
chaotic orbits by assigning them a penalty in the objective
substantiallychangetheirshapeineithercase,neitherindi-
function.Thedifficultyliesinthepoorlydefineddistinction
vidually,nor on average.
between regular and sticky orbits, and a small fraction of
Uptonowweconsideredtheevolutioninshapeduring (weakly) chaotic orbits in the model is unavoidable.
a period of a fixed number of dynamical times, i.e. a time Second,wedeliberately chosearathersmall intervalof
whichvariesstronglywithenergy.Wemayinsteadfocuson integration (100T ) for SM, which is not enough to en-
dyn
thechangesthatmayoccurtothemodelduringafixedphys- surethatchaoticorbitsuniformlyfilltheirallowedregionof
icaltime(say,1000dimensionlesstimeunits,whichroughly phasespace.Indeed,amethodproposedbyPfenniger(1984)
corresponds to the Hubble time and to the time of the N- consists in an adaptive selection of the integration time for
body simulation which we will discuss in Sect. 5). For the orbitsthatseemtobechaotic,toimprovethecoverageofthe
weak-cusp case the dynamical time of orbits at the ener- available phase space: if the cell occupation numbers differ
gies where changes in shape occur is T & 50; hence the toomuchforthefirstandthesecondhalvesoftheorbit,we
dyn
chaotic diffusion occurring on timescales & 100T is es- continue the integration. The convergence, however, is not
dyn
sentially non-relevant on the timescale of the calculation. guaranteed, and can be rather slow, so one should impose
However, for the strong-cusp case the orbits in the inner anupperlimit forintegrationtime.Capuzzo-Dolcetta et al.
(cid:13)c 2011RAS,MNRAS000,1–14
Chaotic mixing and the evolution of triaxial galaxies 7
x-y y-z x-z γ=1 case
0.1
(D E/E)2 D E/E 0.3
(4,-1,-2) thin orbit
0.2
(initial)
g=2
x-y y-z x-z
0.01 0.1
6:7:8 resonance
g=1
(final)
0.05
x-y y-z x-z γ=2 case
bin
0.001
0.03
1 2 3 4 5 6 7 8 9 10
(initial)
Figure 4. Relative energy change, ∆E/E, after T = 1000, for
particles in 10 energy bins of N-body models (leftmost are the
x-y y-z x-z mostboundparticles).Theredlinescorrespondtoγ=2models:
the solid for N = 106 and the dashed for 2·106. The solid and
dashed blue lines correspond to a γ = 1 model, with N = 106
and5·106,respectively.Thebluedot-dashedlinecorrespondsto
(final) aN =106 modelbuiltbytheiterativemethod(Sec. 7).
Figure 3. An example of shape change for two 100Tdyn seg- body model (NM) by sampling points randomly from tra-
ments(separated by1000Tdyn)ofthesamechaoticorbit. jectories, their number being proportional to the weight of
Toptworows:orbitwithenergyfromthe2ndoutofthe10bins,in theorbitintheSM.OurstandardnumberofbodiesinNM
theweak-cuspmodel.Theinitialandfinalorbitsegmentsarepar- is106,whichmeansthat,onaverage,anorbitfromSM with
entedbydifferentresonancesandhaveratherwell-definedshapes. nonzero weight is sampled with 102 points.
Even though the change in shape is substantial, there may well ∼
We use the N-body code gyrfalcON (Dehnen 2000,
exist another orbit in the self-consistent solution which changes
2002) to evolve NM for T = 1000 N-body time units, cor-
intheoppositesense.
responding roughly to a Hubble time (1010 yr) for a model
Bottom two rows: orbit with an energy from the 9th bin in the
strong-cuspmodel.Theinitialshapeismoreelongatedalongthe scaledtoM =3·1011M⊙,a=5kpc.(Thehalf-mass(long-
xaxis,butthefinalonefillsalmostalltheequipotential surface, axis) radius of the γ = 1 model is 2.4 and the dynamical
whichisclosetoasphere.Therewerenoorveryfewsuchround time for this radius is 23 time units; for γ =2 these values
orbitswhichcouldbecomeelongatedtocounteractthisdiffusion. are 1 and 6.1, respectively). The smoothing length was set
to ǫ=0.01.
First we consider only one variant of the SM, namely
(2007) found that even setting the maximum time equal to with both the fraction of chaotic orbits and the velocity
2THubble andthencontinuingtheintegration oforbitsuntil anisotropy being unconstrained. Other possibilities will be
5THubble,somechangeoftheoverallshapeofdensitydistri- considered in section 6.
bution is observed. Allmodelswereinitiallyinvirialequilibrium(virialra-
Third, part of the problem with chaotic orbits lies tio 2T/W = 1 within a fraction of percent accuracy), and
in their stickiness, which may be reduced at the stage of the kinetic and potential energy remained the same in the
orbit integration by adding a weak noise (random force) course of evolution to within few 10−3(10−2) in the weak
to the equations of motion. Kandrupet al. (2000) and (strong) cusp case. To check the s×tability of our model, we
Siopis & Kandrup (2000) demonstrated that even a weak measured thechange in time of thefollowing parameters:
noise may dramatically increase the rate of chaotic diffu-
sion. mean-squarechangeofenergiesofparticles,toquantify
•
These efforts, however, are pretty much useless when theeffect of two-body relaxation;
wearetryingtosuppressthechaoticdiffusioninanN-body density profile (spherically averaged);
•
model, since we found it very difficult, if at all possible, to axialratiosofequidensitysurfacesdependingonradius;
•
preserve the attribute of chaoticity in transferring orbits to
a slightly different potential.
5.1 Energy change of individual particles
First we focus on the energy change of individual parti-
5 N-BODY EVOLUTION OF MODELS
cles, which can be due to two effects. One is the non-
In order to confirm our expectations about the chaotic dif- stationarity of the gross potential. This, in turn, may be
fusion, and to test the overall stability of triaxial Dehnen eitherduetothemodelbeinginitiallynotinperfectequilib-
models,weconvertourSchwarzschildmodel(SM)toanN- rium,orduetolarge-scaledynamicalinstabilities.Although
(cid:13)c 2011RAS,MNRAS000,1–14
8 E. Vasiliev, E. Athanassoula
0.7 0.7
y n
0.6 od 0.6 ze
b o
0.5 S, ny 0.5 S, fry
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
S, smooth S, smooth
y y
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
1 1
y n
d e
o z
10-1 , nb 10-1 , fro
10-2 Dw 10-2 Dw
10-3 10-3
Dw , smooth Dw , smooth
10-4 10-4
10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 1 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 1
Figure5. CorrespondencebetweenorbitsinSM andinNM,for Figure 6. Correspondence between orbitsinSM (smooth) and
asampleof2000orbitsintheinner40%oftheweak-cuspmodel, frozen-N-bodypotential,forthesamesampleasinFig.5.
whichhadcompletedatleast50periodsduringtherun. The top panel displays a correlation between orbit shapes, and
Top panel: correlation between orbit shapes (more specifically, shows that orbits in the frozen-N-body potential resemble their
flattening along the y axis); the shapes of the orbits in the N- counterparts inthesmoothpotential quitewell.
bodyrunarereasonablysimilartothoseoftheparentSM orbits. Thebottompaneldisplaysasimilarcorrelation,butnowbetween
Bottompanel:correlationbetweenFDR–meanttomeasurethe the FDR. Orbits that were regular in the smooth potential (i.e.
degreeofchaos–forthetwosetsoforbits.Notethedifferencein had∆ω<10−3),havesystematicallyalowerFDRinthefrozen
scalebetweentheordinateandtheabscissa.InNM thechangeof potential, although this rarely gets below 10−3. This is due to
orbit frequencies is caused mostly by fluctuations in energy (see thegraininessofthepotential andnottoanenergyerror,which
Fig.7),andsohasessentiallynocorrelationwitheitherthetrue isstillanorderofmagnitudelower.
degreeofchaos,orthe∆ω oftheparentorbitinSM.Almostall
orbitsinNM have ∆ω>10−3,i.e.above the thresholdusedfor
separatingregularfromchaotic orbitsinthesmoothpotential. for diffusion, apart from a possible initial period of faster
growthifthemodelwasnotinperfectequilibrium).Figure4
showsthat thischangeremains quitelow fortheweak-cusp
theSchwarzschildmodellingtechniqueisaimedatconstruct- model, of order ∆E/E few percent, weakly depending
| | ∼
ing models in equilibrium, the actual accuracy of this equi- ontheinitialenergy.Inthestrong-cuspcase,therelaxation
libriummayvary.Inparticular,wefoundthatthestandard rate is 5 times stronger, although the simulation time is
∼
practiceofassigninginitialconditionsfororbitsfrom agrid still much shorter than the relaxation time for all energies.
of discrete energy levels results in quite large initial fluctu- This, however, leads to an important implication for
ations of energies. We therefore adopted a smooth energy the detection of chaos by the frequency diffusion rate. In
distribution for orbits in the SM. Nevertheless, the change fact, for a perfectly regular orbit ∆ω ∆E/E, and so the
≃
of energy that is due to the fact that the model is not in lowerboundonFDRdependsontheenergyconservationof
perfectequilibriuminitiallyshouldbeconfinedtotheinitial agiven orbit and thenumbersshown in Fig. 4 show that it
times of the simulation and thus should not influence the is indeed quite high. In the initial distribution of FDR (in
following discussion. thesmoothpotentialofSM)therewerenomorethanafew
Thesecondreasonforanenergychangeistheunavoid- percent of orbits with FDR ∆ω > 0.03, while in the NM
able two-body relaxation. It leads to a random-walk in en- they comprise the majority of orbits. It therefore does not
ergy space, with a mean squared deviation of energy ∆E2 makesensetouse∆ωasanindicatorofchaosfororbitsina
growing linearly with time. To estimate the importance of liveN-bodysimulation.Toachieveanacceptablylowenergy
this effect, we measured the squared relative change in en- diffusionrateof∆E/E 10−3,wewouldneedtohaveofthe
ergy, (∆E/E)2, accumulated by the end of the integration orderof109 or1010 part≃icles(basedontheratesforsquared
time and averaged over particles in several energy bins (we energychangefromFig.4beinginverselyproportionaltothe
checked that it indeed grew linearly with time, as expected numberof particles).
(cid:13)c 2011RAS,MNRAS000,1–14
Chaotic mixing and the evolution of triaxial galaxies 9
Not surprisingly, there is essentially no correlation be-
tween theFDR in thesmooth-potential and in the N-body
bits models(Fig.5).Thiscastsashadowontheusefulnessofat-
or temptstocontroltheamountofchaosintheNM bychang-
of
n ing the properties of the SM. The tests in the next section
o
cti confirm that indeed the evolution of the NM is basically
a
fr independent of whether most orbits in SM are regular or
irregular. On the other hand, orbit shapes in the NM are
mostlyclose tothoseoftheirparentorbitsfrom theSM,as
can be seen from the top panel of Fig. 5. Furthermore, we
found no correlation between the change of orbit shape in
the NM (which may be taken as an indicator of its chaotic
10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1 1 behaviour)andits∆ω eitherintheNM,orintheSM (not
Dw or D E/E shown here).
It is interesting to compare the orbits in the live N-
Figure 7. Histograms of FDR (∆ω, solid lines) and of energy body models with the corresponding orbits in the frozen
conservation error (∆E/E, dashed lines) for orbits from Figs. 5
N-body potential (Fig. 6). The latter share with the for-
and 6. The blue linecorresponds to orbits inthe smooth poten-
mer the energy conservation error caused by the tree-code
tial, greenone to orbits inthe frozen-N-bodypotential, andthe
force approximation and by finite-difference integration er-
redlinetoorbitsinthelivesimulation.
rors,butthiscan bekeptaslowas necessary bychoosing a
Forthelivesimulationtheenergyconservationerrorisquitelarge,
between 10−2 and 10−1, and it determines the error in the fre- suitable tree cell opening angle θ and integration timestep.
quency estimation. For the frozen-N-body potential the energy Forourrunsthisenergyerrorisbetween10−4 and10−3 for
conservationisstillfarfrombeingperfect(oforder∼10−4..10−3, all orbits, i.e. much lower than the energy diffusion in the
due to the approximations of the tree-code and the finite num- live simulation which, as shown in Fig. 7, is of the order of
berofparticles),butgoodenoughnottobethemajorsourceof 10−2 or 10−1. Yet ∆ω is an order of magnitude larger than
frequencydiffusion. the energy error and we checked that it does not change
noticeablywhenweincreasedenergyconservationaccuracy.
This shows that the main contribution to ∆ω comes from
graininessofthepotential,andformostorbitsisnotcaused
log(r /r )
exact bythenon-integrability of corresponding smooth potential.
0.4
Note in particular that for the case of a spherically sym-
metric frozen-N-body potential, ∆ω is also between 10−3
0.2 and 10−2 for the majority of orbits, which confirms that
thislowerlimitiscausedbygraininess, notby‘real’chaotic
propertiesofthepotential).Similarvaluesof∆ωwerefound
0 inValluri et al.(2010)fororbitsinsphericalNFWpotential
represented by 106 frozen particles (theirfig. 1).
-0.2
5.2 Evolution of global quantities
-0.4 For all models considered, the velocity distribution (not
log(r)
shown) and the density profile (Fig. 8) do not show sub-
-2 -1.5 -1 -0.5 0 0.5 1 1.5
stantial evolution, except for an unavoidable smoothing of
thecuspatr.ǫ.Themostimportantchangesoccurinthe
Figure 8. Density profile evolution for the γ = 2 model. We axial ratios of themodel (Fig. 9).
plot the difference between density in the model and the exact There are two different methods to measure the
value (log(ρ/ρexact)) as a function of radius (log(r)). The red axial ratio and its dependence on radius, and they
solidcurvecorresponds totheinitialmodel,whichisfairlyclose
give very similar results in our case. The first method
to the exact Dehnen profile and the green dashed curve corre-
(Athanassoula & Misiriotis 2002)istosortparticlesinden-
sponds to the model after it evolved over 10 time units. From
sityandbin themintoN bins,thencalculatetheprincipal
thelatterweseethatthesmoothingreducesthedensityatr.ǫ b
axesoftheinertiatensorofallparticlesbelongingtoagiven
and distributes the excess of mass slightly further from centre.
Thebluedottedlinecorrespondstotheendofsimulation,where bin(ifthedensityisasmoothmonotonicfunctionofradius,
the axial ratios have changed substantially. For comparison, the thesebinsareroughlyellipsoidalshells).Thesecondmethod
browndot-dashedlinegivestheinitialdensityinSM withinitial – a variation of that used by Dubinski& Carlberg (1991) –
conditionsassignedonaregulargridat20fixedenergylevels(see consistsin binningparticles intolayersboundedbyconcen-
Sec. 6) and clearly deviates substantially from the exact profile tric ellipsoids whose axes are determined by the following
atsmallandlargeradii. iterative procedure6. Consider first the inner shell contain-
ing 1/N of the total mass. We search for an ellipsoid with
b
axesa ,b ,c suchthatthemomentofinertiaofallparticles
0 0 0
6 Amoreextended discussionisfoundinZempetal.(2011).
(cid:13)c 2011RAS,MNRAS000,1–14
10 E. Vasiliev, E. Athanassoula
1 1
0.9 0.9 b/a
b/a
0.8 0.8
0.7 0.7
0.6 0.6 c/a
c/a
0.5 0.5
time time
0.4 0.4
0 200 400 600 800 1000 0 200 400 600 800 1000
1 1
0.9 0.9 b/a
b/a
0.8 0.8
0.7 0.7
0.6 0.6 c/a
c/a
0.5 0.5
time time
0.4 0.4
0 200 400 600 800 1000 0 200 400 600 800 1000
Figure 9. Evolution of the b/a and c/a axial ratios, starting Figure 10. Evolution of the b/a and c/a axial ratios, starting
from0.79and0.5,respectively. from0.79and0.5,respectively,forvariantsoftheγ=2model.
The top panel compares two γ = 1 models. The red solid line ThetoppanelshowsNp=5·105 particlesNM createdfromSM
correspondstoaN =106 referencerun,thedottedbluelinetoa having Ns = 20 shells and No = 104 orbits integrated for 100
N =5·106 high-resolutionrunandthegreendash-dottedlineto orbital periods. The green solid lines correspond to the uncon-
aN =106 modelbuiltwiththeiterativemethod. Theyallshow strainedmodel,thebluedashed linetothemodel withaprefer-
little evolution, and the high-resolution run conserves its shape enceofregularorbits,thereddottedlinetothemodelwithapref-
almostperfectly. erenceofchaotic orbits,andthebrowndot-dashed tothemodel
The bottom panel shows that the γ = 2, Np = 106 simulation createdfromgridinitialconditions (theother threemodelshave
(redsolidline)showsasubstantialevolutionofshape.Itiscom- random IC). This last one clearly displays the strongest shape
pared to the evolutions of orbits in the fixed potential (blue – evolution, while there is almost no difference between ‘mostly
smoothDehnen,green–frozenN-bodypotential),whichdemon- regular’,‘mostlychaotic’and’unconstrained’ models.
stratesimilaramountofshapechange exclusivelyduetochaotic ThebottompanelshowstheeffectofparticlenumberinNM and
diffusion. ofnumberoforbitsandshellsinSM.Thegreensolidlineshows
the shape evolution of a model with Ns = 20,No = 2500 and
Np=5·105integratedfor100dynamicaltimes(sameasinthetop
withinthisellipsoidhasthesameaxisratioasthisbounding panel). The brown dot-dashed shows the Ns = 50,No = 15000
ellipsoid. We start from a spherical shell and find the mo- andNp=5·105 modelintegratedfor100dynamicaltimes.The
ment of inertia of particles within this radius, then change red dashed line corresponds to the Ns = 50,No = 12500 and
Np = 106, integrated for 500 dynamical times (same as the red
theaxial ratioof theboundingellipsoid tothesevaluesand
solid curve in the bottom panel of Fig. 9). The blue dotted line
repeat the iteration until convergence is achieved. We then
corresponds to a model similar to that of the red line, but with
remove the inner particles from consideration and repeat Np=2·106.
the procedurefor the next shell. In practice, the numberof Thisfiguredemonstrates thatincreasingthenumberofparticles
shellsshouldberathersmalltoavoid‘shell-crossing’andthe does slow down the shape evolution, and that at fixed Np aSM
divergence of iterations. withmoreshellsandorbitsbehaves better.
In what follows, we usually present the axis ratios for
the2ndoutof4bins(particlesbetween25and50%),which
contains particles close enough to the centre yet not much
affectedbysoftening.(Indeed,fortheγ =2modeltheaxial
ratio changed with almost thesame rate at all radii).
(cid:13)c 2011RAS,MNRAS000,1–14