Actions, angles and frequencies for numerically integrated orbits

2 MB·

by  Jason L. Sanders

Mon.Not.R.Astron.Soc.000,1–12(2014) Printed16June2014 (MNLATEXstylefilev2.2)
Actions, angles and frequencies for numerically integrated orbits
Jason L. Sanders(cid:63) & James Binney
RudolfPeierlsCentreforTheoreticalPhysics,KebleRoad,Oxford,OX13NP,UK
4 16June2014
1
0
2 ABSTRACT
Wepresentamethodforextractingactions,anglesandfrequenciesfromanorbit’stimeseries.
n
u The method recovers the generating function that maps an analytic phase-space torus to the
J torus to which the orbit is confined by simultaneously solving the constraints provided by
each time step. We test the method by recovering the actions and frequencies of tori in a
3
1 triaxialSta¨ckelpotential,anduseittoinvestigatethestructureoforbitsinatriaxialpotential
thathasbeenfittedtoourGalaxy’sSagittariusstream.Themethodpromisestobeusefulfor

] analysingN-bodysimulations.Italsotakesasteptowardsconstructingdistributionfunctions
A
forthetriaxialcomponentsofourGalaxy,suchasthebaranddarkhalo.
G
Keywords: methods:numerical–Galaxy:kinematicsanddynamics–galaxies:kinematics
.
h anddynamics
p
-
o
r
t
s 1 INTRODUCTION tothreefunctionallyindependentisolatingintegrals,Jeanstheorem
a statesthatweshouldbeabletorepresentanequilibriumstellarsys-
[ Althoughnogalaxyiseverinperfectdynamicalequilibrium,equi-
tembythedensityofstarsinathree-dimensionalspaceofintegrals
libriumdynamicalmodelsarecentraltotheinterpretationofobser-
3 ratherthaninfullsix-dimensionalphasespace.Thisreductionin
vationsofbothourGalaxyandexternalgalaxies.Amajorreason
v dimensionalitymakesthesystemverymucheasiertocomprehend
for the importance of equilibrium models is that we can infer a
0 andmodel.
galaxy’s gravitational potential, and thus its dark-matter distribu-
0 Sinceanyfunctionofintegralsisitselfanintegral,infinitely
6 tion,onlytotheextentthatthegalaxyisinequilibrium.Moreover,
manydifferentintegralsmaybeusedasargumentsoftheDF.How-
3 equilibriummodelsarethesimplestmodelsandmorecomplexcon-
ever,theactionintegralsJ standoutasuniquelysuitedtobeused
. figurations,involvingspiralstructureoranon-goingminormerger i
1 asargumentsoftheDF.Whatmakesactionsspecialisthattheycan
forexample,arebestmodelledasperturbationsofanequilibrium
0 becomplementedbycanonicallyconjugatevariables,theanglesθ ,
model. i
4 toformacompletesetofcanonicalphase-spacecoordinates.Ob-
1 Globular clusters are the stellar systems that are most com- viouslytheequationsofmotionoftheactionsaretrivial:J˙ = 0.
: pletelyunderstood,andthetheoryofthesesystemsillustratesthe i
v Moreremarkablytheequationsofmotionoftheanglesarealmost
importanceofequilibriummodels:ateachinstanttheclusterisas-
Xi sumed to be in dynamical equilibrium, so, by Jeans’ theorem, its as trivial: θ˙i = Ωi(J) = constant. Thus the angle variables in-
creaselinearlyintimeandifweuseangle-actioncoordinates,the
distributionfunction(DF)isafunctionoftherelevantisolatingin-
r unperturbedmotionofstarsbecomestrivial.Thisfactmakesangle-
a tegrals, such as stellar energy E, total angular momentum L, or
actioncoordinatesuniquelysuitedtoworkinvolvingperturbation
angularmomentumaboutasymmetryaxis,L .Overmanydynam-
z theory,andindeedtheangle-actioncoordinatesoftheKeplerprob-
icaltimesencountersbetweenstarsandstellarevolutioncausethe
lemwereinventedtoexploretheroleplayedbyplanet-planetinter-
DF to change, but in such a way that the DF continues to satisfy
actionsinthedynamicsoftheSolarSystem.
Jeans’theorem,sotheclusterevolvesthroughaseriesofdynami-
McMillan&Binney(2008)haveshownthatangle-actionco-
calequilibria.
ordinatesmakeitpossibletoidentifystarsneartheSunthathave
N-bodysimulationsofcosmologicalclusteringlikewiseyield
been stripped from an object that was tidally disrupted gigayears
apictureinwhichdark-matterhaloesarefarfromdynamicalequi-
ago,andeventodeterminethedateofthedisruptiontogoodpre-
librium only during short-lived and quite rare major mergers. In
cision. Sellwood (2010) and McMillan (2013) have used angle-
generaladark-matterhalocanbewellapproximatedbyadynami-
actioncoordinatestoidentifystarsneartheSunthatareresonantly
calequilibriumthatismildlyperturbedbyaccretion.
trappedbyspiralstructure.Sanders&Binney(2013)haveshown
ThenaturalwaytomodeladynamicalequilibriumisviaJeans
howangle-actioncoordinatesforthestarsofastreamenableone
theorem, which assures us that the system’s DF can be assumed
toconstrainthegravitationalpotentialinwhichthestreammoves.
tobeanon-negativefunctionofisolatingintegrals.Sinceoneex-
pectsasmoothtime-independentgravitationalpotentialtoadmitup Cosmologicalsimulationshaveshownthattriaxialdarkmatter
haloesaretobeexpected,atleastuptothepointatwhichbaryons
become gravitationally dominant (Valluri et al. 2010). Moreover,
(cid:63) E-mail:[email protected] Law & Majewski (2010) and Vera-Ciro & Helmi (2013) present
c 2014RAS
(cid:13)
2 J.L.Sanders&J.Binney
evidencethatthetidaltailsoftheSagittariusdwarfgalaxycanonly shortaxis,whilstJ quantifiesoscillationparalleltothisaxis.Fora
3
be fitted if the Milky Way has a triaxial dark matter halo. Hence long-axislooporbit,J quantifiescirculationaroundthelongaxis,
3
weneedtobeabletodetermineangle-actioncoordinatesforstars whilstJ quantifiesoscillationparalleltothisaxis.Wechoosethis
2
in triaxial potentials. In this paper we show how to evaluate the definitionsuchthatouractionsmatchJ ,J andJ foraSta¨ckel
λ µ ν
angles and actions of particles in a given triaxial potential. If the potential(deZeeuw1985),andeachclassoforbitoccupiesadis-
potential is axisymmetric, the actions can be evaluated using the tinctregionofactionspace(seeSection4).
algorithmgivenbyBinney(2012a).
InSection2wederivetheequationsthatyieldvaluesofan-
2.1 Toypotentials
gles, frequencies and actions. In Section 3 we test our solutions
oftheseequationsbycomparingtheresultingangles,frequencies 2.1.1 Triaxialharmonicoscillator
andactionsfortwoorbitsinaSta¨ckelpotentialwithanalyticval-
Forboxorbitsweusethepotentialofthetriaxialharmonicoscilla-
ues.InSection4weusetheequationstoexploreaconstant-energy
tor,
surfaceoftheactionspaceofthetriaxialpotentialforourGalaxy
thatLaw&Majewski(2010)fittedtothetidalstreamoftheSagit- 3
(cid:88)
tariusdwarf.InSection5werelateourworktopreviousworkin Φho(x)= 12 ωi2x2i, (1)
thefield,discussapossibleextension,andexplorehowthemethod i=1
copeswithresonanttrapping.Section6sumsupandlookstothe whichhasthreeparameters,ω .Herewehavechosentheprincipal
i
future. axesofthepotentialtoliealongtheCartesianx,y,zdirectionson
theassumptionthatthetimeserieshasalreadybeenrotatedintothe
coordinatesystemthatisalignedwiththeprincipalaxesofthetrue
potential.Theactionsandanglesinthispotentialaregivenby
2 FORMALISM
p2+ω2x2
Anglesandactionscanbeassignedtoorbitsthatare“regular”or Ji = i 2ωi i,
i (2)
quasiperiodicbecausesuchanorbitisconfinedtoatoruslabelled (cid:16) p (cid:17)
bytheactions(Arnold1978).Wewillworkinthreedimensionsso θi =arctan ω xi .
i i
willhavethreeactionsdenotedasJ = (J ,J ,J ).Eachaction
1 2 3
quantifiesthemagnitudeoftheoscillationinasuitablecoordinate.
2.1.2 Isochronesphere
The transformation from ordinary phase-space coordinates
(x,v) to angle-action coordinates (θ,J) is possible analytically Forlooporbitsweusetheisochronepotential,
in only a few cases. McGill & Binney (1990) used one of these
−GM
casesasastartingpointforthenumericalconstructionofmoregen- Φ (x)= √ , (3)
iso
b+ b2+r2
eraltransformationsby“torusmapping”.Thekeypointabouttorus
mappingisthatityieldsorbitswithspecifiedactionsratherthanor- whereristhesphericalradius.Thispotentialhastwofreeparam-
bitswithspecifiedinitialconditions(x,v).WhenanalysinganN- eters:themassM andthescaleradius,b.Theexpressionsforthe
bodymodel,werequireactionsgivenaninitialconditionandnot actionsandanglesinthispotentialaremoreinvolvedthanforthe
viceversa.HereweadapttheapproachofMcGill&Binney(1990) harmonicoscillatorsoarenotrepeatedhere.Readerscanconsult
into a procedure which finds the actions, angles and frequencies Binney&Tremaine(2008)fortheappropriateequations.Thethree
givenaseries ofphase-spacecoordinates(xi,vi)sampledalong actionsintheisochronepotentialaregivenbytheradialactionJr,
anorbitattimesti,where0(cid:54)ti (cid:54)T. thez-componentoftheangularmomentumLzandtheverticalac-
Withthistimeseriesweseekageneratingfunctionthatwill tionJz ≡L−|Lz|,whereListhetotalangularmomentum.With
mapa“toytorus”ofasimple“toypotential”intothe“targettorus” thischoicewemustorientourcoordinatesystem,suchthattheor-
towhichtheorbitisconfined.Thetoypotentialmusthaveanalyt- bitcirculatesaroundthez-axis,beforefindingtheactions.
icallytractableanglesandactionsandpermitorbitsthathavethe
correctgeometry.
2.1.3 Offsets
In the absence of figure rotation, a general triaxial potential
admitstwobasicclassesofnon-resonantorbit:looporbitsandbox Onemightalsoincludetheoffsetofthecentreofthepotentialfrom
orbits (Schwarzschild 1979; de Zeeuw 1985). Loop orbits have a the coordinate centre as a free parameter, but we shall not do so
definite sense of rotation either around the long- or short-axis of here,presuminginsteadthatthetimesamplesx havealreadybeen
i
the potential, whilst a box orbit has no sense of rotation and can adjustedtoberelativetoone’sbestestimateofthecentreofthetrue
reachdowntothecentreofthepotential.Hencetheclassofanor- potential.
bitcanbedeterminedbyinspectionofcomponentsoftheangular
momentum along the orbit: if all components of the angular mo-
2.1.4 Parameterchoice
mentumchangesign,theorbithasnosenseofcirculationandisa
boxorbit;whenacomponentoftheangularmomentumretainsits Onceaclassofpotentialhasbeenchosen,wesettheparametersof
sign,theorbitisalooporbitaroundthecorrespondingaxis(Carpin- thepotentialbyminimizing(McGill&Binney1990)
tero&Aguilar1998).Foreachclassoforbitweuseatoypotential
(cid:88)
thatprovidestoriwiththesamegeometricalstructureasthetoriof χ2 = (Hi−(cid:104)H(cid:105))2, (4)
thegivenorbitclass. i
ForaboxorbittheactionsJ ,J andJ quantifytheoscil- wherethesumisoverthetimes,H isthevalueofthetoyHamil-
1 2 3 i
lation in the x, y and z directions, respectively. For loop orbits, tonianat(x ,v ),and(cid:104)H(cid:105)isthemeanofthesevalues.Themin-
i i
J quantifies oscillation in a generalized radial coordinate. For a imizationofχ2isdoneusingtheLevenberg–Marquardtalgorithm
1
short-axis loop J quantifies the particle’s circulation around the (Pressetal.2002).
2
c 2014RAS,MNRAS000,1–12
(cid:13)
Actionsandanglesforintegratedorbits 3
Theexperimentsdescribedbelowsuggestthatthismethodfor wheretheinnersumisoverthedimensionoftheactionspaceand
selectingtheparametersissub-optimalinthatitleadstoarather thesetNislimitedtoafinitenumberofvectorsn.Wetakethisset
centrally concentrated toy potential being selected. This central tobetheN vectorsthatsatisfythecondition|n| (cid:54) N ,where
max
concentration then leads to high-order Fourier components being N (cid:39)6.
max
requiredinthegeneratingfunction.However,ourattemptstofind WeminimizeEbysettingtozeroitsderivativeswithrespect
abetterprocedureforselectingthetoypotentialhavenotmetwith totheunknowns:
success.
∂E
0=
∂J
k(cid:48)
2.2 GeneratingFunction (cid:88)(cid:16) (cid:88) (cid:17)
=−2 Jk(ti)−Jk(cid:48) −2 nkSn(J(cid:48))cosn·θ(ti)
Withatoypotentialchosen,weconstructthegeneratingfunction
i n N
totransformbetweentheangle-actions(θ,J)ofthetoypotential, ∂E ∈
andthose(θ(cid:48),J(cid:48))ofthetargetpotential.Thegeneratingfunction 0= ∂S
m
forthistransformation,S(θ,J(cid:48)),canbewritten (cid:88)(cid:88)
=−2 2m cosm·θ(t )
k i
S(θ,J(cid:48))=θ·J(cid:48)−i(cid:88)Sn(J(cid:48))ein·θ, (5) i k
(cid:16) (cid:88) (cid:17)
n(cid:54)=0 × Jk(ti)−Jk(cid:48) −2 nkSn(J(cid:48))cosn·θ(ti) .
wherethevectornhasintegercomponents.Thefirsttermonthe n N
∈
rightgeneratestheidentitytransformation,whilstthestructureof (11)
thesecondpartisrequiredbytheperiodicityoftheanglevariables.
McGill&Binney(1990)showthatiftheHamiltonianistime- Tosolvetheseequationswedefineamatrixcnk thathasas
reversible,therealityofthegeneratingfunctionrequirestheS to subscripts the vector n and the integer k = 1,2,3 that selects a
n
satisfy particularspatialdimension.ThisN-by-3matrixis
(cid:0) (cid:1)
S =−S . (6) c (t )≡2n cos n·θ(t ) , (nosumovern). (12)
n n nk i k i
−
For this condition to be satisfied there must exist a point on the Wefurtherdefinetwo(3+N)-vectors
toy torus at which J˙ = 0 – in Appendix A we demonstrate that
(cid:88)
thisistrueforthetoypotentialsoftheprevioussection.Withthis xJ ≡(J(cid:48),Sn), bJ ≡ (J(ti),cn(ti)·J(ti)), (13)
constraint,thegeneratingfunctioncanbewrittenas
i
(cid:88)
S(J(cid:48),θ)=θ·J(cid:48)+2 Sn(J(cid:48))sinn·θ, (7) andthesymmetricmatrix
n∈N (cid:88)(cid:18) I cT(t ) (cid:19)
wheretheintegervectorsnarenowrestrictedtojusthalfofathree- AJ ≡ c(3t ) c(t )·ciT(t ) . (14)
dimensionallattice.WetakethishalftobethesetN={(i,j,k)}, i i i i
whereeither(k > 0),(k = 0,j > 0)or(k = 0,j = 0,i > 0).
Here I is the 3-by-3 identity matrix. With these definitions, the
Symmetries of the target potential require some of the S to be 3
n equations(11)tobesolvedcanbewrittenas
zero.ThisisdiscussedfurtherinAppendixA.
Fromthegeneratingfunction(7)wefindthatthetoyactions A ·x =b . (15)
J J J
are
∂S (cid:88) WesolvetheseequationsforxJ byLUdecomposition(Pressetal.
J = ∂θ =J(cid:48)+2 nSn(J(cid:48))cosn·θ, (8) 2002).
n∈N Asimilarprocedureyieldsthetargetanglesfromequation(9).
andthetargetanglesare Wenotethatattimetitheorbithasθ(cid:48)(ti)=θ(cid:48)(0)+Ω(cid:48)tiwhere
θ(cid:48) = ∂S =θ+2(cid:88) ∂Sn(J(cid:48))sinn·θ. (9) Ωthe(cid:48)iisnitthiaeltpaorginettifnrethqeueonrcbyit,ianntdegθra(cid:48)(ti0o)ni.sTthheeraenlegvleanctosruremspoofnsdqiunagretod
∂J ∂J
(cid:48) n N (cid:48) residualsis
∈
Notethatbythechoiceofourgeneratingfunction,thetargetangle
zero-pointcoincideswiththetoy-anglezero-point. F =(cid:88)(cid:88)(cid:16)θk(cid:48)(0)+Ω(cid:48)kti−θk(ti)−2 (cid:88) ∂∂SJn(J(cid:48))sinn·θ(cid:17)2.
Given the choice of a toy Hamiltonian, we may find the toy i k n N k(cid:48)
∈
actionsandangles(J(ti),θ(ti))ateachtime.Eachtimethenpro- (16)
ducesaseparateequation(8)withcommonunknowns:thetarget
actionsandtheFouriercomponentsofthegeneratingfunction,Sn. The unknowns now are θ(cid:48)(0), Ω(cid:48) and the set of ∂Sn/∂J(cid:48). The
Wecannotsolvetheseequationsexactlybecausewearedeal- requirement of vanishing partial derivatives of F with respect to
ingwithequationsinaninfinitenumberofunknowns.Becausewe theunknownsyieldsthematrixequation,
canincludeonlyafinitenumberoftermsontheright-handsideof
A ·x =b . (17)
eachequation,theright-handsidesshouldnotagreeexactlywith θ θ θ
the left-hand sides, and the correct procedure is to minimize the
ThesesymbolsaredefinedinAppendixB.Thetoyangleswillbe
sum of the squares of the residuals of individual equations. This
sumis 2π-periodic,andwerequirethesameforthetargetanglesθ(cid:48)(0)+
Ω(cid:48)ti.However,inordertosolvethematrixequationwemustfirst
E =(cid:88)(cid:88)(cid:16)Jk(ti)−Jk(cid:48)−2(cid:88)nkSn(J(cid:48))cosn·θ(ti)(cid:17)2, (10) make the θ(ti) from the orbit integration continuously increase,
andthenwesolveforthetargetanglesandtakethe2π-modulus.
i k n N
∈
c 2014RAS,MNRAS000,1–12
(cid:13)
4 J.L.Sanders&J.Binney
2.3 ChoiceofN ,N andT of periods. This is an inevitable drawback of the approach taken
T max
here because we have very little control over the sampling in the
Giventheschemepresentedabove,theonlyquestionsthatremain
toyanglespace.
arehowtoselecttheorbitintegrationtimeT,thenumberoftime
Havingidentifiedamodewhichwillnotbewellconstrained,
samples,N ,touse,andwhatvaluetouseforN ,whichde-
T max onepossibilityistosetS = 0forthismode.However,bydoing
terminesthenumberN ofFouriercomponentswesolvefor.Here n
this,weriskthrowingoutamodewhichissignificant,andthere-
wediscusshowwecanautomaticallychoosetheseparameterssuch
coveryoftheactionsandfrequencieswilldeterioratesoweoptnot
thatwehavegoodrecoveryoftheunknowns.
todothis.
Anecessaryconditionisthatthenumberofunknownsmustbe
lessthanthenumberoftimesamples,N .Fortheactioncalcula- AnotherrequirementisthattheSn (and∂Sn/∂J(cid:48))decrease
T aswegotolargernsuchthatthetruncationatN isvalid.If
tionthenumberofunknownsisapproximatelyN3 /2,whilstfor max
max the S do not decrease with n, this is evidence of aliasing such
the angle-frequency calculation we have ∼ 3N3 /2 unknowns. n
max thatthesehighernmodesarenotwellrecoveredandweexpectthe
Wealsoexpectourabilitytorecovertheunknownstodependupon
actions,anglesandfrequencieswillalsonotbewellrecovered.
thesamplingofthetoyanglespace.
Let us first consider an idealised 1D case. If we were able
to sample uniformly in the toy angle of a 1D system, we would 2.3.1 Procedure
select N points in a single period separated in toy angle by
T
Wewillnowsummarizetheabovediscussionintoaprocedurethat
∆ = 2π/N . With this sampling rate we would be able to con-
T
strain all modes einθ with n∆ (cid:54) π. We can choose to constrain canbeimplemented:
onlytheNmax modeswithn < π/∆asthenwewouldbesuper- • We first select a reasonable Nmax, for instance Nmax = 6 is
samplingthehighestconsideredmodes.Hereweareusingatime usedinthelaterexamples.
seriesthatisaproductofanorbitintegrationsoisnotuniformly • We then integrate for some time T recording at least N =
T
spacedintoyangles–thetoy-angledistributiondependsonthetar- 3N +6timesamples(orN = N +3ifweonlyneedtheac-
T
getHamiltonian,thetoypotentialandthedistributionofsampling tions)suchthatwehaveasmanyequationsasunknowns.Thisis
times.TherecoveryofFouriercomponentsfromnon-uniformsam- alwayssatisfiedifwechoose
plesisdiscussedinMarvasti(2001).Toconstrainmodesfroma1D
9N3
non-uniformsamplingwemustsampleonaverageatorabovethe N =max(200, max). (21)
Nyquistfrequency.Ifwehavetoy-anglesamplesθ werequire T 4
i
• Foreachtimesamplewefindthetoyanglesandcheckthatequa-
n i=(cid:88)NT−1(θ −θ )(cid:54)π, (18) tions(19)and(20)aresatisfiedforeachmode.Ifequation(19)is
NT −1 i+1 i notsatisfied,T ismuchlongerthanthefundamentalorbitalperiods
i=1
andsowerequireafinertimesamplingfromtheorbitintegration.
toconstrainmoden.
If equation (20) is not satisfied then we continue integrating the
Hereweareattemptingtorecovercomponentsfromsamples,
orbituntilthisequationissatisfiedforallthemodes.
θ , in 3D toy-angle space. As we are restricted to using samples
i • We then perform the procedure outlined in Section 2.2 to find
generatedfromanorbitintegration,oursamplingislimitedtosome
the S . We require the S to be decreasing with n such that on
n n
sub-spaceofthefull3Dtoy-anglespace.The3Dsamplingcanbe
theboundariesthevaluesoftheS aresmall.Ifwefindthatthe
n
consideredasaseriesof1Dsamplesinn·θ (wefirstunrollthe
i boundaryvaluesofS arelarge,wehavenotincludedasufficient
n
anglessuchthattheyincreasecontinuously).Inordertorecoverthe
numberofmodesinthegeneratingfunctionsowemustincrease
S fromthistoy-anglesamplingweneedtosatisfytwoconditions:
n N andrepeattheaboveprocedureuntilwearesatisfiedthatall
max
(i) Asinthe1Dcaseweneedtosampleonaverageatorabovethe dominantmodesareincluded.
Nyquistfrequencysuchthat
Aswewillseebelowthisprocedureisveryconservativebutshould
1 i=(cid:88)NT−1 ensure that the recovery of the actions, angles, frequencies and
n·(θ −θ )(cid:54)π. (19) componentsofthegeneratingfunctionareaccurate.
N −1 i+1 i
T
i=1
(ii) Foreveryincludedmode,n,wewouldalsolikeagoodtotal
coverageinn·θ.Wechoosetorequirethatthen·θsamplescover 3 EXAMPLE
thefullrangefrom0to2π:
As a test of the above, let us look at an example. The most gen-
max(n·θ)−min(n·θ)>2π. (20) eralseparabletriaxialpotentialisthetriaxialSta¨ckelpotential(de
Zeeuw1985).Wechoosetoworkwiththeperfectellipsoid,which
If this condition is not satisfied, we are including a mode which hasdensityprofile
willnotbewellconstrainedbythetoy-anglesamplingi.e.theaver-
ρ
ageofcosn·θwillnotbenearzero.Wethereforeexpectthatthe ρ(x,y,z)= 0 , (22)
(1+m2)2
correspondingS willnotbewellrecoveredfromthissampling.
n
ItcouldbethatthisS isnotsignificantsowillnotaffectthere- where
n
coveredactionsandfrequenciessignificantly.However, aconser- x2 y2 z2
vativeapproachwouldensurethatequation(20)issatisfiedforall m2 ≡ + + , a(cid:62)b(cid:62)c(cid:62)0. (23)
a2 b2 c2
includedmodes.
The associated coordinates are confocal ellipsoidal coordinates
The second of these conditions is the stricter. To ensure that the in which the actions can be expressed as one-dimensional inte-
toy-angle sampling satisfies equation (20) when an orbit is near- grals.ThesemaybecalculatednumericallyusingGauss-Legendre
resonant,werequiretimesampleswhichspanaverylargenumber quadrature.Similarlythefrequenciescanalsobedeterminedfrom
c 2014RAS,MNRAS000,1–12
(cid:13)
Actionsandanglesforintegratedorbits 5
20 15 4 4
15 3 3
10
10 2 2
5
5 1 1
c c c c
p p p p
k 0 k 0 k 0 k 0
/ / / /
y −5 z 5 y −1 z −1
10 − 2 2
− 10 − −
15 − 3 3
− − −
20 15 4 4
− 16 8 0 8 16 − 16 8 0 8 16 − 5.0 2.5 0.0 2.5 5.0 − 5.0 2.5 0.0 2.5 5.0
− − − − − − − −
x/kpc x/kpc x/kpc x/kpc
2.0 2.0 2.0 2.0
1.5 1.5 1.5 1.5
π π π π
/ 1.0 / 1.0 / 1.0 / 1.0
2 3 2 3
θ θ θ θ
0.5 0.5 0.5 0.5
0.0 0.0 0.0 0.0
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0
θ1/π θ1/π θ1/π θ1/π
Toyaction Trueaction Toyaction Trueaction
1600 450
1400 400
1 1200 1 350
− −
s 1000 s 300
m m
k 800 k 250
c c
p 600 p 200
k k
/ 400 / 150
J J
200 100
0 50
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
t/Gyr t/Gyr
Figure1.ExamplelooporbitinthetriaxialSta¨ckelpotential–thetoppan- Figure2.AsFig.1butforaboxorbit.
elsshowsinblacktheorbitintegratedinthetestpotential.Thisisashort-
axislooporbitsocirculatesabouttheaxisz=0.Infaintredweshowthe
initialpointintegratedinthebest-fittingisochronepotential.Inthemiddle Forthelooporbitthetrueandrecoveredactionsare
panelsweshowthetoyanglescalculatedateachtimesample.Inthebottom
panelweshowthetoyactionsateachtime-stepasadottedline(blackfor Jtrue =(212.09,1307.54,708.15)kpckms−1
Jth1e,abrlruoewfsomrJar2katnhdeerestdimfoarteJd3a)c.tTiohnes.solidlinesshowthetrueactionsand Jrecov =(213.33,1307.29,709.16)kpckms−1,
andthetrueandrecoveredfrequenciesare
Ωtrue =(21.76474,15.65172,19.33786)Gyr−1
Ωrecov =(21.76508,15.65185,19.33780)Gyr−1.
one-dimensionalintegrals.Equationsforthesequantitiesaregiven
In Fig. 3 we show two cross-sections of n-space showing
indeZeeuw(1985).Hereweworkwiththepotentialwithparam-
etersρ0 = 7.2×108M kpc−3,a = 5.5kpc,b = 4.5kpcand the absolute value of the components of the generating function.
c=1kpc. (cid:12) (For the isochrone potential we use the convention that subscript
1 refers to J , subscript 2 refers to L and subscript 3 refers to
In this potential we examine two orbits – a short-axis r z
J ≡ L−|L |.)Weseethatthetwomostsignificantmodesare
loop orbit with initial condition (x,y,z) = (10,1,8) kpc, z z
(vx,vy,vz) = (40,152,63) kms−1 and a box orbit with ini- n=(−1,2,0),whichcausesamixingbetweentheradialmotion
andazimuthalmotion,andn=(0,0,2).NotethattheS decrease
tial condition (x,y,z) = (0.1,0.1,0.1) kpc, (v ,v ,v ) = n
x y z
(142,140,251) kms−1.Eachorbitwasintegratedforeighttimes towardstheboundarysowearecontentthatwehaveincludedthe
relevantmodes.
thelongestperiod,T .WesetN =6andcalculatedthecorre-
F max
Fortheboxorbitthetrueandrecoveredactionsare
spondingnumberofuniformly-spacedtimesamplesrequiredfrom
equation(21).Weensuredthatequations(19)and(20)weresatis- Jtrue =(336.39,137.78,237.96)kpckms−1
fiedforalltheincludedmodes.InFigs1and2weshowtheorbits
inthe(x,y)and(x,z)planes,thesamplingofthetoy-anglespace Jrecov =(336.85,137.26,238.17)kpckms−1,
andtheresultantactions.Wealsoshow,infaintred,theresultof andthetrueandrecoveredfrequenciesare
integratinginthebest-fittingtoypotential.Thisgivesusanideaof
theworkthatthegeneratingfunctionhastodotodeformthetoy Ωtrue =(39.752,46.409,73.814)Gyr−1
torusintothetargettorus. Ωrecov =(39.750,46.406,73.811)Gyr−1.
c 2014RAS,MNRAS000,1–12
(cid:13)
6 J.L.Sanders&J.Binney
loop InFig.4weshowtwocross-sectionsofn-spaceshowingthe
6 absolutevalueofthecomponentsofthegeneratingfunction.The
4 n3=0 2.50 two most significant modes are n = (2,−2,0), which causes
a mixing between the x motion and the y motion, and n =
2 2.25 (2,0,−2), which mixes the x and z motions. These modes are
2 0 required to distort the rectangular orbits of the triaxial harmonic
n
2.00
oscillatorintothoseboundedbysurfacesofconstantconfocalel-
2 1)
− 1.75 0. lipsoidalcoordinate.NotethattheSndecreasetowardsthebound-
4 + ariesasrequired.AlsothestructureofFig.3ismuchricherthan
−6 1.50 1−| thatofFig.4,signallingthatthegeneratingfunctionhasmanymore
− 6 4 2 0 2 4 6 ms significantterms.
− − −
n1 1.25 ck
6 p
k
4 n2=0 1.00 S/n 3.1 Accuracyofthemethod
(|
2 0.75 log Fig.5showserrorsinJ3(cid:48) andΩ(cid:48)3fortheboxorbitasafunctionof
n3 0 0.50 Nmax forvariouschoicesofthetotalintegrationtimeT.Wehave
linked N to N via equation (21). However we have not en-
2 max T
− suredthatequations(19)and(20)aresatisfiedforeachcase.The
0.25
4 weightofthepointsisproportionaltothelargestgapincoverage
−
fortheN modes.Weseethatingeneralalongerintegrationtime
6
− 6 4 2 0 2 4 6 providesamoreaccurateestimateoftheactionandparticularlythe
− − −
n1 frequency. We can understand this as a longer line segment pro-
vides a better measurement of the gradient for noisy data. From
Figure3.Cross-sectionsoftheSn asafunctionofnforthelooporbit. Fig.5weseethatwhenworkingwithhighNmaxitisnotsufficient
Inthetoppanelweshowthecross-sectionn3 = 0.Themostsignificant to satisfy equation (21). We must also satisfy equation (20) such
modeinthisplaneis( 1,2,0),whichcausesamixingbetweentheradial thatwehaveasufficientsamplingintoy-anglespacetoconstrain
−
motionandazimuthalmotion.Inthelowerplaneweshowthecross-section thesehighermodes.
n2=0,inwhichthemode(0,0,2)isthemostsignificant. ForT =2T equation(20)isnotsatisfiedforN (cid:62)4.For
F max
largeN andT =2T manymodeshaveinsufficientcoverage
max F
andtheresultsareverypoor.Fortheotherthreeintegrationtimes
box equation (20) is not satisfied for N (cid:62) 8. For T = 4T this
6 max F
resultsinanimmediatedeteriorationofthefrequencyrecoveryas
4 n3=0 2.50 wehaveincludedamodewithmax(n·θ)−min(n·θ)≈π/2.
2 2.25 ForT = 8TF andT = 12TF alackofcoveragehasnotaffected
the results apart from for T = 8T and N = 12 where the
2 0 F max
n 2.00 frequencyrecoveryispoorer.Themodewhichisnotwellcovered
2 1) isalsonotwellcoveredforN = 8butweonlyseetheeffects
−4 1.75 +0. ofthislackofcoveragewhenmwaxetrytoincludemoremodes.For
−−6 6 4 2 0 2 4 6 1.50 1ms−| dTes=pit1e2eTqFuabtiootnh(t2h0e)ancotitobneainngdsfaretiqsufieendcywhreecnoNvemryaxar(cid:62)ev8e.ryIngpoaord-
− − − n1 1.25 ck ticularthereisonemodeforwhichmax(n·θ)−min(n·θ)≈4.3.
6 p Itseemsthatthiscoverageissufficienttonotdegradetheresults.In
k
4 n2=0 1.00 (S/n| ccoiensclaunsdioanc,twiohnesnweeqlul,atwiohnil(s2t0w)hisensaittisifisendowtesarteicsfiovedertthheerferecqouveerny-
2 0.75 og deterioratesinsomecases,particularlythatofthefrequency.
l
3 0 Finally,wefindthatwhenwedoublethenumberoftimesam-
n 0.50
plesusedfortheexamplesshowninFig.5theresultschangesig-
2
− nificantlyonlywhenequation(20)isnotsatisfied.Therefore,we
0.25
4 concludethatprovidedwehavemoreequationsthanunknownsand
−
havesatisfiedequations(19)and(20)theactionsandfrequencyre-
6
− 6 4 2 0 2 4 6 coverywillbesatisfactory.
− − −
n1
Figure 4. As Fig. 3 for the box orbit. The most significant mode in the 3.2 Near-resonantorbit
planen3 = 0(top)is(2, 2,0),whichcausesamixingbetweenthex
motionandtheymotion.In−theplanen2=0themostsignificantmodeis To illustrate some of the points discussed we show results for
(2,0, 2). a near-resonant orbit. This orbit is a box orbit with the ini-
−
tial conditions (x,y,z) = (0.1,0.1,0.1) kpc, (v ,v ,v ) =
x y z
(142,150,216.5) kms−1.AgainweintegratefortimeT = 8TF
andsetN =6.TheresultsareshowninFig.6.Thefrequency
max
vectorofthisorbitisnearlyparallelton=(−4,0,2)sothecover-
ageofthismodeisverypoorandmax(n·θ)−min(n·θ)≈1.11
c 2014RAS,MNRAS000,1–12
(cid:13)
Actionsandanglesforintegratedorbits 7
T =2.0TF T =8.0TF 4 3
T =4.0TF T =12.0TF 3 2
2
6 1
1
c c
1| 4 kp 0 kp 0
− / /
J/J3,true003|−202 y −−−−−43215.0−2.5 0.0 2.5 5.0 z −−−−3215.0−2.5 0.0 2.5 5.0
g1 4 x/kpc x/kpc
o −
l
6 2.0 2.0
−
1| 0
− 1 1.5 1.5
ue −
3,tr −2 /π 1.0 /π 1.0
Ω 3 2 3
/ − θ θ
Ω03 −4 0.5 0.5
10|−5
log −76 0.00.0 0.5 1.0 1.5 2.0 0.00.0 0.5 1.0 1.5 2.0
− 2 4 6 8 10 12 θ1/π θ1/π
Nmax Toyaction Trueaction
450
400
aFsigaufruen5ct.ioEnrroofrNinmthaxe raencdovtheeretdotvaalliunetsegorfatJio3(cid:48)natinmdeΩ.W(cid:48)3efowrothrkewboitxhoevrbeint 1− 350
multiplesoftheperiod,TF,correspondingtothelowestfrequency.The ms 300
sizeofthepointsisproportionaltothelargestgapincoveragefortheN k 250
c
modes.Ingeneralalongerintegrationtimeprovidesmoreaccurateactions kp 200
andfrequencies.Whenattemptingtoconstrainhighermodesitisnecessary J/ 150
tointegratetheorbitforalongerperiodtoensurethatthesamplingintoy 100
angleissufficient. 50
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
t/Gyr
forn=(−4,0,2).However,thetrueandrecoveredactionsare
Figure6.Examplenear-resonantorbitinthetriaxialSta¨ckelpotential–the
Jtrue =(301.74,147.63,165.36)kpckms−1 toppanelsshowsinblacktheorbitintegratedinthetestpotential.Thisis
ashort-axislooporbitsocirculatesabouttheaxisz = 0.Infaintredwe
Jrecov =(300.69,147.66,165.89)kpckms−1, showtheinitialpointintegratedinthebest-fittingisochronepotential.In
themiddlepanelsweshowthetoyanglescalculatedateachtimesample.
andthetrueandrecoveredfrequenciesare
Inthebottompanelweshowthetoyactionsateachtime-stepasadotted
Ωtrue =(43.318,50.369,86.724)Gyr−1 line(blackforJ1,blueforJ2andredforJ3).Thesolidlinesshowthetrue
actionsandthearrowsmarktheestimatedactions.
Ωrecov =(43.386,50.371,86.777)Gyr−1.
Asseenbefore,poorcoverageinoneofthemodesisnotdetrimen- withM = 3.4×1010M andc = 0.7kpc;andthetriaxial
bulge
taltotheactionandfrequencyrecovery. (cid:12)
logarithmichalo
(cid:16) z2 (cid:17)
Φ (x,y,z)=v2 log C x2+C y2+C xy+ +r2
halo halo 1 2 3 q2 halo
4 APPLICATION z
(26)
Asabriefapplicationofthemethodoutlinedinthispaperwewill
inspecttheactiondiagramforarealistictriaxialGalacticpotential. withvhalo =121.7kms−1,C1 =0.99kpc−2,C2 =0.53kpc−2,
WetakethepotentialfromLaw&Majewski(2010).Thispotential C3 =0.11kpc−2,qz =1.36andrhalo =12kpc.
wasfoundtoproducethebestfittotheSagittariusstreamdata.This
potentialhasthreecomponents:adiscdefinedbytheMiyamoto-
4.1 Anexampleorbit
Nagaipotential
We inspect a single orbit in this potential in Fig. 7. The
−GM
Φdisc(x,y,z)= (cid:113) dis√c , (24) chosen orbit is a short-axis loop orbit with initial condi-
x2+y2+(a+ z2+b2)2 tion (x,y,z) = (14.69,1.80,0.12) kpc, (v ,v ,v ) =
x y z
(15.97,−128.90,44.68) kms−1. We use different, but overlap-
withMdisc =1011M(cid:12),a=6.5kpcandb=0.26kpc;aspheri- ping,8TF longsegmentsoftheorbitwithNT = 500tocalculate
calbulgedescribedbytheHernquistprofile the actions, angles and frequencies using N = 6. We ensure
max
−GM equation(20)issatisfiedforthesetimesamplings.Thisorbitliesin
Φbulge(r)= r+bculge, (25) thesurfaceofconstantenergyexploredinthenextsection.Wefind
c 2014RAS,MNRAS000,1–12
(cid:13)
8 J.L.Sanders&J.Binney
thattheactionandfrequencyare 15 2.0
J(cid:48) ≈ (160.18,2186.16,36.09)kpc kms−1 10 1.5
1.0
Ω(cid:48) ≈ (27.26,19.12,37.01)Gyr−1. c 5 c 0.5
p p
k 0 k 0.0
The error in the actions and frequencies can be estimated by the / /
spreadoftheestimatesfromeachsegment.Wefind y 5 z −0.5
− 1.0
10 −
∆J(cid:48) ≈ (0.07,0.08,0.03)kpckms−1, − −1.5
15 2.0
∆Ω(cid:48) ≈ (3×10−4,6×10−5,2×10−3)Gyr−1. − 8 0 8 − 8 0 8
− −
For each orbit segment we find θ(cid:48)(0) and these different values x/kpc x/kpc
should all lie along straight lines with gradients given by the de- 1− 0.15 1− 1.0
risivfieedd.frequencies.InFig.7weshowthattheconditioniswellsat- ms 0.10 Gyr 00..68
Usingdifferentorbitsegmentsisperhapstheonlywaytoes- pck 0.05 000 0.4
k 1 0.2
timatetheerrorinanactionorfrequencyfoundusingthepresent / 0.00 /
)i )i 0.0
mheerteh.oHdo.wIteisvesri,mapbleestttetormuseethcoodnsiesctuotiuvseeoorbrbititsesgegmmenentstsassewpearhaatevde J01h−0.05 Ω01h−00..42
− 0.10 −−
buytilaizilnarggethteimesetiimntaetrevdalg.eTnherisatcinagnfbuenacctihoinevteodfimndosatneifnfeitciatilvceolyndbiy- (J01 −−0.150.0 0.4 0.8 1.2 (Ω01 −−00..860.0 0.4 0.8 1.2
tionforasecondorbitintegration.Asimplechoiceistoincrease
t/Gyr t/Gyr
oneofthederivedanglecoordinatesbyπ/2.
2.0 2.0
4.2 Atypicalconstantenergysurface 1.5 1.5
Nowweturntoconstructingtheactiondiagramforthechosenpo- π π
/ 1.0 / 1.0
tential.Foragivenenergy(thatofaparticledroppedat18kpcon θ2 θ3
the intermediate axis) we launched particles at a series of points 0.5 0.5
linearly spaced between 0.2 and 18kpc along the potential’s in-
termediate axis with the velocity vector perpendicular to the axis 0.0 0.0
andinclinedatlinearlyspacedanglestothez-axisbetween0and 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0
π/21. We integrated each initial condition for ∼ 10Gyr saving θ1/π θ1/π
N = 1000 samples. For all orbits the energy was conserved
T
to one part in 106. We set N = 6 and ensured that equa- Figure7.AnexampleorbitintheLaw&Majewski(2010)potential.Itis
max
tions (19) and (20) were satisfied. If equation (19) was not satis- ashort-axislooporbitwithactionsJ(cid:48) ≈ (160,2186,36)kpckms−1.
fied,wehadundersampledtheorbit,sowetookafinersampling. In the top panel we show a 16TF long orbit segment in the (x,y) and
If equation (20) was not satisfied we did not have sufficient cov- (x,z)planes.InthecentraltwopanelsweshowthespreadinJ1(cid:48) andΩ(cid:48)1
erage,sowecontinuedintegratingforanother10Gyr,takingan- calculatedusing500time-samplesfroman8TF orbitsegmentlabelledby
itsinitialtimesample.Inthebottompanelweshowthecalculatedanglesat
other1000samples.Wethencalculatedtheactionsfromthetime
thesetimeswithblackdots.Wealsoshowtheanglesfoundusingθ(cid:48)(0)+
series.Fig.8showseachorbitasapointin3Daction-space2.We Ω(cid:48)ti withoneofthecalculatedfrequenciesandinitialanglesinsmaller
seethatthesurfaceofconstantenergyisatriangle-shapedplanein bluedots.
action-space.Thepointsarecolouredbasedontheirorbitclassifi-
cation.AnequivalentfigureforaSta¨ckelpotentialcanbefoundin
axis,theJ = 0,J = 0orbitisaclosedorbitinthe(x,y)plane
deZeeuw(1985). 1 3
andtheJ = 0,J = 0orbitisaclosedorbitinthe(y,z)plane.
Inatriaxialpotential,thelooporbitscanbedividedintotwo 1 2
Wenotethatneartheinterfacebetweenthedifferentorbitclasses
classes:theshort-axisloopsthatlooparoundtheshortaxis(inour
someregionsoftheplanearedepletedofpoints(ourchoiceofini-
casethez-axis)andthelong-axisloopsthatlooparoundthelong
tialsamplingcausesanincreaseddensityofpointsneartheedges
axis(thex-axis).Alongwiththeboxorbitsthesethreeclassesof
oftheplane).Alsothereissomeoverlapbetweenthedifferentorbit
orbitoccupydistinctregionsontheaction-spaceplaneofconstant
classesintheactionspace.Thesefeaturesareduetothepresence
energy. At each corner of the plane only one action is non-zero
ofresonantislandswithsurroundingchaoticorbitsattheinterface
andthecorrespondingorbitistheparentorbitofeachofthethree
of the regular orbit regions (see Section 5.3). For orbits near the
classes:theJ = 0,J = 0orbitisaradialorbitalongthelong
2 3
box/loop boundary it can take many orbital periods to correctly
identifytheorbitclass(Carpintero&Aguilar1998),andsomemay
1 NotethattheintermediateaxisofthehalomodelproposedbyLaw& bemisclassified.
Majewskiisactuallythez-axis.However,atsmallradii((cid:46) 18kpc)the
intermediateaxisofthefullpotentialisinthe(x,y)planeduetothedisc
contribution,andthez-axisistheshortaxis.
5 DISCUSSION
2 Toproduceacontinuousplaneinaction-spacewemustscalethe‘radial’
actionsofthelooporbits,J1,byafactorof2.J1 foralooporbitcorre- 5.1 Relationtopreviouswork
spondstoasingleoscillationfromminimumtomaximumcoordinateand
back,whilstforaboxorbitasingleoscillationcoverstheinterval0tomax- TheproblemaddressedheregoesbacktoBinney&Spergel(1982,
imumcoordinatefourtimes. 1984),whoFouriertransformedthetimeseriesx(ti)ofindividual
c 2014RAS,MNRAS000,1–12
(cid:13)
Actionsandanglesforintegratedorbits 9
theconstructionofanglecoordinates.Inboththesestudiesangle-
Short-axisloops Long-axisloops Boxes
actionvariableswereevaluatedalongnumericallycomputedorbits.
Thecoordinatesevaluatedwerenotthoseofatoypotentialbutof
atrialtorusthathadbeenpreviouslyconstructed:Warnock(1991)
wasrefiningtheFouriercoefficientsS whileKaasalainen&Bin-
n
3000 ney(1994)weresolvingforthe∂ S giventheS .Inboththese
i n n
J2400 studies,severalinitialconditionsfororbitintegrationwerechosen
30/ oneachtorustoovercometheproblemthatwithasingleshortinte-
k
p1800 grationaresonantorbityieldsahighlynon-uniformdistributionof
c
km1200 samplepointsonthetorus.Sincewedonothaveagoodrepresen-
−s tationofthetargettorusuntiltheequationshavebeensetupand
1600 solved,wecannottakeadvantageofthispossibility.
Warnock(1991)solvedforthediscreteFouriertransformsof
0
thenS ratherthanfortheS becausethematrixthatthenhasto
0 600 600 0 beinvenrtedisnearlydiagonalnwhenthetoyandtargettoriareclose
J20/kp1c2k0m018s−0021400300030002400J1801/0k01p2c0k0m s −1 tidoniftofhenereesnaptna,ocatenhdoefrittaoinsydhatanhrgedlestsoa.maScpihnleiceepvooeiunartustonpyirfooavnrmiddetsaaarmgnepetlationrlrgyiocreafgntoubylea-raqngugriitldee
space,wehavenotusedWarnock’stechnique.
5.2 PossibilityofusingSta¨ckeltori
3000
Wehaveusedcompletelydifferenttoypotentialsforeachclassof
24001 orbit,anditisnaturaltoaskwhetheritwouldnotbeadvantageous
−
s tousealwaysaSta¨ckelpotentialsincesuchapotentialhastoriof
1800m
k everytype.Wehavenotpursuedthisoptionfortworeasons.First,
c
p theactionsandanglesofSta¨ckelpotentialsrequiretheevaluation
1200k
J/03 of integrals whereas the potentials we have used yield algebraic
600 expressionsforanglesandactions.Secondly,andmorefundamen-
tally, when integrating an orbit that lies close to the box/loop in-
J10/3k0p20c401k08m0102s0−060010030002400J1208/0k0pc12k0m0s6−010 0 0 tttehorerfuaascc.etBi,oyintsuwsoiofnugtlhdpeobttaeerngnteoitanlo-srtrbtihivtaitahlasudtoptpheoenrstsuaormneleythgoaentoemthteeytprtyeoyaosftotthoruerustsaw,rgwitehet
are assured from the outset that this condition is satisfied. How-
Figure8.Twoprojectionsofasurfaceofconstantenergyinthe3Daction ever,thisrestsonourcorrectidentificationoftheorbittypefrom
spaceofthepotentialproposedbyLaw&Majewski(2010).Blackcircles
thetimeseries.AswesawwiththeLaw&Majewskipotential,in
showshort-axislooporbits,redcrossesshowlong-axislooporbitsandblue
somemarginalcasesitmaytakemanyorbitalperiodstocorrectly
trianglesshowboxorbits.
identifytheorbit.
coordinatesandassignedtoeachlineintheresultingspectrumap-
(cid:80) 5.3 Resonancesandchaos
propriateintegersn sothatωtcouldbeidentifiedwith n Ω t.
j j j j
Oncethisidentificationhadbeensuccessfullyaccomplished,Ω t Wehavefocusedhereonorbitsthatarenon-resonantmembersof
j
couldbereplacedwithθ toyieldtheorbit’sanglerepresentation. themajororbitalfamilies.Inrealgalacticpotentialsoneencounters
j
Thisapproachisinferiortothatintroducedhereinseveralrespects: orbitsthatareeitherresonantlytrappedorchaotic(e.g.§3.7Binney
(i)Whereasthegeneratingfunctionisascalar,astar’slocationis &Tremaine2008).Chaoticorbitscanbethoughtofassequences
describedbyavector,soitiswastefultoconstructtheanglerep- ofsectionsofresonantlytrappedorbits,sothesetwotypesoforbit
resentationsofallthreecoordinatesratherthantheanglerepresen- raisesimilarissues.
tationofthegeneratingfunction:Binney&Spergel(1984)failed In a generic integrable potential, the frequencies Ω depend
i
totakeadvantageofthestrongrestrictionsontorithatarisefrom on the actions, so on some tori a resonant condition n·Ω = 0
angle-actioncoordinatesbeingcanonical.(ii)Itisnotstraightfor- issatisfied.Consequently,individualorbitsontheseresonanttori
wardtomeasurecorrectlythecomplexamplitudesAfromthedis- donotcovertheentiretorussincetheconditionn·θ = constant
creteFouriertransformofatimeseriessuchasx(t )becausethe constrainstheanglevariables.Thislackofcoveragemakesitim-
i
requiredamplitudewillingeneralnotlieatoneofthediscretefre- possibletodeterminesomeoftheFouriercoefficientsS .
n
quenciessampled.(iii)Whenanorbitisnear-resonantthereisoften Whenthepotentialisstrictlyintegrable,orbitsontorithatare
dangerousambiguityintheintegersn thatshouldbeassignedto adjacent to a resonant torus completely cover their tori although
j
aparticularline.Withthepresenttechniqueweworkfromtheout- theytakealongtimetodoso.Inagenericpotential,however,such
setwithperiodicfunctionsandtheirFourierseriessotheissueof orbitsmoveoveraseriesoftoriwithoutcoveringanyofthem,as
how frequencies fall on a discrete grid does not arise. Moreover, theylibratearoundthestrictlyresonantorbit.Consequently,these
theassignmentofintegersn toFouriertermsisunambiguous. orbitshavesomeofthecharacteristicsofastrictlyresonanttorus.
j
The method described here has significant overlap with the When the present technique is used on a resonantly trapped or-
workofWarnock(1991)ontheconstructionofmagneticcoordi- bit,thegeneratingfunctionwillmapthetoytorusintoacloseap-
natesandtherelatedmethodofKaasalainen&Binney(1994)for proximationtothestrictlyresonanttorus,soinanN-bodymodel
c 2014RAS,MNRAS000,1–12
(cid:13)
10 J.L.Sanders&J.Binney
thedensityofstarsonthistoruswillseemtobelargerthanitre-
allyis.Hencewiththepresenttechnique,resonantlytrappedorbits
will give rise to apparent crowding in action space that is analo-
goustothesignatureofresonanceswhenparticlesaremappedinto
2.0
frequencyspacebydeterminingorbitalfrequenciesbyFourierde-
compositionofcoordinates(Dumas&Laskar1993):whenthera-
tiosΩ /Ω andΩ /Ω areusedtoplaceorbitsinfrequency-ratio
2 1 3 1
space,theexistenceofresonantlytrappedorbitsleadstoacrowding 1.8
ofpointsalongthestraightlinesassociatedwithcertainresonance
conditionsn·θ=constant(Binney&Tremaine2008,§3.7.3(b)).
Chaoticorbitscanbeconsideredasmovingthroughaseries 1
Ω
of quasi-periodic orbits. Therefore the recovered actions and fre- / 1.6
3
quenciesfromourmethodwillbeafunctionofthetotalintegration Ω
time.Weseethattheregionoftheconstantenergysurfaceoccu-
piedbytheboxorbitsinFig.8hasconsiderablecrowdingandthe
1.4
regular grid of initial conditions is not visible. This is indicative
ofchaoticorbitswhichhavebeenallocatedverydifferentactions
fromoneinitialconditiontothenext.
In Fig. 9 we perform the same procedure as outlined in
1.2
§3.7.3(b)ofBinney&Tremaine(2008)toinspecttheratiooffre-
quenciesplaneofalogarithmicpotential.Weusethepotential
0.9 1.0 1.1 1.2 1.3 1.4 1.5
1 (cid:16) y2 z2 (cid:17)
Φ(x)= 2log x2+ q2 + q2 +rc2 , (27) Ω2/Ω1
y z
Figure 9. Frequency ratios in the triaxial logarithmic potential extracted
withq =0.9,q =0.7andr2 =0.1.Wedropaseriesoftestpar-
y z c fromorbitaltimeseriesusingthemethodpresentedinthispaper.Eachpoint
ticlesonthesurfaceΦ(x) = 0.5regularlyspacedinthespherical corresponds to an initial condition for a particle dropped on the surface
polarcoordinatesφandcosθ,andintegrateeachinitialcondition Φ(x)=0.5.
foratimeT = 200extractingN = 2048samples.Wethenuse
T
ourmethodtofindthecorrespondingorbitalfrequenciesandplot
their ratio in Fig. 9. As noted in Binney & Tremaine (2008) the An important application is to the analysis of N-body sim-
top-rightcornerofthisplaneshowstheregularspacingoftheini- ulations. A single N-body snapshot consists of 3D positions and
tialconditionswhilstthelower-leftcornershowsamoreirregular velocitiesfor∼ 109 particles.Lettingthesimulationevolvefora
distribution with no evidence of the regular grid of initial condi- fewtimestepsproducesanothersnapshotwithacompletelydiffer-
tionsusedtoproduceit.Also,wefindthatthereareoverdensities entsetof109 positionsandvelocities.Thustheparticles’phase-
alonglinescorrespondingtoresonances.Ourplotisverysimilarto spacecoordinatesconstituteahighlydegenerateandnon-compact
thatshowninBinney&Tremaine(2008).However,thestructure representationofthesimulation.Effectiveanalysisofthesimula-
oftheirregularbottom-leftregiondiffers.Thisistobeexpectedas tionshouldstartbycondensingthecoordinatesintoasmallersetof
itistheseorbitswhichareirregular,andhowoneassignsregular numbers.Thiscanbedonebyreplacingthe6NT numbers(xi,vi)
propertiestothemdependsonthemethodemployed. withjustthreenumbersJi andplottingeachparticleasapointin
3Dactionspace.Thesimulationthenbecomesadensityofparticles
ina3Dspace.Thisrepresentationwillgreatlyfacilitatethecom-
parisonofdifferentN-bodymodels.Alsoitmayprovepossibleto
6 CONCLUSIONS
findgoodfitstothestardensityintermsofanalyticfunctions,as
We have presented a method for finding actions, frequencies and Pontzen&Governato(2013)havedonefornumericaldark-matter
angles from numerically integrated orbits in a general potential3. haloes and appears to be possible for the Galactic discs (Binney
ThemethodreliesonestimatingtheFouriercomponentsofthegen- 2012b;Binneyetal.2014).Wehopetoreportonanapplicationof
eratingfunctionthatmapsatoytorusintothetorusonwhichthe thismethodtoanN-bodysimulationsoon.
computed orbit lies by solving systems of linear algebraic equa- ItshouldbenotedthatitisnotadvisabletotaketheNT time
tions.Thismethodenablesonetodeterminetheangle-actioncoor- samples of a given orbit directly from the simulation. Rather at
dinates(θ,J)ofagivenphase-spacepoint(x,v).Ithasnumerous some time t the potential should be computed on a spatial grid
possibleapplicationsinastronomy. (e.g.Magorrian2007),andtheequationsofmotioninthispoten-
Oursisthefirstmethodpresentedintheliteratureforfinding tialshouldbeintegratedforNT timestepsstartingfromthephase-
the actions in a general triaxial potential. Triaxiality is an essen- space location of each particle at time t. These integrations in a
tial ingredient of dark-matter distributions, and a realistic Galac- fixedpotentiallendthemselvestomassiveparallelization,forex-
ticmodelwhichshouldincludenon-axisymmetricfeaturessuchas ampleonaGraphicalProcessorUnit(GPU)soitshouldbepossi-
the bar, and the potentially triaxial halo. This method is a neces- bletocomputeangle-actioncoordinatesforverylargenumbersof
sary first step towards constructing distribution functions, f(J), particles.
forthesemorecomplexGalacticcomponents. Here we discussed time-reversible triaxial potentials. In this
casewecandetermineapriorithephasesofthetermsinthegen-
eratingfunction.Rotationofthefigureofthepotentialdestroysthe
3 We will make the code developed for this paper available at time-reversibilityoftheHamiltonianandwelosetheabilitytoset
https://github.com/jlsanders/genfunc. thephasesapriori.Intheworstcase,theSninequation(5)become
c 2014RAS,MNRAS000,1–12
(cid:13)