Dynamical friction in flattened systems: A numerical test of Binney's approach

0.7 MB·English

by  Jorge Penarrubia

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed2February2008 (MNLATEXstylefilev1.4)
Dynamical friction in flattened systems: A numerical test
of Binney’s approach
Jorge Pen˜arrubia1 , Andreas Just1, Pavel Kroupa2,3
1Astronomisches Rechen-Institut, M¨onchhofstrasse 12-14, D-69120 Heidelberg, Germany
2Institut fu¨r Theoretische Physik und Astrophysik, Universit¨at Kiel, D-24098 Kiel, Germany
3Heisenberg Fellow
4
2February2008
0
0
2
ABSTRACT
n
Wecarryoutasetofself-consistentN-bodycalculationstoinvestigatehowimportant
a
thevelocityanisotropyinnon-sphericaldark-matterhalosisfordynamicalfriction.For
J
this purpose we allow satellite galaxies to orbit within flattened and live dark-matter
9
haloes (DMHs) and compare the resulting orbit evolution with a semi-analytic code.

This code solves the equation of motion of the same satellite orbits with mass loss
1
v and assumes the same DMH, but either employs Chandrasekhar’s dynamical friction
9 formula,whichdoesnotincorporatethevelocityanisotropy,orBinney’sdescriptionof
5 dynamicalfriction in anisotropicsystems. In the numericaland the two semi-analytic
1 models the satellites are given different initial orbital inclinations and orbital eccen-
1 tricities,whereasthe parentgalaxyiscomposedofaDMH withaspectratioqh =0.6.
0 We find that Binney’s approach successfully describes the overall satellite de-
4 cay and orbital inclination decrease for the whole set of orbits, with an averaged
0
discrepancy of less than 4 per cent in orbitalradius during the first 3 orbits. If Chan-
/
h drasekhar’s expression is used instead, the discrepancy increases to 20 per cent. Bin-
p ney’s treatment therefore appears to provide a significantly improved treatment of
- dynamical friction in anisotropic systems.
o
The velocity anisotropy of the DMH velocity distribution function leads to a
r
t significant decrease with time of the inclination of non-polar satellite orbits. But at
s
a the same time it reduces the difference in decay times between polar and coplanar
: orbits evident in a flattened DMH when the anisotropic DMH velocity distribution
v function is not taken into account explicitly. Our N−body calculations furthermore
i
X indicate that polar orbits survive about 1.6 times longer than coplanar orbits and
thattheorbitaleccentricityeremainsclosetoitsinitialvalueifsatellitesdecayslowly
r
a towardsthegalaxycentre.However,orbitsofrapidlydecayingsatellitesmodelledwith
the semi-analytic code showa strongorbitalcircularisation (e˙ <0)not presentin the
N-body computations.
Key words: stellar dynamics – methods: N-body simulations– methods: analytical
– galaxies: kinematics and dynamics – galaxies: haloes – galaxies: dwarf
1 INTRODUCTION Matter model (CDM). In this framework, aspherical bound
darkmatterhaloes(DMHs)formasaresultofgravitational
clustering. The inclusion of gas dynamics in the CDM sim-
According to current ideas galaxy formation began with
ulations (Udry & Martinet 1994, Dubinsky 1994) results
small-amplitude Gaussian fluctuations at the early stages
to a Gaussian distribution of DMH density aspect ratios,
of the Universe. In hierarchical cosmological models, these
q c/a>0, where c and a are the minor and major axes
fluctuationsdecreasewithincreasingscales,resultingfirstin h ≡
of an oblate spheroid, of mean <q >=1/2 and dispersion
theformation oflow-massobjectsthatmaymerge,building h
equal to 0.15 (Dubinsky 1994). The degree of asphericity
up ever more massive structures. The shape and morphol-
dependsonthedarkmatternature.Forexample,numerical
ogy of these objects are strongly dependent on the cosmo-
computations based on the ΛCDM predict halo axis-ratios
logical models, as one can conclude from N-body computa-
ofq 0.7(Bullock2001),HotDarkMattermodelsleadto
tions, although none of them predict spherical structures. h ≃
haloes as round as q =0.8 (Peebles 1993), whereas candi-
The most successful hierarchical theory is the Cold Dark h
(cid:13)c 0000RAS
2 J. Pen˜arrubia, A. Just & P. Kroupa
datessuchascoldmoleculargas(Pfenniger,Combes&Mar- andWhade&Donner(1996)findthatdynamicalfriction is
tinet 1994) and massive decaying neutrinos (Sciama 1990) artificially increased due to numerical noise if the particle
may produce halo profiles as flattened as q =0.2. numberissmall. Semi-analyticmethodsthatincludeChan-
h
Observationally, measuring galaxy axis-ratio is compli- drasekhar’s dynamical friction have been demonstrated to
cated: (i) Stellar kinematics: Olling & Merrifield (2000) ob- reproduceaccuratelytheoverallevolutionofsatellitegalax-
tain an axis-ratio of q 0.8 for our Galaxy. This method ies(e.gVel´azquez&White1999,Taylor&Babul2001)and,
h
≈
hasthedisadvantage of havingaccess to information of our therefore, represent a useful tool in order to carry out ex-
Galaxyonlyonsmallscales.(ii)Theflyinggaslayermethod tensivestudies with a large parameter space.
(Olling 1996, Becquaert, Combes & Viallefond 1997): as- It is, however, still unclear how the inhomogeneity of the
suming that the HI emission comes from gas in hydro- system distribution affects the satellite dynamics. Whereas
staticequilibriumingalacticpotentialsaxis-ratiosaslowas DelPopolo(2003),DelPopolo&Gambera(1998)andMaoz
q 0.3 are obtained for the galaxies NGC 891 and 4244, (1993)showthatdynamicalfrictionincreasesthesteeperthe
h
≈
(iii) Warping gas layer: Hofner & Sparke (1994) find axis- density profile is, Just & Pen˜arrubia (2002), following the
ratios of approximately 0.7 for NGC 2903 and of q 0.9 theoretical scheme proposed by Binney (1977), hereinafter
h
≈
for NGC 2841, 3198, 4565 and 4013, (iv) X-ray isophotes: B77, finda negligible effect on thesatellite orbit dueto the
Boute & Canizares (1998) measure values of q 0.5 for symmetryoftheinhomogeneoustermsofdynamicalfriction.
h
≈
NGC3923,1332and720,(v)Polarringgalaxies:Arnaboldi Although Chandrasekhar’s formula, which assumes an
et al. (1993) and Sackett et al. (1994) find an axis-ratio isotropic velocity distribution, reproduces accurately dy-
of q 0.3 for NGC 4650A, 0.5 for the galaxy A0136- namical friction in spherical systems, an analytical study
h
≈
0801 and 0.6 for AM2020-504, (vi) Precessing dusty discs: ofStatler(1991)showsthatinthecaseofSt¨ackelpotentials
Steinman-Cameron, Kormendy & Durisen (1992) measure the velocity anisotropy produce strong effects on the satel-
an axis-ratio of 0.9 for thegalaxy NGC4753. lite orbit. The N-body computations of PKB confirm that
Another method, which we focus on here, is the anal- Chandrasekhar’s formula cannot account for the resulting
ysis of satellite dynamics. There are two main different ap- satellitedecayandevolutionifthehaloisflattened.Theaim
proaches to infer the halo shape from satellites. First, one ofthispaperistoimplementasemi-analyticschemecapable
may attempt to reproduce the observed tidal streams of of trackingthedynamicalevolution of substructureswithin
Milky Way satellites as done, for instance, by Ibata et al. flattened as well as spherical DMH’s. With this purpose in
(2001)andLawetal.(2003)who,frommeasurementsofve- mind,weimplement in ourcodetheanalytic expressionsof
locity,positionandstructureoftheSagittariusdwarfgalaxy, dynamical friction in systems with anisotropic velocity dis-
constraintheinitialparameterspaceand,subsequently,cal- persionssuggestedbyBinney(B77)andalsousedbyStatler,
culate in detail the satellite mass loss. They find that the whichreproducesChandrasekhar’sfornullanisotropyinve-
Milky Way halo potential cannot be more flattened than locity space. Weas well compare theresults of usingChan-
0.8, otherwise tidal streams would betoo spread out and drasekhar’s formula in axi-symmetric systems to determine
≈
thick compared to the observations due to orbital preces- theeffects of the velocity anisotropy on satellite decay.
sion. The second approach is a statistical study of satellite Section2introducesthemodels.InSection3weprovide
distributionaroundspiralgalaxies. Holmberg(1969),Zarit- the code and galaxy parameters. We outline the dynamical
sky & Gonz´alez (1999), Prada et al. (2003) and Sales & frictionapproachesinSection4whereasinSection6wepro-
Lambas (2003) point out that satellites around disc galax- pose a simple techniqueto fit a free parameter (in our case
ies are found more often aligned with the poles of the host theCoulomblogarithm)totheN-bodydata.InSection7we
galaxy, the so-called ’Holmberg effect’, whereas Quinn & study how flattened DMHs affect satellite decay and calcu-
Goodmann (1986) findin their N-body study that thedisc late the degree of accuracy of Binney’s formula. The paper
alonecannotaccountfortheoriginal statisticaldistribution concludes with Section 8.
of Holmberg’s data. A remedy may be sought in the form
ofan extendednon-sphericalDMH.Ananisotropic velocity
(andmass)distributionwillcauseasatellite’s orbittoalign
with the axes of the velocity ellipsoid of the host galaxy
2 GALAXY AND SATELLITE PARAMETERS
(Pen˜arrubia, Kroupa & Boily 2001, hereinafter PKB). In
both schemes, a large number of numerical calculations is OurDMHmodelisthatusedbyPKBtofacilitate aninter-
needed.Intheformerapproachoneshouldintegrateseveral comparison of the disc and bulge effects on satellite decay
initialorbitalparameterstofindthebestfittotheobserved (Pen˜arrubia 2003). Here we do not add the bulge and disc
satellite characteristics, whereas in the latter approach an componentsin ordertodistill thedependencyofdynamical
initialsatellitesampleshouldstatisticallyreproducethedis- friction on thevelocity anisotropy in a spheroidal DMH.
tributions expected from cosmological models of DMH for- In order to minimise computational time when con-
mation.Sofarthisisprohibitivelytime-expensivebymeans structing flattened DMHs, we apply a highly-efficient tech-
of any of thepresent N-bodyalgorithms. nique using multi-pole potential expansions to tailor the
Several studies of satellite decay have shown that, local velocity ellipsoid to the required morphology (Boily,
in spherical systems, Chandrasekhar’s dynamical friction Kroupa&Pen˜arrubia2001). TheMaGaliecodescaleslin-
(Chandrasekhar 1943) is accurate enough if the Coulomb earlywithparticlenumberandhencewecanconstructflat-
logarithm remains as a free parameter to fit to the N-body tenedDMHsconsistingof >106particlesormore,inashort
data(e.g.vandenBoschetal.1999,Colpietal.1999)since computational time. ∼
it also depends on the code parameters and the number of Following PKB, the flattened DMH is described by a
particlesemployed.Forinstance,Prugniel&Combes(1992) non-singular isothermal profile which, in cylindrical coordi-
(cid:13)c 0000RAS,MNRAS000,000–000
Dynamical friction in flattened systems: A numerical test of Binney’s approach 3
nates (for m2[u=0] see eq. 7), can bedescribed as
Symbol Value(ph.u) Value(m.u)
ρh[m2(0)]= 2πM3/h2rαcutexp(cid:2)m−2(m0)2(+0)γ/2rc2ut(cid:3), (1) (DHM2)H Nh 1400000
with m2(0)=R2+z2/q2 Mh 7.84×1011M⊙ 14.00
h
qh 0.60 0.60
Mh being the DMH mass, rcut the cut-off radius and γ the γ 3.5kpc 1.00
core radius, and rcut 84.00kpc 24.00
α 1 √πβexp(β2)[1 erf(β)] −1 = (2) Satellite Ns 40000
≡{ − − } (S1)
1+√πβ+(π−2)β2+O(β3) Ms 5.60×109M⊙ 0.10
where β = γ/rcut<1/24 in our calculations. For β = 1/24 Ψ(0)/σ02 5.00 5.00
we find α 1.076∼ 1 already and hence thereafter we set rc 0.67kpc 0.19
≃ → rt 7.24kpc 2.07
α=1 in our analysis.
c 1.03 1.03
Weuseself-consistent Kingmodels(King1966) torep-
<r> 1.64kpc 0.47
resentourdwarfgalaxies.Thesemodelsfitearly-typedwarf σ0 60.30kms−1 0.23
galaxies (Binggeli et al. 1984). For a comparison with the
workofPKBweadoptc=log (r /r )=0.8,wherer and
10 t c c
r are the core and tidal radii, respectively.
t Table1.Primarygalaxyandsatellitemodels.TheDMHhasan
To construct the models we choose the satellite mass aspect ratio qh = 0.6. For the satellite model: Ψ(0) = Φ(rt)−
Ms and rt. The tidal radius is determined by computing Φ(0), Φ(0) are the central potential and Φ(rt) the potential at
the density contrast, ρs(rt)/ρg(ra) 3, at the apo-centric the tidal radius (following the notation of Binney & Tremaine
∼
distance (ra = 55 kpc) at t = 0, ρg(r) being the averaged 1987); σ0 is the velocity dispersion at the centre, and < r >
densityofthegalaxy(sameprocedureasVel´azquez&White theaverageradiusofthesatellite.Theunitsaresuchthat Ph.u.
1999). This guarantees that all satellite particles are bound means’physicalunits’andm.u.’modelunits’.NhandNsarethe
numberofparticlesusedtorepresenttheDMHandthesatellite,
at t=0.
respectively.
We employ the system of units of PBK, which refers
to the parameters of the Milky Way disc. Adopting M =
d
R = G = 1 and according to Bahcall, Smith & Soneira
d Myrwhichisabout1/40ththedynamicaltimeofoursatel-
(v1e9lo8c2i)t,yMundi=ts5a.r6e,×r1es0p10eMcti⊙vealyn,d1R.3d =1307.5ykrpacn,dth2e6t2imkmesa−n1d. lite.Wehavethreeresolutionzones,eachwith643grid-cells:
× (i) The inner grid covers out to 3 radial halo scale-lengths,
The values of the galaxy and the satellite parameters
providingaresolutionof350pcpergrid-cell.(ii)Themiddle
can be found in Table 1. The parent galaxy corresponds to
grid covers the whole galaxy, with an extension of 24 halo
the model G2 of PKB, with the difference that we remove
scale-lengths (84 kpc), giving a resolution of 2.8 kpc per
thedisc and bulgecomponents here.
grid-cell. The satellite always orbits within this grid except
attheverylatestagesofitsevolution,avoidingcross-border
effects (see Fellhauer et al. 2000) and (iii) The outermost
gridextendsto348kpcandcontainsthelocaluniverse,ata
3 NUMERICAL CALCULATIONS
resolutionof11.6kpc.Asforthesatellitegrid-structure,the
3.1 Code parameters resolutions are 816 pc per grid-cell for the inner grid that
extends to 24.48 kpc, 1.2 kpc per grid-cell for the middle
The numerical experiments were carried out by using the
grid which extendsto36 kpc,and 11.6 kpcpergrid-cell for
code Superbox which is a highly efficient particle mesh-
theoutermost grid that covers thelocal universe.
algorithm based on a leap-frog scheme (for a detailed de-
The selection of grid parameters ensures the conserva-
scription see Fellhauer et al. 2000). Superbox has already
tion of energy and angular momentum for satellites in iso-
been implemented in extensive studies of satellite disrup-
lation over times as long as our calculations to better than
tionbyKroupa(1997),Klessen&Kroupa(1998)andPKB.
1% for all themodels.
Theprogram calculatestheaccelerations usingahighorder
NGP (‘nearest grid point’) force calculation scheme based
on the second derivatives of the potential. A self-consistent
3.2 Orbital parameters
system of several galaxies can be treated by forming sub-
grids (3D ’boxes’) which follow the motion of each galaxy. Theparameterspaceofsatellitegalaxies isextremelylarge.
Each sub-grid has three levels of resolution, the two finest A complete survey should account for different satellite
levels co-move with the galaxies allowing a high-resolution masses, apo-galacticon distances, orbital eccentricities and
calculationoftheforcesactingontheparticles,whereasthe inclinations...etc. In this paper, we carry out a set of cal-
third one covers the local ’Universe’. The finest levels are culations selecting those parameters that best reflect the
centred on thedensity maximum of thegalaxy, which is re- effects of the velocity dispersion anisotropy on satellite de-
computed at every time-step. cay. These parameters are: (i) the initial orbital inclination
ThecodeparametersarethoseofPKB.Inthatpapera (i),definedastheanglebetweentheinitialangularmomen-
detailed description of the system and the grid structure is tum vector of a satellite and the axis perpendicular to the
presented, whereas here we merely give a brief description axi-symmetry plane (selected as the z-axis). We expect or-
oftheN-bodyparameters.Ourintegrationtimestepis0.39 bitalinclination todecreaseintimeaspredictedbyBinney.
(cid:13)c 0000RAS,MNRAS000,000–000
4 J. Pen˜arrubia, A. Just & P. Kroupa
where(v ,v )arethecoordinatesofthesatellitevelocityin
R z
Name Gal. Sat. i e rp ra thisframe.
model model [kpc] [kpc]
As Binney shows, a body with mass M will suffer a
s
H2S100 H2 S1 0◦ 0.5 18 55 decreaseofitsorbitalinclinationwheneverBz >BR (oblate
H2S130 H2 S1 30◦ 0.5 18 55 halo).Iftheorbitiseithercoplanarorpolar,theinclination
H2S145 H2 S1 45◦ 0.5 18 55 remains constant since, respectively, the perpendicular and
H2S160 H2 S1 60◦ 0.5 18 55 the planar component of v is zero. It is straight-forward to
H2S190 H2 S1 90◦ 0.5 18 55
show that this expression reproduces Chandrasekhar’s for
H2S100c H2 S1 0◦ 0.3 30 55 ev =0,
HH22SS114350cc HH22 SS11 4350◦◦ 00..33 3300 5555 Fch =−4πGMsρh[m2(0)]lnΛ erf(X)− √2Xπe−X2 vv3s, (4)
H2S160c H2 S1 60◦ 0.3 30 55 (cid:2) (cid:3) s
H2S190c H2 S1 90◦ 0.3 30 55 where X = v /√2σ.
s
| |
H2S100e H2 S1 0◦ 0.7 10 55 One important aspect to note is that both expressions
H2S130e H2 S1 30◦ 0.7 10 55 of dynamical friction include here an anisotropic halo den-
H2S145e H2 S1 45◦ 0.7 10 55 sity, denoted by ρh[m2(0)] = ρh[R2+z2/qh2] in cylindrical
H2S160e H2 S1 60◦ 0.7 10 55 coordinates. This is the “local approximation” made here
H2S190e H2 S1 90◦ 0.7 10 55 butwenotethat thederivation of eq.4assumes auniform,
infinitelyextendedbackgroundmedium.Inpractice,thelo-
Table2.Thenumericalexperiments.Theperi-andapo-galactica cal approximation implies the only difference between both
are rp and ra, respectively, and e = (ra−rp)/(ra +rp) is the expressions to be the anisotropic terms of the velocity dis-
orbitaleccentricity. tribution.
We note that all the calculations proceed with the same
5 THE SEMI-ANALYTICAL CODE
orbital sense, but this is not important since the halo is
non-rotating. (ii) The satellite’s initial orbital eccentricity,
Inordertoanalyse theaccuracy oftheanalyticexpressions
definedas e=(r r )/(r +r ),wherer ,r aretheapo-
a p a p a p of dynamical friction we have constructed a semi-analytic
−
and peri-galacticon, respectively.
algorithmthatsolvestheequationofmotionofapoint-mass
Theparametersofthenumericalexperimentsarelisted
satellite
in Table 2.
d2x
dt2 =Fg+Fdf, (5)
where F is the specific force from the parent galaxy and
4 HALO DYNAMICAL FRICTION g
Fdf that dueto dynamical friction (eqs. 3 or 4).
Chandrasekhar’sexpressioncannotexplainsomeeffectsob- If the parent galaxy follows the fixed density profile of
served in N-body calculations of satellite decay within flat- eq. (1), the specific force can be written in Cartesian coor-
tened haloes (PKB). Our aim is to check Binney’s approx- dinates (Chandrasekhar 1960) as
imation (B77) for systems with anisotropic velocity disper-
∞
du
sion. F = 2πGx ρ [m2(u)], (6)
g,i − iZ (1+u)2(1+e2 +u)1/2 h
For simplicity, we reproduce here the analytic for- 0 h
mulæ employed throughout this study (for a detailed anal- ∞ du
F = 2πGz ρ [m2(u)],
ysis of the friction force see Pen˜arrubia 2003 and Just & g,z − Z (1+u)(1+e2 +u)3/2 h
0 h
Pen˜arrubia 2002).
If the distribution function in velocity space is axi- where xi =x,y, the ellipticity of the galaxy is e2h ≡1−qh2,
symmetric, thezeroth order specific friction force is (B77) and
R2 z2
Fi =−2√2πρh[m2(0)σ]RG2σ2zMs√1−e2vlnΛBRvi, (3) m2(u)= 1+u + 1−e2h+u. (7)
2√2πρ [m2(0)]G2M √1 e2lnΛ The algorithm employed to solve eq. (5) is based on
Fz =− h σR2σz s − v Bzvz, tPhreesBsueltirsacl.h-1S9t8o6e)r,mwehtihchodp(rfoovridaescohmigphl-eatcecudreastceripsotilountiosnees
where i = x,y and (σR,σz) is the velocity dispersion el- withminimalcomputationaleffort.Thismethodisbasedon
lipsoid of aSchwarzschild distribution in cylindrical coordi- an adaptive step-size scheme, thus, being ideal for systems
nateswithconstantellipticitye2v =1−(σz/σR)2,lnΛisthe withnon-smoothpotentials,asmaybethecaseforsatellites
Coulomb logarithm of thehalo and on highly eccentric orbits when disc and bulge components
B = ∞dqexp(−vR21/+2qσR2 − 1v−z2/e22vσ+R2q), are inFcolrudcaedlc.ulating Fdf (eqs. 3 and 4) the Coulomb log-
R Z (1+q)2(1 e2+q)1/2 arithm is treated as a free parameter to be fitted to the
0 − v numerical orbits. The satellite mass M (t) is obtained from
s
∞ exp( vR2/2σR2 vz2/2σR2 ) theN-bodydata(seeSection7.1)andistreatedasaninput
B = dq − 1+q − 1−e2v+q , function. Numerical tests of orbits in a Keplerian potential
z Z (1+q)(1 e2+q)3/2 show that energy and angular momentum are conserved to
0 − v
(cid:13)c 0000RAS,MNRAS000,000–000
Dynamical friction in flattened systems: A numerical test of Binney’s approach 5
10−8 and10−9 perorbit,respectively,afterchoosingtheac-
0.08
curacy of the Bulirsch-Stoer algorithm to be EPS= 10−5.
This value remains fixed for the calculations presented in 0.06
this paper. An extensive description of our semi-analytic
0.04
code can befound in Pen˜arrubia (2003).
0.02 H2S160
0
0.08
6 DETERMINING THE COULOMB
LOGARITHM 0.06
TheCoulomblogarithmisusuallyfittonumericaldatawith 0.04
theaimofreproducingtheoverallorbitalevolutionbymeans
0.02 H2S145
ofsemi-analyticalgorithms(e.g.Fellhaueretal.2000).This
procedure may actually be considered as a “calibration” of 0
the semi-analytic code, which must be done carefully if a 0.08
detailed inter-comparison between different schemes of dy-
0.06
namical friction is desired. For that reason we present in
whatfollowsamethodtodescribetheaccuracyofthesemi- 0.04
analytic scheme.
0.02 e=0.3 H2S130 e=0.5 e=0.7
Onepossible way toquantifysimilarity of twoorbits is
via thequantity 0
1 2 3 1 2 3 1 2 3
2k
1
χ2 = 2k (ri−ri,n)2+σ2(r0)(ti−ti,n)2 , (8) Figure 1. The parameter χ for different orbital eccentricities
Xi=1(cid:2) (cid:3) and inclinations. Dotted lines denote the first 4 orbits, whereas
solidlinesthefirst3orbits,inbothcasesusingBinney’sexpres-
r being the satellite position vector at the peri- and apo-
sions for dynamical friction.Dashed and dot-dashed lines repre-
galactica and t the time at which the satellite passes by
sent the results using Chandrasekhar’s expression for k =4 and
these points. The subindex n denotes the numerical values 3 orbits, respectively. The χ values are normalised to the initial
and σ(r0) the velocity dispersion of the DMH at the initial apo-galacticondistancer0=55kpc.
galacto-centre distance. The sum is over a given number of
orbits k.
The definition of χ measures the divergence of the nu-
mericalandsemi-analyticalsatelliteposition.Thetermσ∆t (lnΛ [2.3,2.5]), which shows that this scheme provides a
allows for thepossibility that theorbital periods differdur- much∈betterdescriptionoftheeffectsofanisotropicvelocity
ing the evolution, which would lead to a secular deviation . dispersion on satellite decay, independently of the orbital
Bydefinition,χisequivalenttothediscrepancybetweenthe inclination. The variation range of the Coulomb logarithm
numerical and semi-analytical position evolution per orbit. becomes larger for a wider inclination spread, but χ barely
Theselectionofthemaximumandminimumgalacto-centre presents a dependence on the satellite eccentricity if Bin-
distances for comparison permits a direct control over the ney’sformula is applied.
orbitaleccentricityevolution,althoughthemeasureofχcan
IfChandrasekhar’s approach is used,theCoulomb log-
beextended to other pointswithout loss of generality.
arithm that leads to the best fit becomes smaller as the in-
The value of k depends on the objectives of the study. clination increases. Since dynamical friction is proportional
For instance, if the aim is to find the best calibration for a
tolnΛ,theuseoftheaveragevalueimpliesanoverestimate
large period of time, as it may betoreproduce thesatellite of the force for low inclinations and vice versa. The final
decay in spiral galaxies, the k-valuemust include anumber average,
ofperiodslargeenough,sothattheoverallevolutionofsev-
eral satellite orbits can be reproduced with a single value 1
of lnΛ. In this paper, we limit our fit to the first satellite (χ/r0)= NlnΛ Xχi/r0, (9)
lnΛ
periods, namely, k=3,4, for which thedifferences between
both approaches of dynamical friction can beclearly seen. overtheNlnΛ numericalexperimentsofTable2isplottedin
In Fig. 1 we plot χ for some of the experiments, con- Fig. 2. This figure shows the large discrepancies produced
cretely,thosewithinclinations60◦,45◦ and30◦(rows),with byChandrasekhar’s expression if thefit is for a large range
eccentricities0.3,0.5and0.7(columns).Foreachmodel,the of orbital inclinations and eccentricities, as expected. The
semi-analytic code is employed to generate the satellite or- minima of the curves determine the values of lnΛ that lead
bitusingChandrasekhar’s(dashedanddotted-dashedlines) tothebest fits,which wesummarise in Table 3.Thevalues
and Binney’s (full and dotted lines) formula to reproduce of χmin denote the average error during the first k orbits
dynamical friction. This figure shows that Chandrasekhar’s associated with the fit.
formula poorly describes thedependenceof thesatellite or- For the following analysis the value of the Coulomb
bitwiththeinitialinclination,leadingtoawiderdispersion logarithm implemented in our semi-analytic code will be
of the Coulomb logarithm values (for this range of inclina- found in Table 3. Looking at Fig. 1, we expect that Chan-
tions, between 30◦ and 60◦, lnΛ [0.9,2.8]). If Binney’s drasekhar’s friction will lead to more accurate fits for low
expression is used, the variation o∈f lnΛ is highly reduced (i 30◦)thanforhighinclinedorbitsafterfixinglnΛ=2.2.
≃
(cid:13)c 0000RAS,MNRAS000,000–000
6 J. Pen˜arrubia, A. Just & P. Kroupa
thepotentials
Gm
Binney k=3 Φ = s , (10)
i
−Xi6=j |ri−rj|2+ǫ2
Binney k=4 p
Φext = Φg(rs),
| |
Chandr. k=3 the softening being ǫ 0.23mu = 0.8kpc, which is the
≃
resolution of the inner grid focused on the satellite centre-
Chandr. k=4 of-densityrs,and Φg thegalaxy potentialat thispoint,ne-
glectingthetidalterms.Thus,alltheparticlesofasatellite
are assumed to feel the same external potential, which is a
useful and sufficiently accurate approximation, taking into
account that most of the bound particles are located very
close tothis point.
Particles with E > 0 are removed and the procedure
i
repeated until only negative energy particles are left. As
Johnstonetal(1999)show,theenergycriteriumpermitsone
todistinguishthoseparticlesthat,thoughunbound,remain
in orbits inside the tidal radius, which will escape from the
satellite after some orbital periods.
The mass is calculated each ∆t = 0.312 Gyr, so that
the semi-analytic code interpolates the value for intermedi-
ate points at each time-step. The error is of the order of
∆M(t)/∆t, going linearly with the mass loss. This means
Figure2.Averageofthefittingparametersoverthecalculations that theinterpolation might introduce not negligible differ-
ofTable2. ences at times where the mass loss is significant (i.e late
times of the satellite evolution). In this study we are not
concerned in detail with thelate phases of evolution.
In the right columns of Fig. 3a, b and c we plot the
Friction k lnΛ χmin/r0 χmin/(kr0) mass evolution for different orbits. Most of the satellites
Binney 3 2.4 0.024 0.0080 reach the inner most regions of the parent galaxy with a
4 2.4 0.038 0.0095 substantial fraction of their initial mass. A comparison of
these curves with those of PKB (where a bulge and a disc
Chandrasekhar 3 2.1 0.147 0.049
component were included) shows that the baryonic subsys-
4 2.2 0.193 0.048
tems of a galaxy induce a larger mass loss through tidal
heating (e.g, Taylor & Babul 2001, Pen˜arrubia 2003). PKB
Table3.Resultsofthefittingprocedureappliedtothenumerical
observe in their numerical experiments that all satellites
calculations of Table 2 for both formulæ of dynamical friction.
Thefifthcolumngivesrelativedeviationperorbit. with Ms = 0.1Md and r0 = 55 kpc are destroyed before
theremainingboundpartofthesatellitereachesthecentral
region of thegalaxy.
7 THE VELOCITY ANISOTROPY EFFECTS
7.2 Satellite decay
In thisSection, we discuss in more detail various aspects of
One of the most important effects of dynamical friction is
satellite evolution.
the monotonic reduction of the orbital angular momentum
and energy during the satellite’s evolution that leads to a
progressivedecreaseoftheaveragedgalacto-centredistance.
ThecomputationscarriedoutbyPKBshowastrongdepen-
7.1 Satellite mass loss
denceof the decay time on theinitial inclination that must
Tidal forces induce satellite mass loss. The satellite mass becompared to analytic estimates.
plays an important role in determining the ultimate fate In Fig. 3a, b and c we plot the radius evolution (left
of its evolution and survival (e.g. Vel´azquez & White 1999, columns) for those models with initial e = 0.5,0.7 and 0.3,
Klessen & Kroupa 1998, Johnston et al. 1999). The value respectively.Fromthisfigure,weconcludethatBinney’sex-
of M (t)is treated as an input bythe semi-analytical code, pression clearly produces more accurate results than Chan-
s
and is calculated using theself-consistent Superbox code. drasekhar’s forthewhole range of orbital inclinations. This
The mass remaining bound to the satellite, M (t), is result is not surprising due to the small dependence of the
s
determinedintheSuperboxcalculationsbycomputingthe Coulomb logarithm on the inclination and eccentricity as
potential energy Φ < 0 of each satellite particle presumed shown in Fig. 1. Additionally, the value of lnΛ that pro-
i
bound to the satellite, and its kinetic energy (T ) in the duces the best fit for the first few orbits also succeeds in
i
satellite frame. Following PKB, particles with E = T + reproducing thedecay time of thesatellite.
i i
ms(Φi+Φext) > 0 are labelled unbound, where ms is the PKB observe that coplanar satellites suffer larger fric-
mass of one satellite particle (all have the same mass) and tion than those following polar orbits, leading to survival
(cid:13)c 0000RAS,MNRAS000,000–000
Dynamical friction in flattened systems: A numerical test of Binney’s approach 7
60 60
1 1
40 40
0.5 0.5
20 20
600 0 600 0
1 1
40 40
0.5 0.5
20 20
600 0 600 0
1 1
40 40
0.5 0.5
20 20
600 0 600 0
1 1
40 40
0.5 0.5
20 20
600 0 600 0
1 1
40 40
0.5 0.5
20 20
0 0 0 0
0 2 4 6 0 2 4 6 0 1 2 3 4 5 0 1 2 3 4 5
t (Gyr) t (Gyr) t (Gyr) t (Gyr)
Figure3.a:RadiusandmassevolutionforthemodelsofTable2 Figure3–continuedb:AsFig.3aforthosesatellitesofTable2
withinitiale=0.5.Opencirclesdenotethenumericalevolution, withinitiale=0.7.(Notethatthetime-axishaschangedscale.)
whereasfullanddottedlinesrepresentthedataobtainedfromthe
semi-analyticcode using, respectively, Binney’s (lnΛ=2.4) and
Chandrasekhar’s(lnΛ=2.2)expressionstoreproducedynamical 60
friction. 1
40
0.5
20
times 70% shorter. Due to the presence of a disc in their
galaxymodel,thecontributionofthediscanisotropyonthe
600 0
decay differentiation as a function of the inclination can- 1
not be directly measured. The calculations presented here 40
(wherethediscandbulgeareremoved)showsurvivaltimes 0.5
20
that range from 3.7 Gyr (coplanar orbits) to 6 Gyr (polar
orbits),using thesame orbitalparameters and haloflatten- 600 0
ing as PKB. This implies a decay time difference of around 1
40
60% between polar and coplanar satellites, which indicates
that the disc contribution might be of the order of 10%. 0.5
20
The effects of the disc on the satellite orbit can be found
in Pen˜arrubia (2003) and will be addressed in a following 600 0
1
paper.
40
Depending on the symmetry of the halo distribution,
0.5
one can observethe following effects: 20
Spherical mass distribution and isotropic veloc- 600 0
•
ity distribution:Satellitesorbitingsystemswithaspheri- 1
40
caldistributionfunctionmoveonorbitsthatdonotdepend
on theirorientation with respect to thesymmetry axis. 0.5
20
Flattenedmass distribution and isotropic veloc-
•
ity distribution:Bymeansofthelocalapproximationthe 0 0
0 2 4 6 8 0 2 4 6 8
value of dynamical friction is determined by the properties
t (Gyr) t (Gyr)
of the halo at the satellite’s position, being reproduced by
Chandrasekhar’sformula.Ourresultsindicatethatthespa- Figure3–continuedc:AsFig.3aforthosesatellitesofTable2
tialasphericityleadstoastrongdifferentiation ofthesatel- withinitiale=0.3.(Notethatthetime-axishaschangedscale.)
lite decay as a function of theorbital inclination.
Flattened mass distribution and anisotropic ve-
•
(cid:13)c 0000RAS,MNRAS000,000–000
8 J. Pen˜arrubia, A. Just & P. Kroupa
locity distribution: The main influence of the velocity
1
anisotropy on the satellite orbit is the secular reduction of 80
0.8
the orbital inclination (see Subsection 7.3) and the reduc- 60
0.6
tion of the spatial anisotropy effects, which is equivalent 40 0.4
to BR < Bz in Binney’s formula (eq. 3, oblate systems). 20 0.2
Taking into account the velocity anisotropy explicitly thus 0 0
1
reduces the difference in orbital decay times between satel- 80
0.8
60
lites on polar and coplanar orbits. As Fig. 3 shows, assum- 0.6
40
inganisotropicdistributioninvelocityspace(BR =Bz)or, 0.4
20
equivalently, using Chandrasekhar’s formula to reproduce 0.2
dynamical friction, leads to an overestimation of dynamical 0 01
80
frictionforlowinclinedsatellitesandtoanunderestimation 0.8
60
for those satellites following highly inclined orbits. This is 0.6
40
duetothefact that theanisotropyleads toareducedeffec- 0.4
20
0.2
tiveX forhighinclinedorbitsyieldinganenhancedfriction
0 0
forceandviceversaforlow inclination.Thiscompeteswith 1
80
the density effect due to the flattening, because highly in- 60 0.8
clined orbits have lower mean density along the orbit com- 40 0.6
0.4
pared to lower inclination. 20
0.2
0 0
The evolution of the orbital radius of satellites with 1
80
initial e = 0.3,0.7 is plotted in Fig. 3b and c. As we can 0.8
60
observe, Binney’s approximation reproduces accurately the 0.6
40
0.4
overallradiusevolutionindependentoftheinitialeccentric- 20
0.2
ity and orbital inclination. 0 0
0 2 4 6 0 2 4 6
t (Gyr) t (Gyr)
7.3 Evolution of the orbital inclination and Figure4.a:Inclinationandeccentricityevolutionforthemodels
eccentricity of Table 2 with e = 0.5. Left column: Dotted lines represent
the N-body inclination evolution, full and dashed lines denote
Orbits in non-spherical systems have inclinations (i theuseofBinney’sandChandrasekhar’sequations,respectively.
≡
arccos[L /L]) that do not remain constant but suffer pe- Right column: Numerical (dotted lines) against semi-analytical
z
riodicaloscillations duetonutation.Inaddition,theorbital eccentricity evolution. Solid circles and open squares denote the
plane precesses at a constant rate. Once the initial condi- implementation of Binney’s and Chandrasekhar’s expressions in
thesemi-analyticcode,respectively.Wenotethat,forclarity,the
tions are fixed, the amplitude and frequency of nutation
semi-analyticvalues ofeareonlyplottedattheapo-centres.
and theprecession rateremain constant if thefriction force
is removed from the equations of motion, whereas, if im-
plemented, nutation and precession vary according to the
angular momentum and radial distance evolution. Our in- urealsoconfirmsthattheorbitalinclinationofcoplanarand
terest focuses now on the effects induced by the velocity polar satellites remains constant.
anisotropy on thesatellite inclination duringthe orbit. If dynamical friction is modelled by Chandrasekhar’s
Binney (B77) predicts a progressive reduction of i due formula, i.e the velocity distribution is assumed to be
to dynamical friction if the velocity dispersion ellipsoid is isotropic,theaveragedvalueofidoesnotchangeduringthe
axi-symmetric (σ ,σ ) and σ > σ . By symmetry, the orbit, which contradicts the numerical results. Thus, while
R z R z
inclination decrease will not occur if the orbits are either thedifferencebetweenthesurvivaltimesofpolarandcopla-
coplanar (i=0◦) or polar (i=90◦). nar orbits is larger for flattened DMHs treated with Chan-
Theinclinationevolutionofmodelswithe=0.5isplot- drasekhar’s dynamical friction formula (Fig. 3), taking into
ted in Fig. 4 (left column), where dotted lines denote the accounttheDMHanisotropic velocitydistribution function
numericaldataandsolidanddashedlinesthesemi-analytic explicitly via Binney’s formulae leads to an increase with
evolution if dynamical friction is modelled by Binney’s and time of the anisotropy of the satellite distribution due to
Chandrasekhar’s formulæ, respectively. This Figure shows the kinematical coupling of the satellites to the DMH ve-
thereductionofthemeanvalueofipredictedbyBinneyand locity field.Thisclearly agrees with thefully self-consistent
observed byPKB intheirnumerical experiments.Afterthe N body computations reported here.
−
satellitehassunktotheinnermostregionofthehalo,thein- In Fig. 4b and 4c (left columns) we plot the compar-
clinationsareaslowas10◦ 20◦ independentoftheirinitial ison for models with e = 0.7,0.3, respectively. The results
−
value.ThislargedecreaseofiiswellreproducedbyBinney’s show barely adependenceon theeccentricity.It isinterest-
expression, although the nutation process shows discrepan- ingto notethat, independentlyof e,orbits that are neither
cies with the numerical result, which is connected with the coplanar nor polar present large drops of the mean value
notexactreproductionoftheorbit.Despitetheaccuratede- of i. After the satellite sinks to the centre, the final orbital
scription oftheoveralldecayprocess,thisorbitalmismatch inclination lies in between 10-20◦ for all the models.
is also observed when applying Chandrasekhar’s expression We must emphasise the accuracy of Binney’s formula
in spherical systems (see Just & Pen˜arrubia 2002). The fig- in describing correctly the inclination decrease that satel-
(cid:13)c 0000RAS,MNRAS000,000–000
Dynamical friction in flattened systems: A numerical test of Binney’s approach 9
lites suffer in oblate systems. This is crucial for simulating
1
80 properlysatellitemotionsandforinvestigatingsatellitedis-
0.8
60 tributionsaround spiral galaxies.
0.6
40 0.4 Liketheorbitalinclination,theeccentricityisoneofthe
20 0.2 orbitalparametersthatcanbeindirectlymeasuredfromob-
0 0 servationstodetermineasatellite’smotionaroundagalaxy.
1
80
IntherightcolumnsofFig.4weshowthecomparisonofthe
0.8
60
0.6 numericaleccentricityevolutionwithbothsemi-analyticap-
40
0.4 proaches.Theanalyticformulæofdynamicalfrictionleadto
20
0.2 a larger eccentricity decrease, which occurs mostly at late-
0 01 times of the orbit, the so-called orbital circularisation, and
80
0.8 becomes stronger for low inclined orbits, those that suffer
60
0.6 higherdynamicalfriction. Fig.4bandFig. 4cindicatethat
40
0.4 the circularisation is more pronounced if the initial orbital
20
0.2
eccentricityishigheranddecreasesformorecircularorbits.
0 0
1 Bothdynamicalfriction expressionslead toasimilareccen-
80
60 0.8 tricityevolutionforthefirstfeworbitalperiods,howeverat
40 0.6 late-times the eccentricity exhibits a reduction not present
0.4
20 in the numerical calculations that can be as high as 80%
0.2 ∼
0 0 forlowinclinedsatellitesonhighlyeccentricorbits(Fig.4b,
80 1 model H2S100e).
0.8
60
0.6
40
0.4 7.4 Energy and angular momentum evolution
20
0.2
0 0 A flattened system possesses two analytic constants of mo-
0 1 2 3 4 5 0 1 2 3 4 5
tion,theenergy and thecomponent of theangular momen-
t (Gyr) t (Gyr)
tumperpendiculartotheaxi-symmetryplane(whichwede-
Figure4–continuedb:AsFig.4aformodelswithe=0.7.Note note as Lz). The total angular momentum L2 = L2R +L2z
thatthetime-scalehasadifferentvalue. is, however, not constant during a satellite orbit (see e.g
Binney & Tremaine 1987), but has periodic variations that
correspondtoaprecessionandnutationoftheorbitalplane
1 around thez-axis.
80
0.8 Sincethedynamicalfrictionforcehasanoppositesense
60
0.6 with respect to the satellite velocity, it decreases the angu-
40
0.4 lar momentum and the energy which induces a monotonic
20
0.2
sinking into the inner regions of the halo potential. The re-
0 0
1 ductionofangularmomentum,therefore,impliesanincrease
80
60 0.8 ofthebindingenergy(inabsolutevalue),sincethepotential
40 0.6 grows with decreasing radius. Due to the small magnitude
0.4
20 of dynamical friction compared to the mean field force, we
0.2
0 0 expect an easier comparison between numerical and semi-
80 1 analytic data by the slow variation of Lz and E during the
0.8
60 orbit.
0.6
40 In Fig. 5 we plot the changes of specific E and L due
0.4 z
20 to dynamical friction for the models with e = 0.5. The re-
0.2
0 0 sults are equivalent to those of the radial evolution. Dueto
1
80 ourselection oflnΛastheaverageoverthoseCoulomb log-
0.8
60 arithms that lead to the best fit for each particular model,
0.6
40 0.4 Chandrasekhar’s formula overestimates dynamical friction
20 0.2 for low inclined orbits and underestimates it for highly in-
0 0 clined orbits. For orbits with i < 30◦, this appears as a
1
80
0.8 stronger reduction of the z-component of angular momen-
60
0.6 tumand,equivalently,alargeincreaseofthebindingenergy.
40 0.4 The effect is contrary for satellites with i>30◦.
20
0.2 This figure illustrates how the kinetic energy of the
0 0 satelliteislost viafriction, beingtaken-uppartiallybyhalo
0 2 4 6 8 0 2 4 6 8
t (Gyr) t (Gyr) particles and also being deposited in unboundsatellite par-
ticles. At the end of the simulation the angular momentum
Figure 4 – continued c: As Fig. 4a for models with initially has a null value, i.e the satellite remains in the inner most
e=0.3.Notethatthetime-scalehasbeenrescaled. part of thegalaxy.
It isinteresting tonotethat thenumericalevolution of
energy and angular momentum presents small oscillations
(cid:13)c 0000RAS,MNRAS000,000–000
10 J. Pen˜arrubia, A. Just & P. Kroupa
apo-centresforagivennumberkoforbits.Ifdynamicalfric-
1 tion is modelled by Binney’s equation, this quantity shows
3
2 discrepancies of approximately χmin=0.009r0 per orbit af-
0.5 ter averaging over the set of experiments and for the first
1 threeorbits,whileChandrasekhar’sformulaproducesvalues
0 0 of around χ=0.05r0 per orbit (Table 3).
1 3 We conclude that Binney’s expression faithfully re-
produces the process of dynamical friction in anisotropic
2
0.5 systems. The fit is as accurate as that employing Chan-
1
drasekhar’s formula in isotropic systems (see Just &
0 0 Pen˜arrubia 2002).
1
3 The comparison of orbits resulting from Chan-
2 drasekhar’s and Binney’s expression of dynamical friction
0.5
1 gives us the possibility to asses the effects of the DMH ve-
locity anisotropy on satellite dynamics. We have demon-
0 0
strated that, (i) if the density profile is in both equations
1
3
ρ = ρ(R,z), where R,z are the cylindrical coordinates of
2 the satellite position vector, the orbits generated by Chan-
0.5
1 drasekhar’s formula overestimate the decay time for polar
orbits and underestimate it for coplanar ones for an over-
0 0
1 all best-fit Coulomb logarithm (lnΛ = 2.2). One effect of
3
the velocity anisotropy is then to reduce the interval of de-
2
0.5 cay times as a function of the orbital inclination. Binney’s
1 expression has been shown to reproduce accurately the nu-
0 0 merical results independently of the initial eccentricity and
0 2 4 6 0 2 4 6
with lnΛ = 2.4. (ii) Dynamical friction in systems with
t (Gyr) t (Gyr)
anisotropic velocity distribution leads toa markeddecrease
oftheorbitalinclination iwhich iswell reproduced byBin-
Figure5.Specificenergyandangularmomentumevolutiondur-
ney’s expression. After the satellite sinks to the inner most
ing the orbits with e = 0.5. The numerical evolution is denoted
bydottedlines,whereasthesemi-analyticdataisrepresentedby region of the galaxy, i lies within 10-20◦, independently of
solid and dashed lines if dynamical friction is modelled by Bin- the initial value unless i 90o (polar orbits). (iii) The en-
≈
ney’sandChandrasekhar’sformulæ,respectively.Thequantities ergy and angular momentum evolution as a function of the
E and Lz are normalised to the initial value. Note that for the orbital inclination confirm theresults of (i) and (ii).
casei=90◦ onehasLz =0. Thesemi-analyticeccentricityevolution,eitheremploy-
ing Chandrasekhar’s or Binney’s formula, shows the so-
called circularisation process, defined as the progressive re-
during the orbit. This behaviour is due to the self-response
duction of e during the orbit. This variation is stronger
ofthehalotothesatellitemotion.Sincesuperboxpreserves
for increasing friction (as for coplanar orbits or during the
thetotalenergyandangularmomentum,thehaloalsomoves
late-times of the evolution) and barely takes place in the
around the centre-of-mass of the system. Due to the com-
self-consistent N body calculations. The small N body
plexityofthefeedback,itcannotbereproducedanalytically − −
circularisation agrees with the results of van den Bosch
(the halo centre-of-mass is fixed at the coordinate origin in
et al. (1999). A possible solution may be sought in the
thesemi-analytic code).
position-dependenceoftheCoulomblogarithm,asproposed
by Hashimoto, Funato & Makino (2003). However, despite
the improvement in the description of the orbit at early-
times, the scheme of Hashimoto et al. overestimates the
8 CONCLUSIONS
satellite decay time for all the experiments (see Just &
To asses the accuracy of Binney’s equations (B77) and in Pen˜arrubia 2002). We must conclude that the reason for
ordertoreproducethedecayofsatellitesinflattenedDMHs circularisation in the semi-analytic orbits is not yet fully
with a semi-analytical code, we have performed a set of understood.
self-consistent numerical experiments for different orbital
inclinations and eccentricities of the satellite. The semi-
analyticcodeincorporatesdynamicalfrictioneitherinterms
ofChandrasekhar’sexpressionorasBinney’sformulae.Both
treatmentsincludetheasphericaldensityprofilebymeansof
thelocalapproximation.Thismeansthatthedifferenceson
thesatellitemotioninducedbyeachtreatmentofdynamical
9 ACKNOWLEDGEMENTS
friction comes from the anisotropy in velocity space, which
is implemented in theanalysis of B77. JP acknowledges support through a SFB 439 grant at the
The accuracy of Binney’s and Chandrasekhar’s for- UniversityofHeidelberg.PKacknowledgessupportthrough
mulæincomparisonwiththenumericalorbitsisdetermined DFGgrant KRI635/4. Wethanktheanonymousrefereefor
bytheparameterχ2 = 1 (∆r2+σ2∆t2)attheperi-and his/her advice and useful comments.
2k
P
(cid:13)c 0000RAS,MNRAS000,000–000