Recovery of the internal orbital structure of galaxies
aa r X i v : . [ a s t r o - ph ] J a n Mon. Not. R. Astron. Soc. , 1–33 (2007) Printed 2 December 2018 (MN L A TEX style file v2.2)
Recovery of the internal orbital structure of galaxies
G. van de Ven , , ⋆ † , P. T. de Zeeuw , , R. C. E. van den Bosch Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands Department of Astrophysical Sciences, Peyton Hall, Princeton, NJ 08544, USA Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA European Southern Observatory, D-85748 Garching bei M¨unchen
Accepted 0000 Month 00. Received 0000 Month 00; in original 0000 Month 00
ABSTRACT
We construct axisymmetric and triaxial galaxy models with a phase-space distribution func-tion that depends on linear combinations of the three exact integrals of motion for a separa-ble potential. These Abel models, first introduced by Dejonghe & Laurent and subsequentlyextended by Mathieu & Dejonghe, are the axisymmetric and triaxial generalisations of thewell-known spherical Osipkov–Merritt models. We show that the density and higher ordervelocity moments, as well as the line-of-sight velocity distribution (LOSVD) of these mod-els can be calculated efficiently and that they capture much of the rich internal dynamicsof early-type galaxies. We build a triaxial and oblate axisymmetric galaxy model with pro-jected kinematics that mimic the two-dimensional kinematic observations that are obtainedwith integral-field spectrographs such as
SAURON . We fit the simulated observations with ax-isymmetric and triaxial dynamical models constructed with our numerical implementation ofSchwarzschild’s orbit-superposition method. We find that Schwarzschild’s method is able torecover the internal dynamics and three-integral distribution function of realistic models ofearly-type galaxies.
Key words: stellar dynamics – celestial mechanics – galaxies: elliptical and lenticular, cD –galaxies: kinematics and dynamics – galaxies: structure
The equilibrium state of a collisionless stellar system such as anelliptical or lenticular galaxy is completely described by its distri-bution function (DF) in the six-dimensional phase space of posi-tions and velocities. The recovery of the DF from observations isdifficult, as for galaxies other than our own, we can usually onlymeasure the projected surface brightness and the line-of-sight ve-locity distribution (LOSVD) of the integrated light as a functionof position on the plane of the sky. Moreover, we generally do notknow the intrinsic shape of the galaxy, nor the viewing direction, orthe contribution to the gravitational potential provided by a supermassive central black hole and/or an extended halo of dark matter.By Jeans (1915) theorem, the DF is a function of the isolating inte-grals of motion admitted by the potential, but it is not evident howto take advantage of this property other than for the limiting case ofspherical systems. Orbits in axisymmetric geometry have two exactintegrals of motion, the energy E and the angular momentum com-ponent L z parallel to the symmetry z -axis, but the third effectiveor non-classical integral I obeyed by all regular orbits is gener-ally not known in closed form. In stationary triaxial geometry E ⋆ Hubble Fellow † E-mail: [email protected] is conserved, but regular orbits now have two additional effectiveintegrals of motion, I and I , which are not known explicitly.Schwarzschild (1979, 1982) devised a numerical methodwhich sidesteps our ignorance about the non-classical integrals ofmotion. It allows for an arbitrary gravitational potential, which mayinclude contributions from dark components, integrates the equa-tions of motion for a representative library of orbits, computesthe density distribution of each orbit, and then determines the or-bital weights such that the combined orbital densities reproducethe density of the system. The best-fitting orbital weights representthe DF (cf. Vandervoort 1984). Pfenniger (1984) and Richstone& Tremaine (1984) included kinematic moments in this method,and Rix et al. (1997) showed how to include observed LOSVDs.A number of groups have developed independent numerical imple-mentations of Schwarzschild’s method for axisymmetric geometrywhich fit the projected surface brightness and line-of-sight velocitydistributions of early-type galaxies in detail (van der Marel et al.1998; Cretton et al. 1999; Gebhardt et al. 2000, Valluri, Merritt &Emsellem 2004; Thomas et al. 2004; Cappellari et al. 2006). Ap-plications include the determination of central black hole masses(see also van der Marel et al. 1997; Cretton & van den Bosch1999; Verolme et al. 2002; Cappellari et al. 2002; Gebhardt et al.2003; Copin, Cretton & Emsellem 2004), accurate global dynam-ical mass-to-light ratios (Cappellari et al. 2006), as well as darkmatter profiles as a function of radius (Cretton, Rix & de Zeeuw c (cid:13) G. van de Ven et al. E and L z only (van der Marel et al. 1998; Cretton et al. 1999; Verolme& de Zeeuw 2002; Valluri et al. 2004; Cretton & Emsellem 2004;Thomas et al. 2004; Krajnovi´c et al. 2005).One could construct a numerical galaxy model withSchwarzschild’s method itself, compute the observables, and thenuse these as input for the code and determine how well it recoversthe input model. This is useful, but does not provide a fully inde-pendent test of the software. An alternative is to consider the specialfamily of models with gravitational potential of St¨ackel form, forwhich all three integrals of motion are exact and known explicitly.These separable potentials have a core rather than a central cusp, sothe corresponding models cannot include a central black hole, andare inadequate for describing galactic nuclei. However, they can beconstructed for a large range of axis ratios (Statler 1987), and theirobserved kinematic properties are as rich as those seen in the mainbody of early-type galaxies (Statler 1991, 1994; Arnold, de Zeeuw& Hunter 1994).A small number of analytic DFs have been constructed fortriaxial separable models. The ‘thin-orbit’ models (Hunter & deZeeuw 1992) have the maximum possible streaming motions, buttheir DF contains delta functions, and they are therefore not par-ticularly useful for a test of general-purpose numerical machinery.Dejonghe & Laurent (1991, hereafter DL91) constructed separa-ble triaxial models in which the DF depends on a single param-eter S = E + wI + uI , which is a linear combination of thethree exact integrals E , I and I admitted by these potentials, andis quadratic in the velocity components. For a given radial den-sity profile, the DF follows by simple inversion of an Abel integralequation. These so-called Abel models have no net mean streamingmotions, and are the axisymmetric and triaxial generalisations ofthe well-known spherical Osipkov–Merritt models (Osipkov 1979;Merritt 1985), for which the observables can be calculated easily(Carollo, de Zeeuw & van der Marel 1995). Mathieu & Dejonghe(1999, hereafter MD99) generalised the results of DL91 by includ-ing two families of DF components with net internal mean motionsaround the long and the short axis, respectively, and compared theresulting models with observations of Centaurus A. Although theAbel character of the non-rotating components is no longer con-served, the expressions for the velocity moments in these more gen-eral models can still be evaluated in a straightforward way. Whenthe entire DF depends on the same single variable S the famousellipsoidal hypothesis (e.g., Eddington 1915; Chandrasekhar 1940)applies, so that self-consistency is only possible in the sphericalcase (Eddington 1916; Camm 1941). This does not hold for Abelmodels with a DF that is a sum of components for which the vari-able S has different values of the parameters w and u . Such multi- component Abel models can provide (nearly) self-consistent mod-els with a large variety of shapes and dynamics.Here, we show that for Abel models, in addition to the velocitymoments, the full LOSVD can be calculated in a simple way. Next,we construct axisymmetric and triaxial Abel models to test our nu-merical implementation of Schwarzschild’s method. We assume aconvenient form for the gravitational potential, and construct theDF that reproduces a realistic surface brightness distribution. Wecompute the LOSVDs of the models and derive two-dimensionalmaps of the resulting kinematics. We show that, despite the simpleform of the DF, these models display the large variety of featuresobserved in early-type galaxies with integral-field spectrographssuch as SAURON (Emsellem et al. 2004). By fitting axisymmetricand triaxial three-integral Schwarzschild models to the simulatedobservables we find that Schwarzschild’s method is able to recoverthe internal dynamics and three-integral DF of early-type galaxies.In this paper we fix the mass-to-light ratio and viewing directionto those of the Abel models, while in our companion paper vdB07we investigate how well these global parameters can be determinedby Schwarzschild’s method, along with a full description of ournumerical implementation in triaxial geometry.This paper is organised as follows. In Section 2 we summarisethe properties of the triaxial Abel models of DL91 and MD99and present the intrinsic velocity moments in a form which fa-cilitates their numerical implementation. We describe the conver-sion to observables in Section 3, including the computation of theLOSVD. In Section 4 we construct a specific triaxial galaxy modeland in Section 5 we fit the simulated observables with our triax-ial Schwarzschild models to investigate how well the intrinsic mo-ments and three-integral DF are recovered. In Section 6 we con-sider Abel models in the axisymmetric limit and construct a three-integral oblate galaxy model to test our axisymmetric implementa-tion of Schwarzschild’s method. We summarise our conclusions inSection 7. In Appendix A, we describe the simpler Abel models forthe elliptic disc, large distance and spherical limit, and link them tothe classical Osipkov–Merritt solutions for spheres. Readers whoare mainly interested in the tests of the Schwarzschild method mayskip Sections 2 – 4 and 6.1 – 6.3.
The triaxial Abel models introduced by DL91 have gravitationalpotentials of St¨ackel form, for which the equations of motion sep-arate in confocal ellipsoidal coordinates. We briefly describe thesepotentials, and refer for further details to de Zeeuw (1985a). Wethen make a specific choice for the DF, for which the velocity mo-ments simplify.
We define confocal ellipsoidal coordinates ( λ, µ, ν ) as the threeroots for τ of x τ + α + y τ + β + z τ + γ = 1 , (2.1)with ( x, y, z ) the usual Cartesian coordinates, and with constants α, β and γ such that − γ ≤ ν ≤ − β ≤ µ ≤ − α ≤ λ . From theinverse relations x = ( λ + α )( µ + α )( ν + α )( α − β )( α − γ ) , (2.2) c (cid:13) , 1–33 ecovery orbital structure and similarly for y and z by cyclic permutation of α → β → γ → α , it follows that a combination ( λ, µ, ν ) generally corre-sponds to eight different points ( ± x, ± y, ± z ). In these coordinates,the St¨ackel potentials have the following form (Weinacht 1924) V S ( λ, µ, ν ) = U ( λ )( λ − µ )( λ − ν ) + U ( µ )( µ − ν )( µ − λ )+ U ( ν )( ν − λ )( ν − µ ) , (2.3)where U ( τ ) is an arbitrary smooth function ( τ = λ, µ, ν ) . Theright-hand side of eq. (2.3) can be recognised as the second orderdivided difference of U ( τ ) . Henceforth, we denote it with the cus-tomary expression U [ λ, µ, ν ] , which is symmetric in its arguments(see Hunter & de Zeeuw 1992, eqs 2.1–2.3, 2.13 and 2.14). Addi-tion of a linear function of τ to U ( τ ) does not change V S .The density ρ S that corresponds to V S can be found from Pois-son’s equation πGρ S ( λ, µ, ν ) = ∇ V S ( λ, µ, ν ) , (2.4)or alternatively by application of Kuzmin’s (1973) formula (see deZeeuw 1985b). This formula shows that, once we have chosen theconfocal coordinate system and the density along the short axis,the mass model is fixed everywhere by the requirement of separa-bility . For centrally concentrated mass models, V S has the x -axisas long-axis and the z -axis as short-axis. In most cases this is alsotrue for the associated density (de Zeeuw, Peletier & Franx 1986). The Hamilton-Jacobi equation separates in ( λ, µ, ν ) for the poten-tials (2.3), so that every orbit has three exact integrals of motion(cf. de Zeeuw & Lynden-Bell 1985) E = ` v x + v y + v z ´ + U [ λ, µ, ν ] ,I = T L y + L z + ( α − β ) v x +( α − β ) x U [ λ, µ, ν, − α ] , (2.5) I = L x + (1 − T ) L y + ( γ − β ) v z +( γ − β ) z U [ λ, µ, ν, − γ ] , where v x , v y and v z are the velocity components in the Carte-sian coordinate system, and L x = yv z − zv y , the componentof the angular momentum vector parallel to the x -axis. The othertwo components, L y and L z , follow by cyclic permutation of x → y → z → x and v x → v y → v z → v x . Furthermore, T is a triaxiality parameter defined as T = ( β − α ) / ( γ − α ) , (2.6)and U [ λ, µ, ν, σ ] is the third-order divided difference of U ( τ ) . Allmodels for which U ′′′ ( τ ) > have a similar orbital structure andsupport four families of regular orbits: boxes with no net rotation,inner and outer long-axis tubes with net rotation around the x -axis,and short-axis tubes with net rotation around the z -axis (Kuzmin1973; de Zeeuw 1985a; Hunter & de Zeeuw 1992). A third method for the calculation of the density is to use πGρ S = H [ λ, λ, µ, µ, ν, ν ] , where the fifth-order divided difference is of the func-tion H ( τ ) = 4 a ( τ ) U ′ ( τ ) − a ′ ( τ ) U ( τ ) with a ( τ ) = ( τ + α )( τ + β )( τ + γ ) and U ( τ ) defines the potential as in eq. (2.3). This result was ob-tained by Hunter in 1989 (priv. comm.) and by Mathieu & Dejonghe (1996).Similar expressions exist for the related families of potential-density pairsintroduced in de Zeeuw & Pfenniger (1988). According to Jeans (1915) theorem the phase-space distribu-tion function (DF) is a function f ( E, I , I ) of the isolating inte-grals of motion (cf. Lynden-Bell 1962; Binney 1982). The velocitymoments of the DF are defined as µ lmn ( λ, µ, ν ) = ZZZ v lλ v mµ v nν f ( E, I , I ) d v λ d v µ d v ν , (2.7)where l , m and n are non-negative integers, and v λ , v µ and v ν arethe velocity components in the confocal ellipsoidal coordinate sys-tem. Many of the velocity moments vanish due to the symmetry ofthe orbits in these coordinates. The zeroth-order velocity momentis the mass density that corresponds to the DF ρ ⋆ ( λ, µ, ν ) = µ ( λ, µ, ν ) . (2.8)In self-consistent models, ρ ⋆ must equal ρ S given in eq. (2.4), themass density that is related to the potential V S by Poisson’s equa-tion. Following DL91, we choose the DF to be a function of the threeintegrals of motion E , I and I as given in eq. (2.5) through onevariable f ( E, I , I ) = f ( S ) , with S = − E + w I + u I , (2.9)and w and u are two parameters . This choice for the DF isequivalent to the celebrated ellipsoidal hypothesis (e.g., Edding-ton 1915; Chandrasekhar 1940). Self-consistency is only possiblein the spherical case (Eddington 1916; Camm 1941). On the otherhand, these DFs can produce realistic (luminous) mass densities ρ ⋆ ,which differ from the (total) mass density ρ S , as in galaxies withdark matter (see also § w and u .)DL91 and MD99 divided the DF into three types of compo-nents. The non-rotating (NR) type is made of box orbits and tubeorbits with both senses of rotation populated equally. The two rotat-ing types, LR and SR, consist of tube orbits, and have net rotationaround either the long axis or the short axis. Due to the choice (2.9) of the DF, the general expression (2.7) forthe velocity moments can be simplified, as shown by DL91 forthe non-rotating components and by MD99 for the rotating compo-nents. We recast their expressions into a different form to facilitatethe numerical implementation. The resulting velocity moments aregiven by µ lmn ( λ, µ, ν ) = s l + m + n +3 H l +1 µν H m +1 νλ H n +1 λµ × S max Z S min T lmn [ S top ( λ, µ, ν ) − S ] ( l + m + n +1) / f ( S ) d S, (2.10)and set to zero at positions for which S max ≤ S min . The terms H µν , H νλ and H λµ in the square root in front of the integral aredefined as H στ = 1 + ( σ + α )( τ + α ) γ − α w + ( σ + γ )( τ + γ ) α − γ u, (2.11) In contrast with DL91 and MD99, we choose V S ≤ and E ≤ ,consistent with e.g. de Zeeuw (1985a).c (cid:13) , 1–33 G. van de Ven et al. [t]
Figure 1.
The limiting value S lim of the variable S = − E + w I + u I as function of the parameters w and u . The physical region is bounded bythe relations (2.12), indicated by the thick solid lines. The dashed curvesdivide this region into three parts, each with a different expression for S lim .The relations for the separatrices L and L are given in eq. (2.13). with σ, τ = λ, µ, ν . Orbits are confined to the region of space forwhich all three terms are non-negative. In general, this conditionwill not be satisfied for all points, so that the Abel componentshave finite extent. From the requirement that at least the origin ( λ, µ, ν ) = ( − α, − β, − γ ) should be included, we find the fol-lowing limits on w and uw ≥ − β − α and u ≤ γ − β . (2.12)The factor T lmn in the integrand as well as the upper limit S max of the integral are different for each of the three Abel componenttypes NR, LR and SR, and are discussed in §§ S min has to be at least as large as thesmallest value possible for the variable S . This limiting value S lim depends on the choice of the DF parameters w and u in (2.9), as isshown in Fig. 1 (cf. Fig. 7 of DL91). The boundaries follow from(2.12) and the separatrices L and L are given by L : w = U [ − α, u − γ, − γ ]( β − α ) (2.13) L : w = u − ( γ − α ) u . At a given position ( λ, µ, ν ) , orbits with different values of the in-tegrals of motion E , I and I , and hence different values of S ,can contribute to the integral (2.10). The restriction to bound or-bits ( E ≤ ) together with the requirement that v λ , v µ and v ν allthree have to be non-negative determines the part of the integralspace that is accessible by orbits that go through ( λ, µ, ν ) . An ex-ample of the resulting tetrahedron in the ( E, I , I ) -space is shownin Fig. 2 (cf. Fig. 1 of MD99). The largest possible value of S is [t] Figure 2.
The tetrahedron shows all accessible points in integral space ( E, I , I ) for a given position ( λ, µ, ν ) . The tetrahedron is bounded bythe planes for which v λ = 0 , v µ = 0 , v ν = 0 and E = 0 , respectively.The two shaded planes, which are given by v λ = v µ = 0 at λ = µ = − α and v µ = v ν = 0 at µ = ν = − β , divide the tetrahedron into theparts corresponding to the four general orbit families in a triaxial separablepotential: box (B) orbits, inner (I) and outer (O) long-axis tube orbits andshort-axis (S) tube orbits. given by the top of this tetrahedron S top ( λ, µ, ν ) = − U [ λ, µ, ν ] − w ( λ + α )( µ + α )( ν + α ) γ − α U [ λ, µ, ν, − α ] − u ( λ + γ )( µ + γ )( ν + γ ) α − γ U [ λ, µ, ν, − γ ] , (2.14)which is thus a function of the position ( λ, µ, ν ) . At the origin S top ( − α, − β, − γ ) = U [ − α, − β, − γ ] , which is the central valueof the potential V S . In what follows, we normalise V S by setting U [ − α, − β, − γ ] = − , so that ≤ S top ≤ . Since the non-rotating component type can exist everywhere in theaccessible integral space (the tetrahedron in Fig. 2), we simply havethat S max = S top ( λ, µ, ν ) . Spatially the NR components are thusbounded by the surface S top ( λ, µ, ν ) = S min .The factor T lmn follows from the cross section of the S -planewithin the tetrahedron and can be written in compact form as T NR lmn = B ( l +12 , m +12 , n +12 ) , (2.15)where B is the beta function of three variables . Since T NR lmn is in-dependent of S it can be taken out of the integral (cf. eq. [3.10] of The beta function of k variables in terms of the complete gamma function Γ is defined as B ( β , . . . , β k ) = Γ( β ) · · · Γ( β k ) / Γ( β + · · · + β k ) .c (cid:13) , 1–33 ecovery orbital structure DL91), which then becomes of Abel form. Unfortunately, the inver-sion of eq. (2.10) for any chosen moment µ lmn ( λ, µ, ν ) , includingthe case l = m = n = 0 , is generally impossible, as the left-handside is a function of three variables, while the DF depends on onlyone variable, S . The density ρ ⋆ specified along any given curve willdefine a different f ( S ) . A case of particular interest is to choose thedensity along the short axis to be ρ ⋆ (0 , , z ) = ρ S (0 , , z ) . Thisdefines a unique f ( S ) , and hence gives ρ ⋆ everywhere. Kuzmin’sformula applied to ρ S (0 , , z ) similarly defines the density ρ S ev-erywhere. For single Abel DF components these will not be thesame, except in the spherical limit (see Appendix A3).Since the orbits have no net rotation, the velocity moments µ NR lmn are only non-zero when l , m and n are all three even, andvanish in all other cases. The long-axis rotating component type only exists in the part ofthe integral space that is accessible by the (inner and outer) long-axis tube orbits. Within the tetrahedron for all orbits this is theregion for which v ν ≥ at ν = − β . It follows that S max = S top ( λ, µ, − β ) ≤ S top ( λ, µ, ν ) .The term T lmn follows from the cross section of the S -planewithin the tetrahedron and with the above boundary plane v ν = 0 at ν = − β . Without any further constraint this results in zero netrotation, because each orbit with positive rotation around the longaxis with v ν > , is balanced by an orbit with opposite direction ofrotation with v ν < . Therefore, we restrict to orbits with v ν ≥ ,resulting in maximum streaming around the long axis for each LRcomponent. This reduces the accessible integral space, and thusalso the term T lmn , by a factor of two, so that the latter becomes T LR lmn = 2 ( − ( l + m ) / q a l +10 b m +10 M LR0 ( s + 1)( s − . . . ( s + 1 − ( l + m )) , (2.16)with s = l + m + n , the parameters a and b defined as a = ( λ + β ) H µν [ S top ( λ, µ, − β ) − S ]( λ − ν ) H µ ( − β ) [ S top ( λ, µ, ν ) − S ] , (2.17) b = ( µ + β ) H νλ [ S top ( λ, µ, − β ) − S ]( µ − ν ) H ( − β ) λ [ S top ( λ, µ, ν ) − S ] , which for S ≤ S max = S top ( λ, µ, − β ) are non-negative, and M LR0 = ( M ( s, l , m ; a , b , π ) , a ≤ b , M ( s, m , l ; b , a , π ) , a > b . (2.18)The function M ( s, i, j ; a, b, φ ) is defined in Appendix B, wherewe evaluate it in terms of elementary functions (odd s ) and ellipticintegrals (even s ).The LR components have maximum streaming around thelong axis, but the motion parallel to the intermediate axis and shortaxis cancels. As a result, the velocity moments µ LR lmn vanish when l or m are odd . Multiplying µ LR lmn with ( − n results in maximumstreaming in the opposite direction. By choosing different weightsfor both senses of rotation, we can control the direction and theamount of long-axis streaming motion for each LR component. Since l + m is even, the factor ( − ( l + m ) / in eq. (2.16) is always real. The short-axis component type reaches the part of integral spaceaccessible by the short-axis tube orbits. Within the tetrahedron forall orbits this is the region for which v µ ≥ both at µ = − β and µ = − α (Fig. 2). The latter requirement is equivalent to I ≥ . Inthis case, S max = S top ( λ, − α, ν ) ≤ S top ( λ, µ, ν ) .The form of the term T lmn depends on the cross section ofthe S -plane within the tetrahedron and with the above two bound-ary planes. In case each SR component has maximum streamingaround the short axis ( v µ ≥ ), it is given by T SR lmn = 2 ( − ( l + n ) / P i =1 q a l +1 i c n +1 i M SR i ( s + 1)( s − . . . ( s + 1 − ( l + n )) . (2.19)The parameters a and c follow from a and b defined in (2.17)by interchanging µ ↔ ν , and in turn a and c follow from a and c by interchanging α ↔ β . For the terms M SR i we have twopossibilities, I and II, M SRI = M ( s, l , n ; a I , c I , θ I ) , a I ≤ c I , M ( s, n , l ; c I , a I , π ) −M ( s, n , l ; c I , a I , π − θ I ) , a I > c I , (2.20) M SRII = M ( s, l , n ; a II , c II , π ) −M ( s, l , n ; a II , c II , θ II ) , a II ≤ c II , M ( s, n , l ; c II , a II , π − θ II ) , a II > c II , (2.21)where M is given in Appendix B, and θ I and θ II follow from tan θ I = c II ( a I − a II ) a II ( c II − c I ) , tan θ II = c I ( a II − a I ) a I ( c I − c II ) . (2.22)For the assignment of the labels I and II , we discriminate betweenfour cases a ≤ a , c ≥ c : I → , II → ,a ≥ a , c ≤ c : I → , II → , (2.23) a ≤ a , c ≤ c : I → , θ I = π/ , C SR2 = 0 ,a ≥ a , c ≥ c : I → , θ I = π/ , C SR1 = 0 . The SR components have maximum streaming around the shortaxis, so that the velocity moments µ SR lmn vanish when l or n areodd. Multiplying µ SR lmn with ( − m results in SR components withmaximum streaming around the short axis in the opposite direction. Until now, we have chosen the Abel DF to be a function of a sin-gle variable S = − E + wI + uI , and we have separated it inthree component types, NR, LR and SR, but we have not madeany assumption about the form of the DF (apart from the obviousrequirement that it has to be non-negative everywhere and that itdecreases to zero at large radii). Following MD99, we choose theDF to be a linear combination of basis functions of the form f δ ( S ) = „ S − S min − S min « δ , (2.24)which, like the velocity moments (2.10), are non-vanishing as longas S lim ≤ S min ≤ S ≤ S max ≤ S top ≤ . The exponent δ is a(non-negative) constant.Once the St¨ackel potential (2.3) is known by defining the func-tion U ( τ ) , we can use the above relations ( § c (cid:13) , 1–33 G. van de Ven et al. f δ ( S ) the velocity moments (2.10) for the NR, SR and LR compo-nents in an efficient way, where at most the integral over S has tobe evaluated numerically. For the NR components this integral caneven be evaluated explicitly, resulting in µ NR lmn,δ ( λ, µ, ν ) = s [2( S max − S min )] l + m + n +3 H l +1 µν H m +1 νλ H n +1 λµ × „ S max − S min − S min « δ B ( l +12 , m +12 , n +12 , δ +1) , (2.25)where S max = S top ( λ, µ, ν ) (cf. eq. 2.14).Each DF component and corresponding velocity momentsthus depend on the choice of the DF parameters w , u and δ , thetype of component, and for the rotating components (LR and SR),they also depend on the sense of rotation around the axis of sym-metry. By summing a series of DF basis functions over w , u and δ , one might even expect to cover a large fraction of all physicalDFs. Due to the different values of w and u , such a sum of DFcomponents is no longer a function of the same, single variable S , so that the ellipsoidal hypothesis does not apply. Consequently,it becomes possible to construct (nearly) self-consistent dynamicalmodels, with the (combined) luminous mass density ρ ⋆ equal (orclose) to the mass density ρ S associated to the potential. We describe how the intrinsic velocity moments can be convertedto projected velocity moments on the plane of the sky. Alterna-tively, these line-of-sight velocity moments follow as moments ofthe LOSVD, which we show can be calculated in a straightfor-ward way for Abel models. Parameterising the LOSVD as a Gauss-Hermite series, we obtain the observable quantities: the surfacebrightness, the mean line-of-sight velocity V , velocity dispersion σ , and higher-order Gauss-Hermite moments h , h , . . . In order to calculate line-of-sight velocity moments, we introducea new Cartesian coordinate system ( x ′′ , y ′′ , z ′′ ) , with x ′′ and y ′′ inthe plane of the sky and the z ′′ -axis along the line-of-sight. Choos-ing the x ′′ -axis in the ( x, y ) -plane of the intrinsic coordinate sys-tem (cf. de Zeeuw & Franx 1989 and their Fig. 2), the transforma-tion between both coordinate systems is known once two viewingangles, the polar angle ϑ and azimuthal angle ϕ , are specified. Theintrinsic z -axis projects onto the y ′′ -axis, which for an axisymmet-ric galaxy model aligns with the short axis of the projected surfacedensity Σ . However, for a triaxial galaxy model the y ′′ -axis in gen-eral lies at an angle ψ with respect to the short axis of Σ . Thismisalignment ψ can be expressed in terms of the viewing angles ϑ and ϕ and the triaxiality parameter T (defined in eq. 2.6) as follows(cf. eq. B9 of Franx 1988) tan 2 ψ = − T sin 2 ϕ cos ϑ sin ϑ − T ` cos ϕ − sin ϕ cos ϑ ´ (3.1)with sin 2 ψ sin 2 ϕ cos ϑ ≤ and − π/ ≤ ψ ≤ π/ . A ro-tation over ψ transforms the coordinate system ( x ′′ , y ′′ , z ′′ ) to ( x ′ , y ′ , z ′ ) , with the x ′ -axis and y ′ -axis aligned with respectivelythe major and minor axis of Σ , whereas z ′ = z ′′ is along the line-of-sight.The expressions in § ( λ, µ, ν ) . The conver-sion to line-of-sight quantities can be done by four successive ma-trix transformations. First, we obtain the velocity components in thefirst octant of the intrinsic Cartesian coordinate system ( x, y, z ) via Q , of which the first element is given by (cf. DL91) Q = sign( λ + α ) s ( µ + α )( ν + α )( λ + β )( λ + γ )( α − β )( α − γ )( λ − µ )( λ − ν ) , (3.2)and the other elements follow horizontally by cyclic permutationof λ → µ → ν → λ and vertically by cyclic permutation of α → β → γ → α . The second matrix uses the symmetries of theorbits to compute the appropriate signs of the intrinsic Cartesianvelocities in the other octants. The result depends on whether ornot the orbit has a definite sense of rotation in one of the confocalcoordinates. For the three types of Abel components this results inthe following matricesNR : S = diag[sgn( x ) , sgn( y ) , sgn( z )] LR : S = diag[sgn( xyz ) , sgn( z ) , sgn( y )] (3.3)SR : S = diag[sgn( y ) , sgn( x ) , sgn( xyz )] Finally, the conversion from the intrinsic to the observer’s Carte-sian velocities involves the same projection and rotation as for thecoordinates. We represent these two coordinate transformations re-spectively by the projection matrix P = − sin ϕ cos ϕ − cos ϑ cos ϕ − cos ϑ sin ϕ sin ϑ sin ϑ cos ϕ sin ϑ sin ϕ cos ϑ , (3.4)and the rotation matrix R = cos ψ − sin ψ ψ cos ψ
00 0 1 . (3.5)In this way, we arrive at the following relation v x ′ v y ′ v z ′ = M v λ v µ v ν , with M ≡ RPSQ , (3.6)where the full transformation matrix M is thus a function of ( λ, µ, ν ) , the constants ( α, β, γ ) and the viewing angles ( ϑ, ϕ, ψ ) . We can now write each velocity moment in the observer’s Carte-sian coordinate system ( x ′ , y ′ , z ′ ) as a linear combination of thevelocity moments in the confocal ellipsoidal coordinate system µ ijk ( x ′ , y ′ , z ′ ) = X l,m,n c l,m,n ( s ) µ lmn ( λ, µ, ν ) , (3.7)with s = i + j + k = l + m + n . The coefficients c l,m,n ( s ) areproducts of elements of the transformation matrix M in eq. (3.6).They can be obtained with the following recursive algorithm c l,m,n ( s ) = c , , ( s ) c l − ,m,n ( s − , if l > ,c , , ( s ) c ,m − ,n ( s − , if m > ,c , , ( s ) c , ,n − ( s − , if n > , (3.8)with the first order expressions given by c , , ( s ) = M e s , c , , ( s ) = M e s , c , , ( s ) = M e s , (3.9) c (cid:13)000
00 0 1 . (3.5)In this way, we arrive at the following relation v x ′ v y ′ v z ′ = M v λ v µ v ν , with M ≡ RPSQ , (3.6)where the full transformation matrix M is thus a function of ( λ, µ, ν ) , the constants ( α, β, γ ) and the viewing angles ( ϑ, ϕ, ψ ) . We can now write each velocity moment in the observer’s Carte-sian coordinate system ( x ′ , y ′ , z ′ ) as a linear combination of thevelocity moments in the confocal ellipsoidal coordinate system µ ijk ( x ′ , y ′ , z ′ ) = X l,m,n c l,m,n ( s ) µ lmn ( λ, µ, ν ) , (3.7)with s = i + j + k = l + m + n . The coefficients c l,m,n ( s ) areproducts of elements of the transformation matrix M in eq. (3.6).They can be obtained with the following recursive algorithm c l,m,n ( s ) = c , , ( s ) c l − ,m,n ( s − , if l > ,c , , ( s ) c ,m − ,n ( s − , if m > ,c , , ( s ) c , ,n − ( s − , if n > , (3.8)with the first order expressions given by c , , ( s ) = M e s , c , , ( s ) = M e s , c , , ( s ) = M e s , (3.9) c (cid:13)000 , 1–33 ecovery orbital structure and c , , = 1 . The index e s is the s th element of the vector e = [3 , .., , , .., , , .., , where the number of integers 3 ( )is equal to the value of the velocity moment index k , and similarly j and i .The line-of-sight velocity moments now follow from (numer-ical) integration of µ k along the line-of-sight µ k ( x ′ , y ′ ) = Z ∞−∞ µ k ( x ′ , y ′ , z ′ ) d z ′ , (3.10)which are thus functions of position on the sky plane. Using the definition of the intrinsic velocity moments of the DF(eq. 2.7) and rearranging the sequence of integration, we rewriteeq. (3.10) for the line-of-sight velocity moments as µ k ( x ′ , y ′ ) = Z ∞−∞ v kz ′ L ( x ′ , y ′ , v z ′ ) d v z ′ , (3.11)where we have introduced the LOSVD L ( x ′ , y ′ , v z ′ ) = ZZZ f ( E, I , I ) d v x ′ d v y ′ d z ′ . (3.12)Although the integral over z ′ in general can only be evaluated nu-merically, we show that for the choice (2.9) of the DF, the doubleintegral over the velocities can be simplified significantly.Our analysis generalises the results for the well-known spher-ical Osipkov–Merritt models. We describe the spherical limit to-gether with the elliptic disc and large distance limit in Appendix A,while we present axisymmetric Abel models in § Substituting the expressions (2.5) for the integrals of motion in S = − E + wI + uI , we obtain S = S top ( λ, µ, ν ) − ` H µν v λ + H νλ v µ + H λµ v ν ´ , (3.13)where the expression for S top ( λ, µ, ν ) is given by eq. (2.14) andthe terms H µν , H νλ and H λµ are defined in eq. (2.11). Defining X = H µν S top ( λ, µ, ν ) − S ] v λ , (3.14)and similarly Y and Z by cyclic permutation of λ → µ → ν → λ ,we can write the expression for S as X + Y + Z = 1 . (3.15)For a given position ( λ, µ, ν ) , each value of S thus defines the sur-face of the unit sphere in the variables ( X, Y, Z ) . In these variables,we can write the integral of the DF over velocities, i.e., the stellarmass density, as ρ ⋆ = ZZ Z f ( S ) d v x ′ d v y ′ d v z ′ , = S max Z S min s S top − S ] H µν H νλ H λµ f ( S ) d X d Y d Z X + Y + Z =1 d S. (3.16)This is the same expression as for the zeroth-order velocity momentof the DF, µ , in eq. (2.10), where T is equal to the integralbetween square brackets.The matrix M in eq. (3.6) provides the conversion from thevelocity components in the confocal ellipsoidal coordinate system, ( v λ , v µ , v ν ) , to those in the observer’s Cartesian coordinate system, ( v x ′ , v y ′ , v z ′ ) . Hence, for a given line-of-sight velocity v z ′ , we find e X + e Y + e Z = v z ′ / g ( S ) . (3.17)The coefficients e , e and e are defined as h e = p H νλ H λµ M ,h e = p H λµ H µν M , (3.18) h e = p H µν H νλ M , and normalised with respect to h given by h = H νλ H λµ M + H λµ H µν M + H µν H νλ M . (3.19)These coefficients are functions of the position ( λ, µ, ν ) , the con-stants ( α, β, γ ) and the viewing angles ( ϑ, ϕ, ψ ) through the com-ponents of the matrix M , and also depend on the DF parameters w and u through the terms H λµ , H µν and H νλ . It follows that g ( S ) = h s S top − S ] H λµ H µν H νλ , (3.20)which is a function of the variable S .We thus find that each combination of values of S and v z ′ results in the cross section of the surface of the unit sphere ineq. (3.15) with the plane in eq. (3.17), i.e., a circle, in the vari-ables ( X, Y, Z ) . We rotate the latter coordinate system such thatthe normal vector ( e , e , e ) of the plane of the circle coincideswith the Z ′ -axis of the system given by XYZ = cos Φ sin Φ cos Θ sin Φ sin Θ − sin Φ cos Φ cos Θ cos Φ sin Θ0 − sin Θ cos Θ
1A 0@ X ′ Y ′ Z ′ . (3.21)where the rotation angles Φ and Θ follow from tan Φ = e e , tan Θ = p e + e e . (3.22)In these coordinates the circle is conveniently parameterised as X ′ = p − Z ′ cos ξ ′ , Y ′ = p − Z ′ sin ξ ′ , (3.23)where Z ′ = v z ′ /g ( S ) . We can now rewrite the integral betweensquare brackets in eq. (3.16) as ZZ ˛˛˛˛ ∂ R ∂ξ ′ ∧ ∂ R ∂Z ′ ˛˛˛˛ d ξ ′ d Z ′ = 1 g ( S ) ZZ d ξ ′ d v z ′ , (3.24)where the vector R = ( X, Y, Z ) and ∧ indicates the cross product.The integral over ξ ′ is the length of the part of the circle, ∆ ξ ′ ,for which the corresponding integral space is accessible by orbits,and hence is in general a function of S and v z ′ and differs for thedifferent types of Abel components as we show below.Inserting eq. (3.24) in eq. (3.16), we obtain ρ ⋆ = 1 h S max Z S min g ( S ) Z − g ( S ) f ( S ) ∆ ξ ′ ( v z ′ , S ) d v z ′ d S, = 1 h g ( S min ) Z − g ( S min ) S up ( v z ′ ) Z S min f ( S ) ∆ ξ ′ ( v z ′ , S ) d S d v z ′ , (3.25)where after changing the order of integration in the last step, theupper limit of S is given by S up = min[ G ( v z ′ ) , S max ] , with G ( v z ′ ) = S top ( λ, µ, ν ) − H µν H νλ H λµ v z ′ h . (3.26) c (cid:13) , 1–33 G. van de Ven et al.
Comparing the first line of eq. (3.16) with the second line ofeq. (3.25), we see that the choice of the Abel DF, f ( E, I , I ) = f ( S ) , indeed reduces the triple integration (3.12) for the LOSVDto a double integral: L ( x ′ , y ′ , v z ′ ) = ∞ Z −∞ h S up ( v z ′ ) Z S min f ( S ) ∆ ξ ′ ( v z ′ , S ) d S d z ′ , (3.27)and vanishes when | v z ′ | exceeds the ’terminal velocity’ v t = g ( S min ) . The expressions for h and S up follow from eqs (3.19)and (3.26), whereas S max and ∆ ξ ′ are different for each of thethree Abel component types and are considered next. As for the intrinsic moments in § S max = S top ( λ, µ, ν ) , and, since the full in-tegral space is accessible, ∆ ξ ′ NR = 2 π , independent of S and v z ′ .In the case of a basis function f δ ( S ) as defined in eq. (2.24),the integral over S can be evaluated explicitly resulting in L NR δ = 2 π ( δ +1)(1 − S min ) δ ∞ Z −∞ h [ G ( v z ′ ) − S min ] δ +1 d z ′ . (3.28) The integral space accessible by the (inner and outer) long-axis tubeorbits is given by v ν ≥ at ν = − β , so that immediately S max = S top ( λ, µ, − β ) , whereas the calculation of ∆ ξ ′ LR is more complex.Since eq. (3.15) must also hold at ν = − β , we find that forLR components, within the unit sphere in the variables ( X, Y, Z ) ,the space is restricted to that within the elliptic cylinder given by X a + Y b = 1 , (3.29)where a and b are defined in eq. (2.17). In the rotated coordi-nate system ( X ′ , Y ′ , Z ′ ) defined in eq. (3.21), at height Z ′ = v z ′ /g ( S ) , the elliptic cylinder results in an ellipse given by d X ′ + d Y ′ + d X ′ Y ′ + d X ′ + d Y ′ + d = 0 (3.30)for ≤ ξ ′ ≤ π and with coefficients d = ( b e + a e ) ,d = e ( b e + a e ) ,d = 2 e e e ( b − a ) , (3.31) d = 2 e e ( e + e ) ( b − a ) ( v z ′ /g ) ,d = 2 e ( e + e ) ( b e + a e ) ( v z ′ /g ) ,d = ( e + e ) ˆ ( b e + a e ) ( v z ′ /g ) − a b ˜ . Because all the above relations only involve the squared values of X , Y and Z , they are independent of the sign of the correspondingvelocities v λ , v µ and v ν (cf. eq. 3.14), which results in zero netrotation. For the LR components, to obtain net rotation around thelong axis, we simply limit the range of v ν values, e.g., requiring Z ≥ , results in maximum streaming motion around the longaxis. This restricts the space in ( X ′ , Y ′ , Z ′ ) , at given height Z ′ = v z ′ /g ( S ) , to one side of the line − ( e + e ) Y ′ + e ( v z ′ /g ) = 0 . (3.32)The restriction to the opposite side of the line inverts the rotation around the long axis. By choosing different weights for both sensesof rotation, we can control the direction and the amount of long-axis streaming motion.For given values of S and v z ′ , the integral space covered bythe LR components is thus the part of the circle in eq. (3.23) thatfalls within the ellipse in eq. (3.30) and that is on the correct sideof the line in eq. (3.32) (see also Fig. 3). The length ∆ ξ ′ LR of thispart thus ranges from zero to a maximum of π when the circleis completely inside the ellipse and on the correct side of the line.To compute this length, we determine the points where the circle(possibly) intersects the ellipse and the line. Substituting the circleparameterisation of eq. (3.23) in the expression for the ellipse ineq. (3.30), we find that the intersections with the ellipse are the(real) zero points of the following fourth order polynomial in u ≡ tan( ξ ′ / d − d R ′ + d R ′ ) u + 2 R ′ ( d − d R ′ ) u + 2 ˆ d + (2 d − d ) R ′ ˜ u + 2 R ′ ( d + d R ′ ) u + ( d + d R ′ + d R ′ ) = 0 , (3.33)where we have introduced R ′ ≡ p − ( v z ′ /g ) . The intersectionswith the line result in the following two solutions u ± = ( e + e ) R ′ ± ˆ − e − ( v z ′ /g ) ˜ e ( v z ′ /g ) , (3.34)for | v z ′ | ≤ g ( S )(1 − e ) , otherwise the line is outside the circle.We thus (numerically) find up to six real zero points u i andcorresponding angles ξ ′ i = 2 arctan( u i ) , sorted from low to high.For the set {− π, ξ ′ , ξ ′ , . . . , π } , we compute the lengths of the se-quential intervals on the circle for which the corresponding valuesfall within the ellipse and on the correct side of the line. This can bechecked by inserting a value from the corresponding interval, e.g.the central value, in eq. (3.23) and substituting the resulting X ′ and Y ′ into eqs (3.30) and (3.32). If the left-hand side is negative (posi-tive), the interval is inside (outside) the ellipse, and (for Z ≥ ) onthe wrong (correct) side of the line. Finally, the sum of the resultinginterval lengths provides ∆ ξ ′ LR . The short-axis tube orbits are restricted to the region of integralspace for which v µ ≥ both at µ = − β and µ = − α , and hence S max = S top ( λ, − α, ν ) . For the calculation of ∆ ξ ′ SR we have thatthe space within the unit sphere in ( X, Y, Z ) is now restricted tothe part that falls within both the elliptic cylinders X a + Z c = 1 and X a + Z c = 1 . (3.35)As in § a and c follow from a and b defined in (2.17) by interchanging µ ↔ ν , and in turn a and c follow from a and c by interchanging α ↔ β .Both elliptic cylinders result in ellipses in the Z ′ -plane, as ineq. (3.30) for the LR components, but now with coefficients d = c i e ,d = c i e e + a i ( e + e ) ,d = 2 c i e e e (3.36) d = 2 c i e e ( e + e ) ( v z ′ /g ) ,d = 2 e ( e + e ) ˆ c i e − a i ( e + e ) ˜ ( v z ′ /g ) ,d = ( e + e ) ˆ ( c i e + a i e ) ( v z ′ /g ) − a i c i ˜ . c (cid:13) , 1–33 ecovery orbital structure Figure 3.
The part of integral space accessible by long-axis (left panel) and short-axis (right panel) tube orbits, for given values of the Abel DF variable S and the line-of-sight velocity v z ′ . This corresponds to the thick part of the circle, which is restricted to be within the dashed ellipse(s) and below [above] thedashed line for long-axis [short-axis] rotating components with maximum streaming (see text for details). The length of the thick part of the circle equals ∆ ξ ′ in the expression (3.27) of the LOSVD. for i = 1 , respectively. The zero points of the correspondingfourth order polynomials in eq. (3.33) are again the intersectionswith the circle in eq. (3.23).Net rotation around the short axis follows by limiting the rangeof v µ values, e.g., Y ≥ yields maximum streaming, which re-stricts the accessible integral space to one side of the line − e X ′ + e e Y ′ + e ( e + e ) ( v z ′ /g ) = 0 . (3.37)The two solutions of the intersection with the circle are u ± = − e e R ′ ± ( e + e ) ˆ − e − ( v z ′ /g ) ˜ e R ′ + ( e + e ) e ( v z ′ /g ) , (3.38)for | v z ′ | ≤ g ( S )(1 − e ) , otherwise the line is outside the circle.The combination of all the (real) zero points provides the (or-dered) set {− π, ξ ′ , ξ ′ , . . . , π } , with at most ten intersections ξ ′ i with the circle given in eq. (3.23). We compute the lengths of thecircle intervals for which the enclosed values fall within both el-lipses and on the correct side of the line. This means, for which thecorresponding X ′ and Y ′ values substituted in eq. (3.30) result in anegative left-hand side for both pairs of a i and b i , and (for Y ≥ )in a positive left-hand side of eq. (3.37). Finally, ∆ ξ ′ SR is the sumof the resulting interval lengths. When considering the LR type of components we make no dis-tinction between inner and outer long-axis tube orbits because theyhave similar dynamical properties. Similarly, the non-rotating boxorbits are part of the NR type of components and are not consideredseparately. Nevertheless, if we are interested in the specific contri-bution of these orbit families to the LOSVD, this can be achievedby a straightforward extension of the above analysis.As can be seen from Fig. 2, the inner and outer long-axis tube orbits are separated by the plane I = 0 , or equivalently the regionsfor which v λ ≥ at λ = − α and v µ ≥ at µ = − α , respec-tively. This is in addition to the restriction v ν ≥ at ν = − β forboth long-axis tube orbits. For the inner long-axis tube orbits thisimplies that S max = S top ( λ, µ, − β ) . The space within the unitsphere in ( X, Y, Z ) is now restricted to the part that falls within theintersection of the elliptic cylinders in eq. (3.29) and Y b + Z c = 1 , (3.39)where b and c follow from a and b defined in (2.17) by inter-changing ν ↔ λ and β ↔ α . In the Z ′ -plane, these two ellipticcylinders result in ellipses as in eq. (3.30), with coefficients respec-tively given in eq. (3.31) and d = c e ,d = c e e + b ( e + e ) ,d = − c e e e (3.40) d = − c e e ( e + e ) ( v z ′ /g ) ,d = 2 e ( e + e ) ˆ c e − b ( e + e ) ˜ ( v z ′ /g ) ,d = ( e + e ) ˆ ( c e + b e ) ( v z ′ /g ) − b c ˜ . As before, ∆ ξ ′ follows from the combination of the real zeropoints of the corresponding fourth order polynomials in eq. (3.33),and of eq. (3.32) in the case of maximum streaming aroundthe long axis. For the outer long-axis tube orbits, S max =min[ S top ( λ, µ, − β ) , S top ( λ, − α, ν )] . The two elliptic cylindersare the one in eq. (3.29) and the second in eq. (3.35), with the co-efficients of the corresponding ellipses in the Z ′ -plane are given ineq. (3.31) and eq. (3.36) ( i = 2 ), respectively.The part of integral space accessible by box orbits is the regionfor which both v µ ≥ at µ = − β and v λ ≥ at λ = − α (Fig. 2).Therefore, S max = S top ( λ, µ, ν ) , and the two elliptic cylinders arethe first in eq. (3.35) and the one in eq. (3.39). The coefficients of c (cid:13) , 1–33 G. van de Ven et al. the corresponding ellipses in the Z ′ -plane are respectively those ineq. (3.36) ( i = 1 ) and in eq. (3.40). We have seen that the line-of-sight velocity moments µ k ( x ′ , y ′ ) can be derived either via line-of-sight integration of the intrinsic ve-locity moments (eq. 3.10) or as moments of the LOSVD (eq. 3.11).The lowest order line-of-sight velocity moments µ , µ and µ pro-vide the surface mass density Σ , the mean line-of-sight velocity V and dispersion σ by Σ = µ , V = µ µ , and σ = µ µ − µ µ , (3.41)all three as a function of ( x ′ , y ′ ) . Whereas Σ , V and σ can be mea-sured routinely, determinations of the higher order moments ( µ , µ , . . . ) are more complicated. Spectroscopic observations of theintegrated light of galaxies provide the LOSVD as function of po-sition on the sky plane. Unfortunately, the wings of the LOSVDbecome quickly dominated by the noise in the observations, andsince the higher order moments significantly depend on the wings,their measurements can become very uncertain. Instead of thesetrue higher-order moments, one often uses the Gauss-Hermite mo-ments ( h , h , . . . ), which are much less sensitive to the wings ofthe LOSVD (van der Marel & Franx 1993; Gerhard 1993).There is no simple (analytic) relation between the true mo-ments and the Gauss-Hermite moments, including the lower ordermoments Σ GH , V GH and σ GH (but see eq. 18 of van der Marel &Franx 1993 for approximate relations to lowest order in h and h ).Nevertheless, we have shown that for Abel models the full LOSVDcan be computed in a efficient way from eq. (3.27), so that by fittinga Gauss-Hermite series to the resulting LOSVD, we can derive theGauss-Hermite moments accurately, all as function of ( x ′ , y ′ ) .Still, the calculation of the line-of-sight velocity momentsthrough the intrinsic moments is useful, e.g., in case of investi-gating a range of viewing directions. The intrinsic moments haveto be computed once, after which only a (numerical) integrationalong the line-of-sight is needed for each viewing direction. Thisis (much) faster than calculating the LOSVD separately at each di-rection. The higher order true moments can even be used to (nu-merically) determine the Gauss-Hermite moments. One way is tofind the Gauss-Hermite LOSVD of which the true moments best-fit those from the Abel model. However, in practise this direct fit-ting of the true moments has several (numerical) problems. Becauseit is a non-linear minimisation problem, the convergence can takelong and may result in a local instead of the global best-fit solution,possibly resulting in Gauss-Hermite moments that are significantlydifferent from their true values. If, instead, we first (re)constructthe LOSVD from the true moments by means of an Edgeworthexpansion (see Appendix C) and then fit a Gauss-Hermite series,the Gauss-Hermite moments can be calculated accurately and effi-ciently. Evidently, once the viewing direction is known, it is morestraightforward to compute the full LOSVD to derive the (higher-order) Gauss-Hermite moments.When we construct a galaxy model consisting of multipleAbel DF components ( § , each multiplied witha constant weight, and then parameterise the resulting combinedLOSVD as a Gauss-Hermite series. Because the mass included ineach DF component is different, in order to obtain the mass frac-tions per DF component, we multiply the latter weights with themass of the corresponding DF component divided by the total (lu-minous) mass. To change the sense of rotation of a rotating DFcomponent (LR or SR), the corresponding observables do not haveto be recomputed, as a change in the sign of the odd velocity mo-ments is sufficient. The surface brightness follows upon integration of the luminos-ity density along the line-of-sight. The luminosity density in turnis related to the mass density ρ ⋆ via the stellar mass-to-light ra-tio M ⋆ /L . With ρ ⋆ the zeroth-order velocity moment of the DF(eq. 2.8), the surface brightness follows as SB( x ′ , y ′ ) = Z ∞−∞ ( M ⋆ /L ) − µ ( x ′ , y ′ , z ′ ) d z ′ . (3.42)In the special case when ( M ⋆ /L ) does not change (e.g., due tovariation in the underlying stellar populations) with position, wecan take it out of the integral and SB = Σ / ( M ⋆ /L ) , where Σ isthe surface mass density defined in eq. (3.41).In addition to the luminous matter, a galaxy may also containdark matter. While in the outer parts of late-type galaxies the pres-ence of dark matter, as predicted by the cold dark matter paradigmfor galaxy formation (e.g., Kauffmann & van den Bosch 2002),was demonstrated convincingly already more than two decades ago(e.g., van Albada et al. 1985), the proof in the outer parts of early-type galaxies remains uncertain, mainly due to a lack of kinematicconstraints. As a consequence, in the outer parts of galaxies, com-monly a simple functional form for the dark matter distributionis assumed, often the universal profile from the CDM paradigm(Navarro, Frenk & White 1997).The dark matter distribution in the inner parts of galaxies isprobably even more poorly understood (e.g., Primack 2004). Forthis reason, in current dynamical studies of the central parts ofearly-type galaxies, it is commonly assumed that both ( M ⋆ /L ) andthe dark matter fraction are constant, i.e., mass follows light. In thiscase the surface brightness also follows from SB = Σ S / ( M/L ) ,where ( M/L ) is the (constant) total mass-to-light ratio and Σ S thesurface mass density, which after deprojection yields ρ S , the massdensity related to the potential V S via Poisson’s equation (2.4). Incase of a St¨ackel potential (2.3), Σ S (and hence the surface bright-ness) has concentric isodensity contours that show no twist (e.g.,Franx 1988). After choosing a St¨ackel potential, we investigate the shape of thedensity generated by the Abel DF components, and use these com-ponents to construct a triaxial galaxy model with three integrals ofmotion. Or, in case the LOSVD is not readily accessible, the true line-of-sightvelocity moments, which are also linear functions of the DF.c (cid:13) , 1–33 ecovery orbital structure There are various choices for the potential that provide useful testmodels for comparison with the kinematics of triaxial ellipticalgalaxies (e.g., Arnold et al. 1994). One option is to consider the so-called perfect ellipsoid, for which Statler (1987) already computednumerical Schwarzschild models and Hunter & de Zeeuw (1992)investigated the maximum streaming thin orbit models. It has adensity distribution stratified on similar concentric ellipsoids, butthe potential function U ( τ ) contains elliptic integrals, which slowsdown numerical calculations. An alternative is to consider the setof models introduced by de Zeeuw & Pfenniger (1988), which havenearly ellipsoidal density figures, and have a potential and densitythat are evaluated easily and swiftly. They are defined by the choice: U ( τ ) = − GM √ τ ( τ + β ) , (4.1)so that the triaxial St¨ackel potential has the elegant form V S ( λ, µ, ν ) = − GM “ √ λµ + √ µν + √ νλ − β ” ( √ λ + √ µ )( √ µ + √ ν )( √ ν + √ λ ) , (4.2)where we set GM = √− γ + √− α so that V S ( − α, − β, − γ ) = − in the centre. In the oblate axisymmetric limit this potentialis that of the Kuzmin-Kutuzov (1962) models of Dejonghe & deZeeuw (1988), and in the spherical limit it reduces to H´enon’s(1959) isochrone. For all these models, V S = U [ − α, − β, τ ] along the short z -axis is identical to the isochrone potential − GM/ ( √ τ + √− α ) . We therefore refer to models with U ( τ ) ofthe form (4.1) as isochrone models. Since the potential falls of as /r at large radii, all these models have finite total mass.The expressions for the integrals of motion are given in (2.5),where U [ λ, µ, ν ] = V S and the third order divided difference U [ λ, µ, ν, σ ] is given by the symmetric expression U [ λ, µ, ν, σ ] = − GM − V S ( √ λ + √ µ + √ ν + √ σ )( √ λ + √ σ )( √ µ + √ σ )( √ ν + √ σ ) . (4.3)These triaxial isochrone models have the convenient property thatthe expressions for the potential and the integrals of motion containonly elementary functions of the (confocal ellipsoidal) coordinatesand have no singularities.The same is true for the associated mass density ρ S , of whichthe expression is given in Appendix C of de Zeeuw & Pfenniger(1988), and a contour plot of ρ S in the ( x, z ) -plane is shown intheir Fig. 2. These authors also derive the axis ratios of ρ S in thecentre (their eq. C7) and at large radii (their eq. C11), in terms ofthe axis ratios ζ and ξ of the confocal ellipsoidal coordinate system,defined as ζ = ( − β ) / ( − α ) , ξ = ( − γ ) / ( − α ) . (4.4)Although ρ S becomes slightly rounder at larger radii, its axis ratiosremain smaller than unity (for ξ < ζ < ) because at large radii ρ S ∼ /r in all directions. Characteristic values for the axis ratioscan be obtained from the (normalised) moments of inertia along theprincipal axes of the density, a = R x ρ ( x, ,
0) d x R ρ ( x, ,
0) d x , (4.5)where the intermediate and short semi-axis length, b and c , of Substituting eq. (4.2) shows that U [ λ, µ, ν, σ ] is in fact fully symmetric: U [ λ,µ,ν,σ ] GM = √ λµν + √ µνσ + √ νσλ + √ σλµ − β ( √ λ + √ µ + √ ν + √ σ )( √ λ + √ µ )( √ µ + √ ν )( √ ν + √ λ )( √ λ + √ σ )( √ µ + √ σ )( √ ν + √ σ ) the inertia ellipsoid follow from the long semi-axis length a byreplacing x with y and z , and at the same time ρ ( x, , with ρ (0 , y, and ρ (0 , , z ) , respectively. Taking for example ζ = 0 . and ξ = 0 . , the semi-axis lengths of the inertia ellipsoid resultin the characteristic axis ratios b S /a S = 0 . and c S /a S = 0 . for the density ρ S . The contours of the projected density are nearlyelliptic with slowly varying axis ratios.For triaxial mass models with a St¨ackel potential V S , deZeeuw, Peletier & Franx (1986) have shown that the correspondingintrinsic mass density ρ S cannot fall off more rapidly than /r ,except along the short z -axis. All models in which ρ S falls off lessrapidly than /r become round at large radii. When ρ S ∼ /r ,as is the case for, e.g., the above isochrone potential and the per-fect ellipsoid (e.g., de Zeeuw 1985a), the model remains triaxial atlarge radii. Moreover, mass models containing a linear combinationof different St¨ackel potentials are possible as long as the associatedconfocal ellipsoidal coordinate systems share the same foci (e.g., deZeeuw & Pfenniger 1988; Batsleer & Dejonghe 1994). This showsthat, although we choose here a (single-component) isochrone po-tential, our method is capable of providing Abel models for a largerange of St¨ackel potentials, with a similarly large range of shapesof the corresponding mass model. The same holds true for the lu-minous mass density, which we consider next. Whereas the shape of the (total) mass density ρ S is fixed by thechoice of the potential V S , and ζ and ξ (eq. 4.4), the shape of the(luminous) mass density ρ ⋆ , which is the zeroth order velocity mo-ment of the DF (eq. 2.8), also depends on the DF parameters w , u and δ , and the type of component. For ζ = 0 . and ξ = 0 . , weshow in Fig. 4 for non-rotating DF components the characteristic(eq. 4.5) axis ratios of the corresponding density, as function of w and u . We have set δ = 1 , but the axis ratios depend only weaklyon it, with ρ ⋆ becoming slightly flatter for increasing δ . The thickcontours are drawn at the levels that correspond to the values ofthe characteristic axis ratios of ρ S , respectively b S /a S = 0 . , c S /b S = 0 . and c S /a S = 0 . . These values are independentof w and u (as well as the other DF parameters).While the intermediate-over-long axis ratio b/a increases withincreasing w , its value is only weakly dependent of u . By contrast,the short-over-intermediate axis ratio c/b mainly increases with in-creasing u . The short-over-long axis ratio c/a is the product of theprevious two axis ratios and thus depends on both w and u . Whenboth w and u are negative, the density ρ ⋆ has its long-axis along the x -axis and its short-axis along the z -axis, in the same way as thepotential V S and the associated density ρ S . Above certain positivevalues of either w or u , the axis ratios become larger than unity,which means that ρ ⋆ is no longer aligned with the underlying co-ordinate system in the same way as V S and ρ S . For example, when ( − α ) w = − . and ( − α ) u = 0 . , b/a < but c/b > , so thatin this case ρ ⋆ has its short axis along the y -axis.A change in the sign of w and u has a strong effect on the ra-dial slope of ρ ⋆ . In Fig. 5, the radial profiles of ρ ⋆ along the princi-pal axes are shown for three combinations of w and u . The densityis normalised to the central value ρ . The profiles along the y -axis(dotted curves) and along the z -axis (dashed curves) are arbitrarilyoffset vertically with respect to the profile along the x -axis (solidcurves) to enhance visualisation. The thin curves are the profiles ofthe (luminous) mass density ρ ⋆ for varying δ , from δ = 0 (dark-est curve) to δ = 4 (lightest curve), in unit steps. The thick black c (cid:13) , 1–33 G. van de Ven et al.
Figure 4.
The characteristic axis ratios b/a , c/b and c/a of the luminous mass density for a non-rotating Abel component, as function of the DF parameters w and u , while δ = 1 . The axis ratios of the confocal ellipsoidal coordinate system are ζ = 0 . and ξ = 0 . , so that cf. (2.12) ( − α ) w ≥ − / ≈ − . and ( − α ) u ≤ / ≈ . . The thick contours are drawn at the levels that correspond to the characteristic axis ratios of the total mass density ρ S ,associated with the underlying isochrone St¨ackel potential (4.2), respectively b S /a S = 0 . , c S /b S = 0 . and c S /a S = 0 . . The intermediate-over-longaxis ratio b/a depends mainly on w , the short-over-intermediate axis ratio c/b depends mainly on u , and c/a is the product of the previous two. Figure 5.
Principal axes profiles of the luminous mass density ρ ⋆ for a non-rotating Abel component, normalised to the central value ρ ⋆, . Each panel isfor a different combination of the DF parameters w and u , while the grey scale indicates variation in δ from zero (darkest curve) to four (lightest curve), inunity steps. The profiles along the y -axis (dotted curves) and along the z -axis (dashed curves) are arbitrarily offset vertically with respect to the profile alongthe x -axis (solid curves) to enhance visualisation. The thick black curves show the profiles for the (total) mass density ρ S , associated with the underlyingisochrone St¨ackel potential (4.2), with ζ = 0 . and ξ = 0 . . When the value of either w or u is positive (right panel), the profiles show a break at r ∼ √− α ,so that these compact components may be used to represent kinematically decoupled components. curves show the profiles for the (total) mass density ρ S , which isindependent of w , u and δ .The profiles of ρ ⋆ steepen for increasing δ and for increasingabsolute values of w and u . In particular, when either w or u be-comes positive (right panel), the profiles suddenly become muchsteeper and drop to zero already at relatively small radii r ∼ √− α .The resulting Abel components are thus compact and, as we sawabove, can be different in shape and orientation from the main bodyof the galaxy model. Therefore, they can be used to represent kine-matically decoupled components. When both w ≤ and u ≤ (left and middle panel), ρ ⋆ falls off much more gently and the Abelcomponents cover a larger region. When w = u = 0 (left panel), so that the DF only depends on energy, the profiles as well as theshape (Fig. 4) of ρ ⋆ can even be flatter than those of ρ S . However,already for small non-zero values of w and u , generally ρ ⋆ ≤ ρ S everywhere in the galaxy model, and ρ ⋆ < ρ S in the outer parts.Although self-consistency ρ ⋆ = ρ S is only possible in the spher-ical case (for fixed values of w and u , see § w , u and δ so that ρ ⋆ ∼ ρ S . At the same time, hav-ing ρ ⋆ < ρ S in the outer parts of the galaxy model, allows for apossible dark halo contribution.The shape of ρ ⋆ can furthermore change due to the additionalcontribution from long-axis rotating and short-axis rotating compo-nents. Although these components have no density along their ro- c (cid:13) , 1–33 ecovery orbital structure tation axis, the behaviour of their overall shape as function of w , u and δ is similar as for the corresponding non-rotating components.The above analysis shows that, given the triaxial isochronepotential (4.2), we can use Abel components to construct a galaxymodel with a realistic density. Depending on the choice of w , u and δ , the galaxy model can contain compact (kinematically de-coupled) components and account for possible dark matter (in theouter parts). Furthermore, we show below that even with a smallnumber DF components, enough kinematic variation is possible tomimic the two-dimensional kinematic maps of early-type galaxiesprovided by observations with current integral-field spectrographs.This means that we can construct simple but realistic galaxy modelsto test our Schwarzschild software ( § As before, we choose the isochrone St¨ackel potential (4.2), we take ζ = 0 . and ξ = 0 . for the axis ratios of the coordinate system(4.4), resulting in a triaxiality parameter (2.6) of about T = 0 . ,and we set the scale length √− α = 10 ′′ . Assuming a distanceof D = 20 Mpc and a total mass of M ⊙ results in a centralvalue for the potential V ∼ . × km s − , which also setsthe unit of velocity. We restrict the number of DF components tothree, one of each type. For the first component of type NR we set w = u = − . / ( − α ) and δ = 1 , so that the shape of the corre-sponding density is similar to that of ρ S , except in the outer partswhere a steeper profile mimics the presence of dark matter (seeFigs. 4 and 5). For the second and third component, respectivelyof type LR and SR, we adopt the same parameters, expect that wetake w = 0 . / ( − α ) and u = − . / ( − α ) for the SR component,which therefore is more compact than the NR and LR component.We set the line-of-sight by choosing ϑ = 70 ◦ and ϕ = 30 ◦ for the viewing angles. After rotation over the misalignment an-gle ψ = 101 ◦ eq. (3.1), we compute for each DF component theLOSVD as a function of the positions on a rectangular grid on thesky plane, illustrated in Fig. 6 for five sky positions. By fitting aGauss-Hermite series to each LOSVD, we obtain the maps of themean line-of-sight surface mass density Σ , velocity V , dispersion σ and higher-order Gauss-Hermite moments h and h , shown inFig. 7. The parameters of each DF component are given on theright. The NR component has zero (green) odd velocity moments.For the LR and SR component, the even velocity moments showa decrease in the centre, because these components have zero den-sity along respectively the intrinsic long and short axis. We add theLOSVDs of the NR, LR and SR components, weighted with massfractions of respectively 80%, 12.5% and 7.5%, and fit a Gauss-Hermite series to obtain maps of Σ , V , σ , h and h . We convert Σ to the surface brightness by dividing by a constant stellar mass-to-light ratio of ( M ⋆ /L ) = 4 M ⊙ /L ⊙ .To convert these ‘perfect’ kinematics to ‘realistic’ observa-tions, similar to those obtained with integral-field spectrographssuch as SAURON (Bacon et al. 2001), we finally apply the followingsteps. We compute the kinematics on a rectangular grid consistingof 30 by 40 square pixels of 1 ′′ in size. Using the adaptive spatialtwo-dimensional binning scheme of Cappellari & Copin (2003), webin the pixels according to the criterion that each of the resulting(Voronoi) bins contains a minimum in signal-to-noise (S/N), whichwe take proportional to the square root of the surface brightness.For the mean errors in the kinematics we adopt the typical valuesof . km s − for V and σ and . for h and h in the kine-matics of a representative sample of early-type galaxies observedwith SAURON (Emsellem et al. 2004). We then weigh these val- ues with the S/N in each bin to mimic the observed variation inmeasurement errors across the field. Finally, we use the computedmeasurement errors to (Gaussian) randomise the kinematic maps.In this way, we include the randomness that is always present inreal observations. The resulting kinematic maps are shown in thetop panels of Fig. 8. Because of the eight-fold symmetry of thetriaxial model, the maps of the even (odd) velocity moments arealways point-(anti)-symmetric, apart from the noise added.
We briefly describe our numerical implementation of Schwarz-schild’s method in triaxial geometry (see vdB07 for a full descrip-tion), which we then use to fit the observables of the triaxial Abelmodel constructed in § The first step is to infer the gravitational potential from the observedsurface brightness. We do this by means of the Multi-Gaussian Ex-pansion method (MGE; e.g., Cappellari 2002), which allows forpossible position angle twists and ellipticity variations in the sur-face brightness. For a given set of viewing angles ( ϑ, ϕ, ψ ) (see § ( M/L ) to get the intrinsic mass density, fromwhich the gravitational potential then follows by solving Poisson’sequation. We calculate orbits numerically in the resulting gravita-tional potential.To obtain a representative library of orbits, the integrals of mo-tion have to be sampled well. The energy can be sampled directly,but since the other integrals of motion are generally not known,we start, at a given energy, orbits from a polar grid in the ( x, z ) -plane, which is crossed perpendicularly by all families of (regular)orbits. We restrict ourselves to the region in the first quadrant thatis enclosed by the equipotential and the thin orbit curves to avoidduplication of the tube orbits. To have enough box orbits to supportthe triaxial shape, we also start orbits by dropping them from theequipotential surface (Schwarzschild 1979, 1993).Assigning a mass weight γ j to each orbit j from the library,we compute their combined properties and find the weighted su-perposition that best fits the observed surface brightness and (two-dimensional) kinematics. However, the resulting orbital weight dis-tribution may vary rapidly, and hence probably corresponds to anunrealistic DF. To obtain a smoothly varying DF, we both ditherthe orbits by considering a bundle of integrated orbits that werestarted close to each other, and we regularise when looking for thebest-fit set of orbital weights by requiring them to vary smoothlybetween neighbouring orbits (in integral space). Finally, the best-fit Schwarzschild model follows from the minimum in the (Chi-squared) difference between (photometric and kinematic) observ-ables and the corresponding model predictions, weighted with theerrors in the observables. In this case, the gravitational potential is known and given by theisochrone St¨ackel potential V S eq. (4.1). However, to closely sim-ulate the Schwarzschild modelling of real galaxies, we infer thepotential from a deprojection of an MGE fit of the surface mass c (cid:13) , 1–33 G. van de Ven et al.
Figure 6.
Line-of-sight velocity distribution (LOSVD) at five different positions ( x ′ , y ′ ) on sky-plane (in arcsec at the top of each column) of three differentAbel components. The isochrone St¨ackel potential (4.2) is used, with ζ = 0 . and ξ = 0 . ( T = 0 . ), and scale length √− α = 10 ′′ . The model is placedat a distance of D = 20 Mpc and the adopted viewing angles are ϑ = 70 ◦ and ϕ = 30 ◦ . From top to bottom the LOSVDs of a non-rotating (NR), long-axisrotating (LR) and short-axis rotating (SR) Abel component are shown, with the corresponding DF parameters w , u and δ given on the right. The height ofeach LOSVD is normalised to unity, and a (dashed) Gaussian distribution with zero mean and the same dispersion as the LOSVD is shown as a reference. density Σ S generated by V S . The resulting potential reproduces V S to high precision, with relative differences less than − . Wecompute a library of orbits by sampling 21 energies E via a log-arithmic grid in radius from 1 ′′ to 123 ′′ that contains ≥ ( x, z ) -plane and drop box orbits from a similar uniform polargrid on the equipotential surface in the first octant. This results ina total of × × × starting positions, from each ofwhich a bundle of orbits are started. Taking into account the twosenses of rotation of the tube orbits, this results in a total orbits that are numerically integrated in the potential.We sum the velocities of each bundle of orbits in histogramswith 401 bins, at a velocity resolution of km s − . We fit theweighted sum of the velocity histograms to the intrinsic massdensity ρ ⋆ , which we obtain from a deprojection of an MGE fitto the observed surface brightness, multiplied with the (constant) ( M ⋆ /L ) = 4 M ⊙ /L ⊙ . Simultaneously, we fit the projected val-ues of the velocity histograms to the observed surface brightnessand higher-order velocity moments. Finally, at the same time, weregularise the orbital weights in E and in the starting positions byminimising their second order derivatives. The strenght of the reg-ularisation is given by the a smoothening parameter (e.g., Crettonet al. 1999), which we set to λ = 0 . (see vdB07).From Fig. 8 it is clear that the (simulated) observables of thetriaxial Abel model (top panels) are very well matched by the best-fit triaxial Schwarzschild model (bottom panels). The signature of the kinematically decoupled component in the maps of the meanline-of-sight velocity V and Gauss-Hermite moment h is accu-rately fitted, as well as the kinematics of the main body up to h within the added noise ( § § We calculate the intrinsic first and second order velocity momentsof the Schwarzschild model by combining the appropriate momentsof the orbits that receive weight in the superposition, and investigatehow well they compare with the intrinsic velocity moments of theAbel model. In general, there are three first h v t i and six secondorder velocity moments h v s v t i ( s, t = x, y, z ). Combining themyields the six dispersion components σ st of the velocity dispersiontensor, where σ st ≡ h v s v t i − h v s ih v t i .We first consider the ( x, z ) -plane, as it is crossed perpendicu-larly by all four (major) orbit families. Because h v x i = h v z i = σ xy = σ yz = 0 , we are left with h v y i perpendicular to the ( x, z ) -plane as the only non-vanishing mean motion and σ zx inthe ( x, z ) -plane as the only non-vanishing cross-term. The averageroot-mean-square velocity dispersion σ RMS is given by σ = c (cid:13) , 1–33 ecovery orbital structure Figure 7.
Maps of the surface brightness (SB; in L ⊙ pc − ), mean line-of-sight velocity V and dispersion σ (both in km s − ), and higher order Gauss-Hermite moments h and h , of the same three Abel DF components as in Fig. 6, obtained by fitting a Gauss-Hermite series to the LOSVDs at each (pixel)position on the plane of the sky. The numerical artifacts at the edges of the compact SR component (third row) disappear when combined with componentsthat extent over the full field-of-view (see e.g. the top row of Fig. 8). ( σ xx + σ yy + σ zz ) / . The ratio h v y i /σ RMS of ordered-over-random motion is a measure of the importance of rotation for thegravitational support of a galaxy. In Fig. 9, the colours represent thevalues of this ratio in the ( x, z ) -plane, for the input triaxial Abelmodel (left panel) and for the best-fit triaxial Schwarzschild model(right panel).In a St¨ackel potential the axes of the velocity ellipsoid arealigned with the confocal ellipsoidal coordinate system (e.g., Ed-dington 1915; van de Ven et al. 2003). As a result, one of the axesof the velocity ellipsoid is perpendicular to the ( x, z ) -plane, withsemi-axis length σ yy . The other two axes lie in the ( x, z ) -plane andhave semi-axis lengths given by σ ± = ( σ xx + σ zz ) ± q ( σ xx − σ zz ) + σ xz . (5.1)The ellipses overplotted in Fig. 9 show the corresponding cross sec-tions of the velocity ellipsoid with the ( x, z ) -plane. The flatteningof the ellipses is thus given by the ratio σ − /σ + , while the angle θ xz of the major-axis with respect to the x -axis is given by tan(2 θ xz ) = 2 σ xz / ( σ xx − σ zz ) . (5.2)In addition, the cross on top of each ellipse represents the ratio In case of alignment with the confocal ellipsoidal coordinate system,this angle is given by the tangent to the curves of constant ( µ, ν ) , i.e., tan θ xz = ( z/x )( λ + α ) / ( λ + γ ) , which indicates approaching align-ment with the polar coordinate system at large radius. σ yy /σ + , i.e., the (relative) size of the velocity ellipsoid in the per-pendicular direction. For an isotropic velocity distribution the el-lipses become circles and the crosses fill the circles. Finally, theblack curves are contours of constant mass density in steps of onemagnitude.The density of the triaxial Abel model (solid curve) is wellfitted by the triaxial Schwarzschild model (dashed curve), with a(biweight ) mean fractional difference below %. In both the Abelmodel and the fitted Schwarzschild model the value of h v y i /σ RMS is relatively low, with a mean value ∼ . , indicating that gravi-tational support is mainly due to random motion. Still, the averagerotation of the long-axis tube orbits (with h v y i < ) due to themaximum streaming LR component in the input Abel model, aswell as, the opposite maximum streaming of the (compact) short-axis rotating component are clearly visible, and well recovered bythe best-fit Schwarzschild model. The average absolute differencein both h v y i and σ RMS is below km s − , and thus well within thetypical error of . km s − assigned to the simulated mean line-of-sight velocity V and velocity dispersion σ of the Abel model (see § h v y i /σ RMS is ∼ . .We see in Fig. 9 that, at larger radii, the ellipses become moreradially elongated and the relative size of the crosses decreases inthe radial direction, but they stil fill the ellipses in the angular direc- The biweight mean (e.g., Andrews et al. 1972; Beers, Flynn & Gebhardt1990) is robust estimators for a broad range of non-Gaussian underlyingpopulations and is less sensitive to outliers than other moment estimators.c (cid:13) , 1–33 G. van de Ven et al.
Figure 8.
Kinematic maps for a triaxial Abel model (top row) and converted to observables with realistic measurement errors added (middle row; see § § V and dispersion σ (both in km s − ),and Gauss-Hermite moments h and h . Isophotes of the surface brightness of the Abel model are overplotted in each map. At the right side of each map, the(linear) scale of the corresponding kinematics is indicated by the colour bar, and the limits are given below. tion. This implies a velocity distribution that becomes increasinglyradially anisotropic outwards, but remains close to isotropic in thetangenetial direction everywhere. This shape and orientation of thevelocity ellipsoid in the input Abel model is well reproduced by thebest-fit Schwarzschild model, with only a (mild) underestimationof the radial anisotropy towards the z -axis. This is likely the re-sult of numerical difficulties due to the small number of (sampled)long-axis tube orbits that contribute in this region. The absolutedifference in the semi-axis lengths σ + , σ − and σ yy of the velocityellipsoid is on average ∼ km s − . This uncertainty includes bothdeviations in shape and orientation of the velocity ellipsoid, and iswihtin the expected range due to the errors in the simulated kine-matics. The corresponding axis ratios σ − /σ + and σ yy /σ + of thevelocity ellipsoid are on average recovered within ∼ %. Away from the ( x, z ) -plane, the average fractional differ-ence in the density between the input Abel model and the best-fitSchwarzschild model stays below %. Fig. 10 compares the in-trinsic moments of the input triaxial Abel model (top) with thoseof the best-fit triaxial Schwarzschild model (bottom) in three di-mensions. The first column shows the (amplitude) of the streamingmotion v str , given by v = v x + v y + v z , and normalised by σ RMS . These quantities are computed on a polar grid ( r, θ, φ ) inthe first octant. The (logarithmic) sampling of the radius r is indi-cated by the black dots between the top and bottom panels, whileeach row is for a different polar angle θ as indicated on the right,and the colours represent the (linear) change in azimuthal angle φ .The limit φ = 0 ◦ (black curves) corresponds to the ( x, z ) -planediscussed above. The resulting ordered-over-random motion V /σ c (cid:13) , 1–33 ecovery orbital structure Figure 9.
The colours represent the mean motion h v y i perpendicular to the ( x, z ) -plane, normalised by σ RMS (excluding the axes to avoid numericalproblems), for the input triaxial Abel model (left) and for the best-fit triaxial Schwarzschild model (right). The ellipses are cross sections of the velocityellipsoid with the ( x, z ) -plane and the crosses represent the (relative) size of the velocity ellipsoid in the perpendicular ( y -axis) direction. The black curvesare contours of constant mass density in steps of one magnitude, for the input Abel model (solid) and for the fitted Schwarzschild model (dashed). See § is well recovered by the Schwarzschild model, apart from the upperpanel, which is likely the result of the above mentioned numericaldifficulties close to the z -axis. Overall, the average absolute differ-ence in both v str and σ RMS is below km s − and the uncertaintyin v str /σ RMS is ∼ . .The second and third column of Fig. 10 show respectively theintermediate-over-major σ b /σ a and minor-over-major σ c /σ a axisratios of the velocity ellipsoid. The velocity ellipsoid of the triax-ial Abel model is aligned with the confocal ellipsoidal coordinatesystem, so that its semi-axis lenghts σ a ≥ σ b ≥ σ c follow directlyfrom σ τ = h v τ i − h v τ i with τ = λ, µ, ν . In general, this is notthe case for the triaxial Schwarzschild model, and instead we diago-nalize the (symmetric) velocity dispersion tensor with components h σ st i ( s, t = x, y, z ). As before, the axis ratios of the velocity ellip-soid are quite well recovered by the best-fit Schwarzschild model,except towards the z -axis (upper panels) where it underestimatesthe anisotropy in the velocity distribution of the input Abel model.Similarly, away from the ( x, z ) -plane ( φ = 0 ◦ , black curves), theSchwarzschild model increasingly overestimates the σ b /σ a ratio,while the σ c /σ a remains well reproduced. It is plausible that therecovery in the ( x, z ) -plane is better, because it is optimally sam-pled as starting space for the numerical orbit calculations, and itis crossed perpendicularly by all four major orbit families. Never-theless, the absolute difference in σ a , σ b and σ c between the in-put Abel model and the best-fit Schwarzschild model is on average ∼ km s − . The axis ratios σ b /σ a and σ c /σ a are on average re-covered within ∼ %. The fitted triaxial Schwarzschild model results in a mass weight γ per orbit. These mass weights are a function of the three integralsof motion ( E, I , I ) . In general, only the energy is exact, but fora separable potential I and I are also known explicitly and givenby eq. (2.5). The orbital mass weights follow from the DF by in-tegrating f ( E, I , I ) over the part of phase-space ( x , v ) that isaccesible by the orbit. Since each orbit is a (unique) delta-functionin integral-space, the resulting orbital mass weights are in principlezero. However, as described in § § γ ( E, I , I ) = ZZZ cell f ( E, I , I ) ∆ V ( E, I , I ) d E d I d I , (5.3)where ∆ V ( E, I , I ) = ZZZ Ω ˛˛˛˛ ∂ ( v x , v y , v z ) ∂ ( E, I , I ) ˛˛˛˛ d x d y d z, (5.4)with Ω the volume in configuration space accessible by the orbit.The multi-component DF of the input triaxial Abel model consistsof basis functions defined in eq. (2.24), with the DF parameters andweights per component given in § ∆ V and the cell in integral space, and then return to the comparison ofthe orbital mass weights. c (cid:13) , 1–33 G. van de Ven et al.
Figure 10.
Intrinsic velocity moments in three dimensions for the input triaxial Abel model (top) and for the best-fit triaxial Schwarzschild model (bottom).The first column shows the (amplitude) of the streaming motion v str , normalised by σ RMS . The second and third column show the axis ratios of the velocityellipsoid, where σ a , σ b and σ c are respectively the semi-lengths of the major, intermediate and minor principal axes. These quantities are computed on apolar grid ( r, θ, φ ) in the first octant. The (logarithmic) sampling of the radius r (in arcsec) is indicated by the black dots between the top and bottom panels.Each row is for a different polar angle θ (in degrees) as indicated on the right, with the top panel close to the z -axis and the bottom panel close to the ( x, y ) -plane. The colours represent the (linear) change in azimuthal angle φ (in degrees), with limits ◦ and ◦ corresponding to the ( x, z ) -plane and ( y, z ) -plane,respectively. See § (cid:13) , 1–33 ecovery orbital structure Figure 11.
Three quantities involved in the calculation of the orbital mass weights for a triaxial Abel model with an isochrone potential. For a given energy E , in each panel, the values of the second and third integral of motion, I and I , indicated by the symbols, correspond to the orbital starting position andvelocities in the triaxial Schwarzschild model that is fitted to the observables of this triaxial Abel model. The solid curves, calculated with the expressionsin § ( x, z ) -plane and the triangles represent the additional set of orbits dropped from the equipotential surface (see § ( x, z ) -plane. The colours in the left panel indicate the value of the DF f ( E, I , I ) for each orbit in unitsof M ⊙ (km/s) − pc − , with the (linear) scale given by the vertical bar on the right. The colours in the middle panel represent the values of ∆ V ( E, I , I ) defined in eq. (5.4) in units of (km/s) − pc − . The area of each Voronoi bin in the right panel , multiplied by the range in energy E , approximates the cell ∆ E ∆( I , I ) in integral space for each orbit. The product of these three quantities yields an estimate of the mass weight γ ( E, I , I ) for each orbit. The expression for ∆ V ( E, I , I ) of a single orbit in a triaxialSt¨ackel potential can be deduced from the relations in § ∆ V ( E, I , I ) = ( γ − α ) ZZZ Ω ( λ − µ )( µ − ν )( ν − λ ) a ( λ ) a ( µ ) a ( ν ) × s λ + β )( µ + β )( ν + β )[ E − V eff ( λ )] [ E − V eff ( µ )] [ E − V eff ( ν )] d λ d µ d ν, (5.5)where a ( τ ) , τ = λ, µ, ν , is defined as a ( τ ) = ( τ + α )( τ + β )( τ + γ ) , (5.6)the effective potential V eff as V eff ( τ ) = I τ + α + I τ + γ + U [ τ, − α, − γ ] , (5.7)and Ω is the volume in configuration space accessible by the orbitin the triaxial separable potential that obeys ( E, I , I ) . The lastterm in eq. (5.7) is equal to the St¨ackel potential (2.3) along theintermediate y -axis.Because of the separability of the equations of motion, eachorbit in a triaxial separable potential can be considered as a sumof three independent motions. Each of these one-dimensional mo-tions is either an oscillation or rotation in one of the three confocalellipsoidal coordinates ( λ, µ, ν ) , such that the configuration spacevolume Ω is bounded by the corresponding coordinate surfaces.The values of ( λ, µ, ν ) that correspond to these bounding surfacescan be found from Table 1 for the four families of regular orbits:boxes (B), inner (I) and outer (O) long-axis tubes, and short-axis (S)tubes. Whereas α , β and γ are the limits on ( λ, µ, ν ) set by the fociof the confocal ellipsoidal coordinate system, the other limits arethe solutions of E = V eff ( τ ) (see Fig. 7 of de Zeeuw 1985a). In the case of the triaxial isochrone St¨ackel potential (4.2), we can writethis equation as a fourth-order polynomial in √ τ . The solutions arethen the squares of three of the four roots of this polynomial (thefourth root is always negative).For each orbit in our Schwarzschild model, we compute ( E, I , I ) by substituting the starting position and velocities of theorbit into the expressions (2.5). From the value of E and the signof I (while always I ≥ ), we determine to which orbit family itbelongs. The corresponding configuration space volume Ω is thengiven by the boundaries for λ , µ and ν in the last three columns ofTable 1. The value of ∆ V ( E, I , I ) follows by numerical evalua-tion of the right-hand side of eq. (5.5).The integrand in eq. (5.5) contains singularities at the integra-tion limits, which can be easily removed for a triaxial isochronepotential. We write the integrand completely in terms of ( √ σ ±√ τ ) / , where σ, τ = λ, µ, ν or a constant value. Suppose nowthat the integral over λ ranges from λ to λ and the terms ( √ λ −√ λ ) / and ( √ λ − √ λ ) / appear in the denominator. The sub-stitution √ λ = √ λ + ( √ λ − √ λ ) sin η then removes bothsingularities since dλ/ [( √ λ − √ λ )( √ λ − √ λ )] / = 4 √ λ dη . We approximate the triple integration over the cell in integral spacein eq. (5.3) by the volume ∆ E ∆( I , I ) . Here ∆ E is the (log-arithmic) range in E between subsequent sets of orbits at differ-ent energies (see § E = 0 . Because we do not directly sample I and I in our implementation of Schwarzschild’s method, as their expres-sions are in general unknown, we cannot directly calculate the area ∆( I , I ) . Instead, we compute the Voronoi diagram of the pointsin the ( I , I )-plane that correspond to the starting position and ve-locities of each orbit, at a given energy E . An example is given in c (cid:13) , 1–33 G. van de Ven et al.
Table 1.
Configuration space volume Ω accessible by the four families of regular orbits.orbit I E λ µ ν
Box orbits < V eff ( − β ) . . . − α . . . λ max − β . . . µ max − γ . . . ν max Inner long-axis tube orbits < V eff ( µ )] . . . V eff ( − β ) − α . . . λ max µ min . . . µ max − γ . . . − β Outer long-axis tube orbits > V eff ( λ )] . . . V eff ( − β ) λ min . . . λ max µ min . . . − α − γ . . . − β Short-axis tube orbits > { V eff ( − β ) , min [ V eff ( λ )] } . . . λ min . . . λ max − β . . . − α − γ . . . ν max the right panel of Fig. 11. The area of the Voronoi bins approxi-mates the area ∆( I , I ) for each orbit.The four families of regular orbits are separated by two linesthat follow from I = 0 and E = V eff ( − β ) . The latter provides thepart of the boundary on I and I for the box orbits. The remainderof this boundary is given by the positivity constraint on I and bythe solution of (cf. eqs 64 and 65 of de Zeeuw 1985a) E = V eff ( κ ) and » d V eff ( κ )d κ – κ = 0 , κ ≥ − β. (5.8)Substituting V eff from eq. (5.7) and using d U [ τ, − α, − γ ] / d τ = U [ τ, τ, − α, − γ ] , we find the solution I = ( κ + α ) ( α − γ ) n E − U [ − α, κ , κ ] o , (5.9)and similarly for I by interchanging α ↔ γ . For − β ≤ κ ≤ − α ,the solution describes the boundary curve for which I ≤ andcorresponds to the thin I tube orbits. For κ ≥ − α , we find theboundary curve for which I ≥ , corresponding to the thin O andS tube orbits.There are limits on the values of κ depending on the valueof E , and sometimes there are no valid solutions for κ , whichimplies that only box orbits contribute at that energy. These limitscan be obtained from the thin orbit curves in the ( x, z ) -plane. With y = v x = v z = 0 , the expressions (2.5) for the integrals of motionreduce in this plane to E = v y + U [ λ, κ, − β ] ,I = x n v y + ( α − β ) U [ λ, κ, − β, − α ] o , (5.10) I = z n v y + ( γ − β ) U [ λ, κ, − β, − γ ] o , with − γ ≤ κ ≤ − α replacing µ and ν respectively above andbelow the focal curve given by z / ( γ − β ) − x / ( β − α ) = 1 .Next, we substitute the expression for E in those for I and I andwe use that ( τ + β ) U [ λ, κ, τ, − β ] = U [ λ, κ, τ ] − U [ λ, κ, − β ] ,respectively for τ = − α and τ = − γ . We find that the thin orbitcurves follow by solving I = 0 and thus E = U [ λ, κ , κ ] forI tubes, and I = 0 and thus E = U [ κ , κ , κ ] , with κ = µ forO tubes and κ = ν for S tubes. In general these equations haveto be solved numerically, but in the case of the triaxial isochronepotential (4.2), they reduce to a second order polynomial in √ κ and the solutions simply follow from the roots of the polynomial. Once we have computed for each orbit the DF f ( E, I , I ) , ∆ V ( E, I , I ) and the cell ∆ E ∆( I , I ) in integral space(Fig. 11), its (approximate) mass weight γ ( E, I , I ) follows bymultiplication of these three quantities. As before, the choice ofmaximum streaming for the (LR and SR) rotating components re- duces the accessible integral space, and thus also the correspondingorbital mass weights, by a factor two.The resulting orbital mass weight distribution of the input tri-axial Abel model is shown in the top panels of Fig. 12, and thatof the fitted triaxial Schwarzschild model in the bottom panels. Theenergy E increases from left to right, which corresponds to increas-ing distance from the centre as is indicated by the radius R E (inarcsec) at the top of each panel. For this representative radius weuse the radius of the corresponding thin (S) tube orbit on the x -axis.The values of I and I , on the horizontal and vertical axes respec-tively, are both normalised with respect to their maximum ampli-tude at the given energy. In each panel, the mass weight values arenormalised with respect to the maximum in that panel. Between thetwo rows of panels, the fraction of the summed values in each panelwith respect to the total mass weight in all panels is given (in %).The panels with R E < ∼ ′′ are best constrained by the kine-matic observables. This takes into account that even orbits that ex-tend beyond the maximum radius covered by the observables cancontribute significantly at smaller radii. In these panels, the mainfeatures of the orbital mass weight distribution of the triaxial Abelmodel are recovered. In the outer parts the Schwarzschild modelis still constrained by the mass model, which extends to a radiusof about 100 ′′ , but the orbital mass weight distribution deviatesfrom that of the input Abel model due to the lack of kinematic con-straints. A point-by-point comparison yields an average fractionalerror of ∼ %, and if we consider in each panel the mass weightsabove the mean value, which together contribute more than halfof the total mass, the fractional error decreases to ∼ %. How-ever, this way of quantifying the recovery is (somewhat) mislead-ing since the relatively large fractional errors are at least partiallycaused by the strong peaks in the orbital mass weight distribution.For example, if in the input Abel model a certain orbit gets a sig-nificant weight, but in the Schwarzschild model, due to numericaluncertainties, this weight is assigned to a neighboring orbit witha (slightly) different value of I , the relative error at each of thecorresponding points in the integral space can be very large.Henceforth, we show in Fig. 13 the orbital mass weights asfunction of each of the three integrals of motion separately by col-lapsing the cube in ( E, I , I ) in the remaining two dimensions.We again use R E as a representative radius for E (first panel), butsince the (range of) values for I and I change with E (see alsoFig. 12), we use their index in the cube instead. In addition to thetotal distribution, we also show the contribution of the three differ-ent NR, LR and SR components separately, as well as for the lattertwo rotating components the contributions from the two directionsof rotation by making the mass weights for one of the directionsnegative. Since the input triaxial Abel model (diamonds connectedby solid curves) is constructed with maximum streaming in one ofthe two directions for both the LR and SR component (see § c (cid:13) , 1–33 ecovery orbital structure Figure 12.
The orbital mass weight distribution for the input triaxial Abel model (top) and for the best-fit triaxial Schwarzschild model (bottom). From leftto right the energy increases, corresponding to increasing distance from the centre, indicated by the radius R E (in arcsec) of the thin short-axis tube orbit onthe x -axis. The vertical and horizontal axes represent respectively the second and third integral of motion, I and I , normalised by their maximum amplitude(for given E ). In each panel, the colours represent the mass weights, normalised with respect to the maximum in that panel, and with the (linear) scale givenby the vertical bars on the right. Between the two rows of panels, the fraction (in %) of the included mass with respect to the total mass is indicated. (crosses connected by dotted curves) in which ∼ % of the totalmass, or ∼ % of the mass of the LR and SR components, iswrongly assigned to the opposite direction. Keeping in mind thatthe orbital mass weights itself are not directly fitted and that thetypical velocity error of 7.5 km s − is more than % of the max-imum in the simulated velocity field (see § E is well recovered, even in the outer parts where (nearly) allthe constraints come only from the mass model. The average abso-lute difference is ∼ . %. Whereas for E the constraints providedby the mass model already seem sufficient, for I and I the kine-matic constraints are essential. Not suprisingly, we then also seethat the recovery is less good with an average absolute differenceof ∼ . % in I and ∼ . % in I . The main contribution is fromthe NR component, while the two rotating components seem to bet-ter recovered. We now consider three-integral galaxy models in the axisymmet-ric limit. As we have seen in the Introduction (Section 1), variousgroups have successfully developed independent axisymmetric im-plementations of Schwarzschild’s method and verified their codesin a number of ways. The published tests to recover a known (ana-lytical) input model have been limited to spherical geometry or toan axisymmetric DF that is a function of the two integrals of motion E and L z only.Here, we present the velocity moments of the three-integral Abel DF in the axisymmetric limit and we choose again theisochrone form in eq. (4.1) for the St¨ackel potential. The proper-ties of the resulting three-integral Kuzmin-Kutuzov models can beexpressed explicitly in cylindrical coordinates. We construct an ax-isymmetric oblate Abel model and fit Schwarzschild models to theresulting observables to test how well the axisymmetric implemen-tation of Schwarzschild’s method, as presented in Cappellari et al.(2006), recovers the intrinsic velocity moments as well as the three-integral DF. When two of the three constants α , β or γ are equal, the confocalellipsoidal coordinates ( λ, µ, ν ) reduce to spheroidal coordinatesand the triaxial St¨ackel potential (2.3) becomes axisymmetric. When β = α = γ (triaxiality parameter T = 0 ), we cannot use µ as a coordinate and replace it by the azimuthal angle φ , defined as tan φ = y/x . The relation between ( λ, φ, ν ) and the usual cylin-drical coordinates ( R, φ, z ) is given by R = ( λ + α )( ν + α ) α − γ , z = ( λ + γ )( ν + γ ) γ − α . (6.1)The St¨ackel potential V S ( λ, ν ) = U [ λ, − α, ν ] is oblate axisym-metric . The corresponding integrals of motion follow by substitu-tion of µ = − β = − α in the expressions (2.5), so that the secondintegral of motion reduces to I = L z .With the choice (2.9) for the DF, the expression for ve-locity moments µ lmn ( λ, ν ) is that of the triaxial case given in c (cid:13) , 1–33 G. van de Ven et al.
Figure 13.
The orbital mass weights (in % of the total mass) for the input triaxial Abel model (diamonds connected by solid curves) and for the best-fit triaxialSchwarzschild model (crosses connected by dotted curves), as function of each of the three integrals of motion. These ’projections’ of the three-dimensionalorbital mass weight distribution shown in Fig. 12 are obtained by collapsing the cube in ( E, I , I ) in two dimensions. As before, we represent the energy E in the first panel by the radius R E (in arcsec) of the thin short-axis tube orbit on the x -axis. For the second and third integral of motion, I and I , we usethe index in the cube, since the (range of) their values changes with E . The total distribution (black colour) is split into contributions from the non-rotating(NR; red), long-axis rotating (LR; green) and short-axis rotating (SR; blue) components. Moreover, for each rotating component the contributions from thetwo directions of rotation are separated by making the mass weights for one of the directions negative. eq. (2.10), but with µ = − β = − α . From Fig. 1, we see thatthe lower limit on w vanishes. For the NR type of components, S max = S top ( λ, µ, − γ ) and the corresponding velocity moments µ NR lmb ( λ, ν ) vanish when either l , m or n is odd. Because the onlyfamily of orbits that exists are the short-axis tube orbits, we canintroduce net rotation (around the short z -axis) by setting the DFto zero for L z < , so that µ SR lmn ( λ, ν ) = µ NR lmn ( λ, ν ) . These SRvelocity moments vanish when either l or n is odd, but are non-zeroif m is odd. They should be multiplied with ( − m for maximumstreaming in the opposite direction. By choosing different weightsfor both senses of rotation, we can control the direction and theamount of streaming motion.In the conversion to observables described in §
3, the ma-trix Q , which transforms the velocity components ( v λ , v φ , v ν ) to ( v x , v y , v z ) , reduces to Q = A cos φ − sin φ − B cos φA sin φ cos φ − B sin φB A , (6.2)where A and B are defined as A = ( λ + γ )( ν + α )( λ − ν )( α − γ ) , B = ( λ + α )( ν + γ )( λ − ν )( γ − α ) . (6.3)Because of the symmetry around the short-axis, the azimuthalviewing angle ϕ looses its meaning and the misalignment angle ψ = 0 ◦ . We are left with only the polar viewing angle ϑ , which iscommonly referred to as the inclination i , with i = 0 ◦ face-on and i = 90 ◦ edge-on viewing. As a consequence, the projection matrix P is a function of i only and follows by substituting ϑ = i and ϕ = 0 in eq. (3.4). The rotation matrix R in eq. (3.5) reduces tothe identity matrix, so that M = PSQ .The expression for the LOSVD follows from that of the tri-axial case in eq. (3.27) by substituting µ = − β = − α . For theNR components, again ∆ ξ ′ NR = 2 π and the simplified expression(3.28) holds in case of a DF basis function as defined in eq. (2.24).To introduce net rotation, we require that ( v µ =) v φ ≥ as in § ∆ ξ ′ SR is the length of the partof the circle between the intersections ξ ± = 2 arctan( u ± ) with theline (with u ± given in eq. 3.38), and which is on the correct side ofthe line in eq. (3.37). This is again similar to SR components in thetriaxial case, but without the restriction to stay within the ellipses. When β = γ = α ( T = 1 ), we replace the coordinate ν bythe angle χ , defined as tan χ = z/y . The resulting coordinates ( λ, µ, χ ) follow from the above coordinates ( λ, φ, ν ) by taking ν → µ , φ → χ , and γ → α → β . The St¨ackel potential V S ( λ, µ ) = U [ λ, µ, − γ ] is now prolate axisymmetric . By substi-tuting ν = − β = − γ in eqs (2.5) and (2.10), we obtain the ex-pressions respectively for the integrals of motion (with I = L x )and for the intrinsic velocity moments µ lmn ( λ, µ ) . From Fig. 1, wesee that now the upper limit on u vanishes. For the NR components, S max = S top ( λ, µ, − γ ) , and since we only have the long-axis tubeorbits, we can introduce net rotation (around the x -axis) by settingthe DF to zero for L x < , so that µ LR lmn ( λ, µ ) = µ NR lmn ( λ, µ ) .These LR velocity moments vanish if either l or m is odd and mul-tiplication with ( − n yields net rotation in the opposite direction.The matrix Q , which transforms ( v λ , v µ , v χ ) to ( v x , v y , v z ) ,in this case reduces to Q = C − D D cos χ C cos χ − sin χD sin χ C sin χ cos χ , (6.4)where C and D are given by C = ( λ + β )( µ + α )( λ − µ )( α − β ) , D = ( λ + α )( µ + β )( λ − µ )( β − α ) . (6.5)In the projection matrix P in eq. (3.4), we substitute ϑ = π/ − i and ϕ = 0 , so that for inclination i = 0 ◦ and i = 90 ◦ , we are c (cid:13) , 1–33 ecovery orbital structure respectively viewing the prolate mass model end-on and side-on.In the rotation matrix R we take ψ = 90 ◦ to align the projectedmajor axis horizontally. The expression for the LOSVD followsfrom eq. (3.27) by substituting ν = − β = − γ , and by requir-ing ( v ν =) v χ ≥ we obtain LR components with maximumstreaming. As for the oblate case and illustrated in the left panelof Fig. 3, ∆ ξ ′ SR is the length of the circle part between the angles ξ ± = 2 arctan( u ± ) (with u ± given in eq. 3.34) which is on thecorrect side of the line in eq. (3.32). In the axisymmetric limit, the form (4.1) for U ( τ ) results in theKuzmin-Kutuzov (1962) potential. We give the properties relevantfor our analysis, while further details can be found in Dejonghe &de Zeeuw (1988), including expressions and plots of the mass den-sity ρ S , its axis ratios, and the two-integral DF f ( E, L z ) consistentwith ρ S [see also Batsleer & Dejonghe (1993), who also correcteda typographical error in f ( E, L z ) ].When β = α , the oblate axisymmetric potential V S ( λ, ν ) = U [ λ, − α, ν ] and the third order divided difference U [ λ, − α, ν, σ ] ,which both appear in the expressions for the integral of motions(2.5), have the simple forms V S ( λ, ν ) = − GM √ λ + √ ν , (6.6) U [ λ, − α, ν, σ ] = GM ( √ λ + √ ν )( √ λ + √ σ )( √ ν + √ σ ) , (6.7)where again GM = √− γ + √− α , so that V S = − in the centre.By means of the relations λ + ν = R + z − α − γ, λν = αγ − γR − αz , (6.8)and ( √ λ + √ ν ) = λ + ν + 2 √ λν and ( √ λ + √ σ )( √ ν + √ σ ) = √ λν + √ σ ( √ λ + √ ν ) + σ , we can write the potential and in-tegrals of motion explicitly as elementary expressions in the usualcylindrical coordinates.When β = γ , the prolate potential V S ( λ, µ ) = U [ λ, µ, − γ ] and the third order divided difference U [ λ, µ, − γ, σ ] follow respec-tively from (6.6) and (6.7) by replacing ν by µ . The above constructed triaxial Abel model ( § ζ approach unity,while keeping ξ = 0 . fixed. Similar to the triaxial case, the DFcontains a NR component with the same parameters, u = w = − . / ( − α ) and δ = 1 , but we exclude the LR component sincelong-axis tube orbits do not exist in an oblate axisymmetric galaxy.We include two SR components, one with the same parameters asthe NR component, and for the other we set w = 0 . / ( − α ) and u = − . / ( − α ) , and we choose the sense of rotation in the oppo-site direction. The latter implies a compact counter-rotating com-ponent, which is clearly visible in the kinematic maps shown in thetop panels of Fig. 14. The inclination is the same value as the po-lar angle ϑ for the triaxial Abel model, i.e. i = 70 ◦ , and the massfractions of the three DF components are respectively 20%, 60%and 20%. Due to axisymmetry, all maps of the even (odd) veloc-ity moments are bi-(anti)-symmetric and the velocity field shows astraight zero-velocity curve. The signatures of the counter-rotationare similar in the velocity field and h (but anti-correlated), andresult in a decrease of σ and an increase of h in the centre. We now describe the application of our axisymmetric implementa-tion of Schwarzschild’s method to the observables of the oblate ax-isymmetric Abel model of § We use the implementation of Schwarzschild’s method in axisym-metric geometry that is described in detail in Cappellari et al.(2006). The main differences with respect to our triaxial implemen-tation are certain simplifications due to the extra symmetry. Thereare no twists in the surface brightness and of the four families ofregular orbits only the short-axis tube orbits are supported. We usethe same set-up as in the triaxial case, but since there are no boxorbits, the additional dropping of orbits from the equipotential sur-face is not needed.Fig. 14 shows that the (simulated) observables of the oblateaxisymmetric Abel model (top panels) are very well matched bythe best-fit axisymmetric Schwarzschild model (bottom panels).The kinematics of the main body as well as the signatures of thecounter-rotating core are accurately fitted within the (added) noise.
It is convenient to analyse the intrinsic velocity moments of (oblate)axisymmetric models in cylindrical coordinates ( R, φ, z ) . Becauseof axisymmetry the models are independent of the azimuthal angle φ , and it is sufficient to consider the meridional ( R, z ) -plane. Theanalysis of the intrinsic velocity moments in the ( R, z ) -plane issimilar to that for the triaxial case in the ( x, z ) -plane ( § h v φ i , perpendicular to the merid-ional plane, is the only non-vanishing first-order velocity moment.In Fig. 15, we compare the values of h v φ i /σ RMS , indicated by thecolours, for the Abel model (left panel) with those for the fittedSchwarzschild model (right panel). The root-mean-square velocitydispersion σ RMS is defined as σ = ( σ R + σ φ + σ z ) / . Theazimuthal axis of the velocity ellipsoid, with semi-axis length σ φ defined as σ φ = h v φ i − h v φ i , is perpendicular to the meridionalplane. The cross sections with the meridional plane are indicatedby the ellipses in Fig. 15, where the semi-axis lengths follow from(5.1) by replacing ( x, z ) with ( R, z ) .As in the triaxial case the density (solid curve) is well fitted bythe axisymmetric Schwarzschild model (dashed curve). The Abelmodel shows a strong gradient in h v φ i /σ RMS , which is correctlyrecovered by the axisymmetric Schwarzschild model. The absolutedifference is on average less than . , except near the symmetry z -axis. This is likely the result of numerical difficulties due to thesmall number of (sampled) short-axis tube orbits that contributein this region. The shape and orientation of the ellipses are nearlyidentical, indicating that the anisotropic velocity distribution of theAbel model is reproduced within the expected uncertainties due tothe errors in the simulated kinematics. The axis ratios σ − /σ + and σ φ /σ + of the velocity ellipsoid are on average recovered within ∼ %. In the oblate axisymmetric case, all (regular) orbits are short-axis tube orbits with I = L z and energy E ranging from c (cid:13) , 1–33 G. van de Ven et al.
Figure 14.
Kinematic maps for an oblate axisymmetric Abel model (top row) and converted to observables with realistic measurement errors added (middlerow; see § § min [ V eff ( λ )] to zero. The expression for ∆ V in eq. (5.5) reducesto ∆ V ( E, L z , I ) = 4 π | L z | ν max Z − γ λ max Z λ min ( ν − λ )( λ + α )( λ + γ )( ν + α )( ν + γ ) × s ( λ + α )( ν + α )[ E − V eff ( λ )] [ E − V eff ( ν )] d λ d ν (6.9)where as before ν max , λ min and λ max are the solutions of E = V eff ( τ ) (see Fig. 23 of de Zeeuw 1985a). The factor in front of thedouble integral includes the factor π from the integration over theazimuthal angle φ .In Fig. 16, we compare the orbital mass weight distributionof the input oblate Abel model (top panels), with that of the best-fit axisymmetric Schwarzschild model (bottom panels). The three- integral mass weight distributions are quite similar, even in the pan-els with a relatively low mass content. The average fractional erroris ∼ %, and if we consider in each panel the mass weights abovethe mean value, which together contribute more than half of thetotal mass, the fractional error decreases to around ∼ %. Be-cause of possible strong point-to-point fluctuations as discussed in § ( E, I , I ) in the remaining two dimensions.Besides the total distribution, we show separately the contributionsfrom the NR component and the two opposite rotating SR compo-nents in the input oblate Abel model (see § ∼ % of the total mass),which also results in an underestimation of the NR component. This c (cid:13) , 1–33 ecovery orbital structure Figure 15.
The mean azimuthal motion h v φ i perpendicular to the meridional plane, normalised by σ RMS , for an oblate axisymmetric Abel model (left) andfor the best-fit axisymmetric Schwarzschild model (right). Parameters as in Fig. 9.
Figure 16.
The mass weight distribution for the input oblate Abel model (top) and for the best-fit axisymmetric Schwarzschild model (bottom). Parameters asin Fig. 12. The second integral of motion I = L z ≥ , where L z is the component of the angular momentum parallel to the symmetry z -axis. is reflected in the average absolute difference in mass as functionof E , which is ∼ . %. As for the triaxial case, the recovery for I and I is less good with average uncertainties of ∼ . % and ∼ . %, respectively.A similar good recovery was found by Krajnovi´c et al. (2005) for the case of a two-integral DF f ( E, L z ) , which implies anisotropic velocity distribution in the meridional plane. Thomas etal. (2004) showed that their independent axisymmetric numericalimplementation of Schwarzschild’s method is similarly able to re-cover an analytical f ( E, L z ) . Our results show that the orbital c (cid:13) , 1–33 G. van de Ven et al.
Figure 17.
The orbital mass weights for the input oblate Abel model (diamonds connected by solid curves) and for the best-fit axisymmetric Schwarzschildmodel (crosses connected by dotted curves). The parameters are as in in Fig. 13, except that rotation can only come from short-axis rotating (SR) components,for which the two directions of rotation are indicated separately. mass weight distribution that follows from a fully three-integral DF f ( E, L z , I ) can be recovered as well. We have extended the Abel models introduced by DL91 and gener-alised by MD99, and shown that, in addition to the intrinsic velocitymoments, the full LOSVD of these models can be calculated in astraightforward way. We have then used the Abel models to con-struct realistic axisymmetric and triaxial galaxy models to test theaccuracy of Schwarzschild’s orbit superposition method.Although Abel models have separable potentials with a centralcore and assume a specific functional form for the (three-integral)DF, they display a large range of shapes and their observables,which can be calculated easily, include many of the features seenin the kinematic maps of early-type galaxies. We have used anisochrone St¨ackel potential that in the axisymmetric limit reducesto the Kuzmin-Kutuzov model and becomes H´enon’s isochrone inthe spherical limit. Because of the simple form of the isochronepotential, the resulting Abel models are ideally suited to test nu-merical implementations of the Schwarzschild orbit superpositionmethod. The calculation of ∆ V , needed when comparing the or-bital mass weight distribution of the Schwarzschild models withthe three-integral DF of the Abel models, simplifies significantlyfor this case.Integral-field observations in principle provide the LOSVDas a function of position on the sky, so that it is a function L ( x ′ , y ′ , v z ′ ) that depends on three variables. The oblate axisym-metric and triaxial galaxy models we have constructed, have a DFwhich is a sum of Abel components f ( S ) = f ( − E + wI + uI ) with different values of the parameters w and u , so that the DF isa function of three variables as well, namely the integrals of mo-tion E , I and I . We have shown that the simulated integral-fieldobservables of these models are matched in detail by the best-fitSchwarzschild model. This does not automatically imply that theintrinsic velocity moments and the three-integral DF — which arenot directly fitted — are also correctly recovered.First consider three-integral oblate models, i.e., with a DF that is a function f ( E, L z , I ) . In the special case that a galaxy hap-pens to be well approximated by a two-integral DF f ( E, L z ) , thedensity ρ ( R, z ) uniquely determines the even part of f ( E, L z ) andthe mean streaming ρ h v φ i in the meridional plane fixes the partof f ( E, L z ) that is odd in L z (Dejonghe 1986). Ignoring non-uniqueness in the deprojection of the surface density Σ (Rybicki1987) and the mean streaming motion V on the plane of the sky,these two quantities define a two-integral DF completely. The ob-served velocity dispersion and higher-order velocity moments ofthe LOSVD then provide additional information, which for exam-ple can be used to constrain the inclination (e.g., Cappellari et al.2006). However, the reliability of the derived inclination, of course,depends on the correctness of the assumption of a two-integral DF.In the more realistic case of a three-integral DF f ( E, L z , I ) , sucha one-to-one relation with (the velocity moments of) the observedLOSVD L ( x ′ , y ′ , v z ′ ) has not been established. Nevertheless, weshowed that, given integral-field observations of the velocity mo-ments of the LOSVD (up to h ), recovery of the full three-integralDF is possible with Schwarzschild’s method, for the correct incli-nation and mass-to-light ratio.In the triaxial case, the DF is again a function of three inte-grals of motion, but the orbital structure in these models is sub-stantially richer than in the oblate axisymmetric models, with fourmajor orbit families, instead of only one. This introduces a funda-mental non-uniqueness in the recovery of the DF. Whereas in theoblate axisymmetric case ρ ( R, z ) uniquely defines the even partof f ( E, L z ) , in the (separable) triaxial case the density ρ ( x, y, z ) does not uniquely determine the even part of f ( E, I , I ) , althoughboth of these are functions of three variables (Hunter & de Zeeuw1992). It is (yet) unknown how much specification of L ( x ′ , y ′ , v z ′ ) can narrow down the range of possible DFs further, even ignor-ing the non-uniqueness caused by the required deprojection of thesurface brightness. Our results show that Schwarzschild’s methodrecovers the correct orbital mass weight distribution associated to f ( E, I , I ) . Given the very large freedom in the orbit choice forthis case, the modest resolution of our orbit library, and the result-ing approximations in the evaluation of the phase space volume,the agreement between the orbital mass weights found in § c (cid:13) , 1–33 ecovery orbital structure further by refining the sampling of the orbits and the regularisationof the orbital mass weights.Our analysis shows that it is clear that Abel models are usefulfor testing orbit-based modelling methods such as Schwarzschild’smethod. In particular the oblate limiting case with a Kuzmin-Kutuzov potential ( § ACKNOWLEDGEMENTS
We sincerely thank the referee for constructive comments and sug-gestions that improved the presentation, and Michele Cappellariand Anne-Marie Weijmans for a careful reading of an earlier ver-sion of the manuscript. GvdV acknowledges support provided byNASA through grant NNG04GL47G and through Hubble Fellow-ship grant HST-HF-01202.01-A awarded by the Space TelescopeScience Institute, which is operated by the Association of Uni-versities for Research in Astronomy, Inc., for NASA, under con-tract NAS 5-26555. RvdB acknowledges support by the Nether-lands Organization for Scientific Research (NWO) through grant614.000.301.
REFERENCES
Andrews D. F., Bickel P. J., Hampel F. R., Rogers W., TukeyJ., 1972, Robust Estimates of Locations: Survey and Advances.Princeton University, PrincetonArnold R., 1990, MNRAS, 244, 465Arnold R., de Zeeuw P. T., Hunter C., 1994, MNRAS, 271, 924Bacon R., et al. 2001, MNRAS, 326, 23Batsleer P., Dejonghe H., 1993, A&A, 271, 104Batsleer P., Dejonghe H., 1994, A&A, 287, 43Beers T. C., Flynn K., Gebhardt K., 1990, AJ, 100, 32Binney J., 1982, MNRAS, 201, 15Blinnikov S., Moessner R., 1998, A&AS, 130, 193Camm G. L., 1941, MNRAS, 101, 195Cappellari M., 2002, MNRAS, 333, 400 Cappellari M., Bacon R., Bureau M., Damen M. C., Davies R. L.,de Zeeuw P. T., Emsellem E., Falc´on-Barroso J., Krajnovi´c D.,Kuntschner H., McDermid R. M., Peletier R. F., Sarzi M., vanden Bosch R. C. E., van de Ven G., 2006, MNRAS, 366, 1126Cappellari M., Copin Y., 2003, MNRAS, 342, 345Cappellari M., Verolme E. K., van der Marel R. P., Verdoes KleijnG. A., Illingworth G. D., Franx M., Carollo C. M., de ZeeuwP. T., 2002, ApJ, 578, 787Carlson B. C., 1977, Special Functions of Applied Mathematics.Academic Press, New York San Francisco LondonCarollo C. M., de Zeeuw P. T., van der Marel R. P., 1995, MN-RAS, 276, 1131Chandrasekhar S., 1940, ApJ, 92, 441Copin Y., Cretton N., Emsellem E., 2004, A&A, 415, 889Cretton N., de Zeeuw P. T., van der Marel R. P., Rix H., 1999,ApJS, 124, 383Cretton N., Emsellem E., 2004, MNRAS, 347, L31Cretton N., Rix H.-W., de Zeeuw P. T., 2000, ApJ, 536, 319Cretton N., van den Bosch F. C., 1999, ApJ, 514, 704de Zeeuw P. T., 1985a, MNRAS, 216, 273de Zeeuw P. T., 1985b, MNRAS, 216, 599de Zeeuw P. T., Franx M., 1989, ApJ, 343, 617de Zeeuw P. T., Hunter C., Schwarzschild M., 1987, ApJ, 317, 607de Zeeuw P. T., Lynden-Bell D., 1985, MNRAS, 215, 713de Zeeuw P. T., Peletier R., Franx M., 1986, MNRAS, 221, 1001de Zeeuw P. T., Pfenniger D., 1988, MNRAS, 235, 949Dejonghe H., 1986, Phys. Rep., 133, 217Dejonghe H., de Zeeuw P. T., 1988, ApJ, 333, 90Dejonghe H., Laurent D., 1991, MNRAS, 252, 606 [DL91]Eddington A. S., 1915, MNRAS, 76, 37Eddington A. S., 1916, MNRAS, 76, 572Edgeworth F. Y., 1905, Cambridge Philos. Soc., 20, 36Emsellem E., et al. 2004, MNRAS, 352, 721Evans N. W., de Zeeuw P. T., 1992, MNRAS, 257, 152Franx M., 1988, MNRAS, 231, 285Gebhardt K., et al. 2000, AJ, 119, 1157Gebhardt K., et al. 2003, ApJ, 583, 92Gerhard O. E., 1993, MNRAS, 265, 213Gradshteyn I. S., Ryzhik I. M., 1994, Table of Integrals, Seriesand Products. Academic Press, San DiegoHenon M., 1959, Annales d’Astrophysique, 22, 126Hunter C., de Zeeuw P. T., 1992, ApJ, 389, 79Jeans J. H., 1915, MNRAS, 76, 70Kauffmann G., van den Bosch F., 2002, Scientific American, 286,36Krajnovi´c D., Cappellari M., Emsellem E., McDermid R. M., deZeeuw P. T., 2005, MNRAS, 357, 1113Kuzmin G. G., 1973, in Proc. All-Union Conf., Dynamics ofGalaxies and Clusters. ed. T. B. Omarov (Alma Ata: Akad. NaukKazakhskoj SSR), 71 (English transl. in IAU Symp. 127, Struc-ture and Dynamics of Ell. Galaxies, ed. P.T. de Zeeuw [Dor-drecht: Reidel], 553)Kuzmin G. G., Kutuzov S. A., 1962, Bull. Abastumani Astroph.Obs., 27, 82Lynden-Bell D., 1962, MNRAS, 124, 1Mathieu A., Dejonghe H., 1996, A&A, 314, 25Mathieu A., Dejonghe H., 1999, MNRAS, 303, 455 [MD99]Mathieu A., Dejonghe H., Hui X., 1996, A&A, 309, 30Merritt D., 1985, AJ, 90, 1027Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Osipkov L. P., 1979, Pis’ma v Astronomicheskii Zhurnal, 5, 77Pfenniger D., 1984, A&A, 134, 373 c (cid:13) , 1–33 G. van de Ven et al.
Primack J. R., 2004, in IAU Symp. 220: Dark Matter in Galaxies,eds. S. D. Ryder, D. J. Pisano, M. A. Walker, and K. C. Freeman,p. 53Richstone D. O., Tremaine S., 1984, ApJ, 286, 27Rix H., de Zeeuw P. T., Cretton N., van der Marel R. P., CarolloC. M., 1997, ApJ, 488, 702Rybicki G. B., 1987, in IAU Symp. 127: Structure and Dynamicsof Elliptical Galaxies, ed. P. T. de Zeeuw, p. 397Schwarzschild M., 1979, ApJ, 232, 236Schwarzschild M., 1982, ApJ, 263, 599Schwarzschild M., 1993, ApJ, 409, 563Statler T. S., 1987, ApJ, 321, 113Statler T. S., 1991, AJ, 102, 882Statler T. S., 1994, ApJ, 425, 458Teuben P., 1987, MNRAS, 227, 815Thomas J., Saglia R. P., Bender R., Thomas D., Gebhardt K.,Magorrian J., Corsini E. M., Wegner G., 2005, MNRAS, 360,1355Thomas J., Saglia R. P., Bender R., Thomas D., Gebhardt K.,Magorrian J., Richstone D., 2004, MNRAS, 353, 391Valluri M., Merritt D., Emsellem E., 2004, ApJ, 602, 66van Albada T. S., Bahcall J. N., Begeman K., Sancisi R., 1985,ApJ, 295, 305van de Ven G., Hunter C., Verolme E. K., de Zeeuw P. T., 2003,MNRAS, 342, 1056van de Ven G., van den Bosch R. C. E., Verolme E. K., de ZeeuwP. T., 2006, A&A, 445, 513van den Bosch R., de Zeeuw T., Gebhardt K., Noyola E., van deVen G., 2006, ApJ, 641, 852van den Bosch R. C. E., van de Ven G., Verolme E. K., CappellariM., de Zeeuw P. T., 2007, MNRAS, submitted [vdB07]van der Marel R. P., Cretton N., de Zeeuw P. T., Rix H., 1998,ApJ, 493, 613van der Marel R. P., de Zeeuw P., Rix H.-W., Quinlan G. D., 1997,Nature, 385, 610van der Marel R. P., Franx M., 1993, ApJ, 407, 525Vandervoort P. O., 1984, ApJ, 287, 475Verolme E. K., de Zeeuw P. T., 2002, MNRAS, 331, 959Verolme E. K., et al. 2002, MNRAS, 335, 517Verolme E. K., et al. 2003, in Galaxies and Chaos, eds. G. Con-topoulos and N. Voglis, LNP vol. 626, p. 279Weinacht J., 1924, Math. Ann., 91, 279
APPENDIX A: LIMITING CASES
When two or all three of the constants α , β or γ that define theconfocal ellipsoidal coordinate system are equal, the triaxial Abelmodels reduce to limiting cases with more symmetry and thus withfewer degrees of freedom. The oblate and prolate axisymmetriclimits are described in § A1 Elliptic disc potential
The two-dimensional analogues of the triaxial Abel models are theelliptic Abel discs with St¨ackel potential V S ( λ, µ ) = U [ λ, µ ] in confocal elliptic coordinates ( λ, µ ) . The relations with ( x, y ) fol-low from those in § z = 0 or equivalently ν = − γ .The two integrals of motion E and I are given by E = ` v x + v y ´ + U [ λ, µ ] , (A1) I = L z + ( α − β ) v x + ( α − β ) x U [ λ, µ, − α ] . A1.1 Velocity moments
Choosing the DF as f ( E, I ) = f ( S ) , with S = − E + w I , thevelocity moments can be evaluated as µ lm ( λ, µ ) = s l + m +2 h l +1 µ h m +1 λ × S max Z S min T lm [ S top ( λ, µ ) − S ] ( l + m ) / f ( S ) d S, (A2)with the terms h µ and h λ defined as h τ = 1 − ( τ + α ) w, τ = λ, µ. (A3)As in the general triaxial case, S min ≥ S lim , where the expressionof the latter is given along the w -axis ( u = 0 ) in Fig. 1. The acces-sible part of the ( E, I ) -integral space is now a triangle, the top ofwhich is S top ( λ, µ ) = − U [ λ, µ ] + w ( λ + α )( µ + α ) U [ λ, µ, − α ] .For the NR components we have that S max = S top ( λ, µ ) and T NR lm = B ( l +12 , m +12 ) . Of the two possible orbit families, the boxorbits have no net rotation and the tube orbits rotate around theaxis perpendicular to the disc (the z -axis). Since this is similar tothe short-axis tube orbits in the general triaxial case, we refer tothe rotating type as the SR type. This SR type reaches the regionof the accessible integral space (the triangle) for which v µ ≥ at µ = − α (or I ≥ ). Therefore, S max = S top ( λ, − α ) and T SR lm = 2 Z arcsin( √ a )0 sin l θ cos m θ dθ, (A4)where a is defined as a = ( λ + α ) h µ [ S top ( λ, − α ) − S ]( λ − µ ) h ( − α ) [ S top ( λ, µ ) − S ] . (A5)The integral (A4) can be evaluated in terms of elementary functions(e.g., Gradshteyn & Ryzhik 1994, relations 2.513 on p.160–162).The NR velocity moments µ NR lm ( λ, µ ) vanish when either l or m isodd, and the SR velocity moments µ SR lm ( λ, µ ) vanish when l is odd.The latter should be multiplied with ( − m for net rotation in theopposite direction.The matrix Q , which transforms the velocity components ( v λ , v µ , v ν ) to ( v x , v y , v z ) , is that for the prolate case given ineq. (6.4), but with χ = 0 substituted. The sign matrix S , projec-tion matrix P and rotation matrix R are the same as for the triaxialcase given in respectively eqs (3.3)–(3.5). The polar angle is theinclination, ϑ = i , and the azimuthal angle ϕ the orientation of theinfinitesimally thin disc ( γ = 0 ) in the plane z = 0 . In the expres-sion (3.1) for the misalignment angle ψ , the triaxiality parameterthus reduces to T = 1 − β/α , with < β < α , bracketing thelimiting cases of a needle and a circular disc. A1.2 Line-of-sight velocity distribution
Starting with the expression for the stellar (surface) mass density Σ ⋆ ( x ′ , y ′ ) = µ ( λ, µ ) from eq. (A2), we derive the Abel LOSVD c (cid:13) , 1–33 ecovery orbital structure for the elliptic disc in a similar way as in § . The cross section of the unit sphere with a plane, reducesto the cross section of the unit circle X + Y = 1 with a line e X + e Y = v z ′ /g ( S ) . Here, the variable X is defined as X = h v λ g ( S ) √ h λ , with g ( S ) = h s S top ( λ, µ ) − S ] h λ h µ , (A6)and the variable Y follows by interchanging λ ↔ µ . The coeffi-cients of the line are e = √ h λ M /h and e = p h µ M /h . Wecan write the normalisation h = h λ M + h µ M , explicitly as h = sin i ˘ [1 − ( λ + α ) w ] ( C cos ϕ + D sin ϕ ) +[1 − ( µ + α ) w ] ( C sin ϕ − D cos ϕ ) ¯ , (A7)with C and D defined in eq. (6.5). In this way we can write Σ ⋆ as a double integral over S and v z ′ , so that at a given line-of-sightvelocity the LOSVD becomes L ( x ′ , y ′ , v z ′ ) = 1 h S up ( v z ′ ) Z S min f ( S ) p G ( v z ′ ) − S ] ∆( v z ′ , S ) d S, = p G ( v z ′ ) − S min ] h η up Z f ( S ) ∆( v z ′ , S ) sin η d η, (A8)which vanishes when | v z ′ | exceeds the ’terminal velocity’ v t = g ( S min ) . The second expression removes the possible singularityat the upper limit of S , given by S up = min[ G ( v z ′ ) , S max ] , with G ( v z ′ ) = S top ( λ, µ ) − h λ h µ v z ′ / (2 h ) . (A9)Hence, η up is given by sin η up = [ S up − S min ] / [ G ( v z ′ ) − S min ] .Depending on the integral space accessible by the orbits, the valueof ∆( v z ′ , S ) is either zero, one or two.For the NR components, S max = S top ( λ, µ ) , and, since thefull integral space is accessible, ∆ NR = 2 , independent of S and v z ′ . In the case of a basis function f δ ( S ) as defined in eq. (2.24),the integral over S can be evaluated explicitly resulting in L NR δ = 4 δ +1 B ( δ + 1 , δ + 1) √ − S min ) δ h [ G ( v z ′ ) − S min ] δ + . (A10)For SR components with maximum streaming, both v µ ≥ at µ = − α and v µ ≥ , which is equivalently to ≤ Y ≤ √ a ,where a is defined in eq. (A5). The intersection of the above unitcircle and line provides the following two solutions Y ± = e [ v z ′ /g ( S )] ± e p − [ v z ′ /g ( S )] . (A11)Given the values of v z ′ and S , ∆ SR is thus equal to 0, 1 or 2 if forrespectively none, one or both of the solutions ≤ Y ± ≤ √ a .The expression for h in eq. (A7) shows that the LOSVD ineq. (A8) is inversely proportional to sin i . For face-on viewing atinclination i = 0 ◦ the LOSVD reduces to Σ ⋆ ( x ′ , y ′ ) δ ( v z ′ ) . Be-cause the velocity perpendicular to the disc v z = 0 , the face-onLOSVD is zero at all line-of-sight velocities, except at v z ′ = 0 when it equals the surface mass density. For edge-on viewing at in-clination i = 90 ◦ , the LOSVD follows upon substituting y ′ = 0 ineq. (A8) and integrating over the line-of-sight z ′ . For i < ◦ , the Alternatively, one can invert the relations S = S top ( λ, µ ) − ( h µ v λ + h λ v µ ) and v z ′ = M v λ + M v µ to find the Jacobian to transformfrom the coordinates ( v λ , v µ ) to ( S, v z ′ ) . Leaving out the integral over v z ′ yields the same expression for the LOSVD as in eq. (A8). latter integration is not needed, since at each position ( x ′ , y ′ ) on theplane of the sky there is only a single (unique) point along z ′ whereit intersects the infinitesimally thin disc. The edge-on LOSVD andalso Σ ⋆ are thus spatially only one-dimensional functions of x ′ , andvanish for non-zero y ′ -values.Further information on elliptic St¨ackel discs can be found inTeuben (1987), de Zeeuw, Hunter & Schwarzschild (1987), andEvans & de Zeeuw (1992). A2 Large distance limit
At large radii, λ → r ≫ − α , so that the confocal ellipsoidalcoordinates of § ( r, µ, ν ) , with r the usual distance to the origin, i.e., r = x + y + z , and µ and ν angular coordinates on the sphere. In these coordinates theSt¨ackel potential is of the form V S ( r, µ, ν ) = V ( r ) + U [ µ, ν ] /r ,where V ( r ) is an arbitrary smooth function of r . The correspondingintegrals of motion are given by E = ` v x + v y + v z ´ + V S ( r, µ, ν ) ,I = T L y + L z + ( α − β ) x r U [ µ, ν, − α ] , (A12) I = L x + (1 − T ) L y + ( γ − β ) z r U [ µ, ν, − γ ] . With the choice (2.9) for the DF, the expression for the velocitymoments becomes µ lmn ( r, µ, ν ) = 1 r m + n +2 s l + m + n +3 F m +1 ν F n +1 µ × S max Z S min T lmn [ S top ( r, µ, ν ) − S ] ( l + m + n +1) / f ( S ) d S, (A13)where F ν and F µ are defined as F τ = 1 r + ( τ + α ) w − ( τ + γ ) uγ − α , τ = µ, ν. (A14)As in the general triaxial case, S min ≥ S lim , where S lim can beobtained from Fig. 1. The expressions of S max and T lmn for theNR, LR and SR types are those given in §§ S top ( λ, µ, ν ) (eq. 2.14) replaced by S top ( r, µ, ν ) = − V S ( r, µ, ν ) − w ( µ + α )( ν + α ) γ − α U [ µ, ν, − α ] − u ( µ + γ )( ν + γ ) α − γ U [ µ, ν, − γ ] , (A15)and the parameters a and b (2.17) reduce to a = S top ( r, µ, − β ) − SS top ( r, µ, ν ) − S , (A16) b = ( µ + β ) F ν [ S top ( r, µ, − β ) − S ]( µ − ν ) F ( − β ) [ S top ( r, µ, ν ) − S ] , which by interchanging ν ↔ µ become a and b , and in turn a and b follow by β ↔ α .In the conversion to observables described in §
3, in the ma-trix Q , which transforms the velocity components ( v r , v µ , v ν ) to ( v x , v y , v z ) , all terms λ + σ ( σ = − α, − β, − γ, µ, ν ) cancel out(cf. eq. 25 of Statler 1994). The expression for the LOSVD followsfrom that of the triaxial case in eq. (3.27) by substituting H µν = 1 , H νλ = r F ν , F λµ = r F µ and S top ( λ, µ, ν ) = S top ( r, µ, ν ) .Suppose now that at large radii r , the function V ( r ) in the c (cid:13) , 1–33 G. van de Ven et al.
St¨ackel potential vanishes and we keep in the above expressionsonly the dominant terms. In this case, F µ , F ν and S top reduce tofunctions of µ and ν only. As a result, the velocity moments (A13)are independent of r , except for the prefactor /r m + n +2 , and there-fore are scale-free. Once we have calculated the velocity momentsat a radius r , those at radius r ′ = qr , with q a constant, follow bya simple scaling, µ lmn ( r ′ , µ, ν ) = µ lmn ( r, µ, ν ) /q m + n +2 . Thesame holds true for the line-of-sight velocity moments µ k ( r, µ, ν ) ,but not for the LOSVD. A3 Spherical potential
When α = β = γ , both µ and ν loose their meaning and wereplace them by the customary polar angle θ and azimuthal angle φ .The expressions for the Abel models in these spherical coordinates ( r, θ, φ ) follow in a straightforward way from those in § A2 for thelarge distance limit in conical coordinates ( r, µ, ν ) .The St¨ackel potential V S = V ( r ) is spherically symmetric.The expressions for the integrals of motion follow from (A12),where for I and I the right-most terms vanish. The triaxialityparameter T is now a free parameter, so that, together with the pa-rameters w and u , we can rewrite S = − E + w I + u I as S = − E + uL x + [(1 − T ) u + T w ] L y + wL z . (A17)This means that with the choice (2.9) for the DF, we cover themost general homogeneous quadratic form in the velocities that isallowed by the integrals of motion in a spherical symmetric po-tential, i.e., the energy E and all three components of the angularmomentum vector L (cf. DL91). These include the models consid-ered by Osipkov (1979) and Merritt (1985) with the DF of the from f ( − E ± L /r a ) and those studied by Arnold (1990) with a moregeneral DF of the form f ( − E ± L /r a ± L z /r b ) . These modelsfollow by setting u = w = ± /r a , and by taking u = ± /r a , w = u ± /r b and T = 0 , respectively. A3.1 Velocity moments
The velocity moments follow from eq. (A13), with F τ = 1 r − ( w + u )+ ( w − u ) h cos θ + T (sin θ sin φ − ± √ Λ i , (A18)where the positive and negative sign are for F µ and F ν , and Λ = ˆ sin θ + T (sin θ sin φ − ˜ + 4 T sin θ cos θ sin φ. (A19)Taking α = β = γ in Fig. 1, we see that the boundaries on w and u both vanish. The separatrices L and L , defined in eq. (2.13),reduce to the negative w -axis and the line w = u , respectively.Furthermore, S max = S top = − V ( r ) , and for T lmn we use theexpression (2.15). The resulting velocity moments µ lmn ( r, θ, φ ) ,which are in general not spherically symmetric, vanish when either l , m or n is odd.The latter implies no net rotation, which is the case whenthe (conserved) angular momentum vectors L for the orbits arerandomly oriented. We may introduce net rotation by assumingthat (a fraction of) the orbits have a preferred sense of rotationaround an angular momentum vector L that points in a specificdirection given by θ and φ . Using the projection matrix P in eq. (3.4) with ϑ = θ and ϕ = φ , we transform to the co-ordinate system ( r ′ = r, θ ′ , φ ′ ) , in which L is aligned withthe z ′ -axis. If we next set the DF to zero for L z ′ < , we find µ ′ lmn ( r, θ ′ , φ ′ ) = µ lmn ( r, θ ′ , φ ′ ) , which does still vanish when l or m is odd, but is non-zero when n is odd, resulting in maxi-mum streaming around the z ′ -axis, and multiplication with ( − n for opposite direction of rotation. With the inverse of the projectionmatrix, we can then transform these velocity moments to the origi-nal coordinates system ( r, θ, φ ) . In this way, we can build sphericalAbel models, which in addition to a non-rotating part consist of acomponent or several components with a preferred rotation axis.Mathieu, Dejonghe & Hui (1996) used this approach to construct adynamical model of Centaurus A, with a spherical potential, but tri-axial luminosity density and DF components with rotation aroundthe apparent long and short axis.From the customary definition of the spherical coordinate sys-tem, x = r sin θ cos φ , y = r sin θ sin φ and z = r cos θ , it followsdirectly that the matrix Q , which transforms the velocity compo-nents ( v r , v θ , v φ ) to ( v x , v y , v z ) , is given by Q = sin θ cos φ cos θ cos φ − sin φ sin θ sin φ cos θ sin φ cos φ cos θ − sin θ . (A20)In case the orbits have no preferred sense of rotation, we may setthe viewing angles ϑ = ϕ = 0 without loss of generality, so thatwith ψ = 0 from eq. (3.1), ( x ′ , y ′ ) = ( y, − x ) on the plane ofthe sky and z ′ = z along the line-of-sight, and similarly for theCartesian velocity components. A3.2 Line-of-sight velocity distribution
There is no obvious further simplification of the LOSVD for ro-tating components. For the NR components, the LOSVD followsfrom eq. (3.27) with S max = S top = − V ( r ) and ∆ ξ ′ = 2 π ,or from eq. (3.28) after substituting the basis function f δ ( S ) fromeq. (2.24). Since the line-of-sight velocity v z ′ = sgn( z )[cos θ v r − sin θ v θ ] , it follows that M = cos θ , M = sin θ and M = 0 in eq (3.19) for h . Moreover, H µν = 1 and H τλ = r F τ ( τ = µ, ν ), with F τ given in eq. (A18).Carollo et al. (1995) compute the LOSVD for Osipkov-Merrittmodels with f ( − E − L / (2 r a )) , which is a special case of the AbelDF f ( S ) , that follows from eq. (A17) by setting u = w = − /r a .Substituting the latter in eq. (A18), we find that H λµ = H λν =( r a + r ) /r a , so that h = ( r a + r )( r a + r − R ′ ) /r a , (A21) G ( v z ′ ) = − V ( r ) − r a + r r a + r − R ′ v z ′ , (A22)with radius R ′ = r sin θ on the plane of the sky. After substitutionin eq. (3.27), and transforming the integral over d z ′ to d r , we findthe following LOSVD L ( R ′ , v z ′ ) = 4 π ∞ Z R ′ rh √ r − R ′ G ( v z ′ ) Z S min f ( S ) d S d r. (A23)This is the same as the (unnormalised) velocity profile in eq. (27)of Carollo et al. (1995), with Φ ∞ , the lower limit of their (relative)potential Φ( r ) = − V ( r ) , equal to S min , which from Fig. 1 in thiscase has S lim = 0 as lower limit. Their function g ( r, R ′ ) and upperlimit Q max in eqs (25) and (26), are equivalent to respectively the c (cid:13) , 1–33 ecovery orbital structure Table B1.
The function M for odd s . s i j M ( s, i, j ; a, b, φ ) φ (4 − a − b ) φ + ( b − a ) sin 2 φ − φ − sin 2 φ − φ + sin 2 φ (24 − a − b + 3 a + 3 b + 2 ab ) φ + ( b − a )(3 − a − b ) sin 2 φ + ( b − a ) sin 4 φ − (6 − a − b ) φ − (3 − a ) sin 2 φ − ( b − a ) sin 4 φ − (6 − b − a ) φ + (3 − b ) sin 2 φ + ( b − a ) sin 4 φ φ + sin 2 φ + sin 4 φ φ − sin 4 φ φ − sin 2 φ + sin 4 φ inverse of h and G ( v z ′ ) in eqs (A21) and (A22) above. The well-known isotropic case follows upon taken the limit r a → ∞ , so that f ( S ) → f ( − E ) , h → and G ( v z ′ ) → − V ( r ) − v z ′ / . APPENDIX B: THE FUNCTION M The function M that appears in the velocity moments of the rotat-ing Abel components is defined as M ( s, i, j ; a, b, φ ) = φ Z „ ∂∂a « i „ ∂∂b « j − q [1 − p ( θ )] s +1 p ( θ ) d θ, (B1)with p ( θ ) ≡ a cos θ + b sin θ . For odd s , corresponding to oddvelocity moments, the integral can be evaluated in a straightforwardway in terms of elementary functions. In Table B1, we give theresulting expressions for s = 1 , , .For even s , the integral can be evaluated in terms of the (in-complete) elliptic integrals. To simplify the numerical evaluationwe use Carlson’s (1977) symmetrical forms R F , R D and R J (forthe relations between both forms see e.g. de Zeeuw & Pfenniger1988). In Table B2, we give the expressions for s = 0 , , , wherewe have introduced the following quantities based on these sym-metric elliptic integrals F = √ − a sin φa R F (cos φ, ∆ , ,D = sin φ √ − a R D (cos φ, ∆ , , (B2) J = ( b − a ) sin φ a √ − a R J (cos φ, ∆ , , p ( φ ) a ) , with ∆ = [1 − p ( φ )] / (1 − a ) , and we have defined the terms A = 1 √ ab arctan r ba tan φ ! ,P = sin φ cos φ p − p ( φ ) , (B3) Q = sin φ cos φ − p − p ( φ ) p ( φ ) . In Fig. B1, we show the M ( s, i, j ; a, b, φ ) as function of φ for thecase that a = 0 . and b = 0 . , up to order s = 5 .We now consider some special cases. When either a or b iszero, the corresponding velocity moments vanish (eqs 2.16 and 2.19), and when a i > b i the arguments of the function M areinterchanged (eqs 2.18, 2.20 and 2.21). This means we only haveto consider the range < a ≤ b , together with < φ ≤ π/ ,since M vanishes when φ = 0 .When a = b , it follows that p ( θ ) = a in eq. (B1), so we canseparate M ( s, i, j ; a, a, φ ) = M ( s, i, j ; a ) M ( i, j ; φ ) , where M ( s, i, j ; a ) = d i + j d a i + j − p (1 − a ) s +1 a (B4) M ( i, j ; φ ) = Z φ cos i θ sin j θ φ Z cos i θ sin j θ d θ. For a = 1 , the expression for M simplifies to ( − i + j ( i + j )! .The integral in the expression for M can be evaluated explicitlyusing e.g. the relations 2.513 of Gradshteyn & Ryzhik (1994). For φ = π/ , it reduces to the beta function B ( i + 1 / , j + 1 / .When a < b = 1 , the elliptic integrals become elementary, sothat the quantities F , D and J in eq. (B2) reduce to F = √ − aa ln » tan „ π φ «– ,D = a − a F − sin φ √ − a , (B5) J = F − √ a arctan r − aa sin φ ! . Although F diverges when φ → π/ , substitution of these reducedquantities in the expressions of M for even s (Table B2), showsthat all terms with F cancel. For φ = π/ , the function M is thuseverywhere finite, with A = π/ (2 √ ab ) and P = Q = 0 . APPENDIX C: EDGEWORTH EXPANSION
For the (re)construction of the LOSVD from its true line-of-sightvelocity moments, one can use the well-known Gram-Charlier se-ries, the terms of which are simple functions of the true moments(see e.g. Appendix B2 of van der Marel & Franx 1993), but it haspoor convergence properties. The terms in the Edgeworth (1905)expansion are also directly related to the true moments, but sinceit is a true asymptotic expansion its accuracy is controlled, so that,unlike the Gauss-Hermite and Gram-Charlier expansions, conver-gence plays no role (see Blinnikov & Moessner 1998 for a compar-ison between the expansions and for further references).The Edgeworth expansion of the LOSVD up to order N isgiven by L ED N ( v ) = Σ e − w √ πσ " N X n =3 D n , (C1)with w = ( v − V ) /σ and D n = X { l i − } H n +2( l − ( w ) n Y i =3 l i − ! „ d i i ! « l i − . (C2)The Hermite polynomials H m are related to those defined by vander Marel & Franx (1993) as H m ( w ) = √ m ! H m ( w/ √ . Wehave defined l = P n − j =1 l j , where the sets { l j } are the non-negativeinteger solutions of the Diophantine equation l j + 2 l j + · · · + ( n − l n − = n − , n ≥ , (C3) c (cid:13) , 1–33 G. van de Ven et al.
Table B2.
The function M for even s . s i j M ( s, i, j ; a, b, φ ) A − F + J A − (1 − a ) F − ( b − a ) D + J − a [ A + Q − (1 + a ) F + (1 − a ) D + J ] − b [ A − Q − F − (1 − b ) D + J ] A + ( b − a ) P − (2 a + ab − a + 3) F + (2 a + 2 b − b − a ) D + J − a ˆ A + a P + Q − (1 + 2 a )(1 − a ) F + (2 a − a − ab + 1) D + J ˜ − b ˆ A − b P − Q − (1 − ab ) F − (2 b − b − ab + 1) D + J ˜ a n A + a p ( φ ) − ab b − a ) p ( φ ) P + a cos φ +3 b sin φ p ( φ ) Q + a − a b +4 a +3 a − ab − b b − a ) F − (2 a +5 a − ab − b )(1 − a )3( b − a ) D + J o ab n A + ab − ab p ( φ )( b − a ) p ( φ ) P + b sin φ − a cos φp ( φ ) Q + a b − ab + a − bb − a F + a b + ab − ab + a + bb − a D + J o b n A + b p ( φ ) − ab b − a ) p ( φ ) P − a cos φ +5 b sin φ p ( φ ) Q − b − a − ab + ab b − a ) F − (2 b +5 b − ab − a )(1 − b )3( b − a ) D + J o Figure B1.
The function M ( s, i, j ; a, b, φ ) defined in eq. (B1) plotted against φ , for a = 0 . and b = 0 . , up to order s = 5 . The curves in the left twopanels are for odd values of s corresponding to the odd velocity moments, whereas the curves in the right two panels are for even values of s . The indices ofthe labels M sij refer to the first three parameters of the function M . Substituting these solutions, we find up to order N = 5 L ED5 ( v ) = Σ e − w √ πσ " H ( w ) d
3! + H ( w ) d H ( w ) 12 „ d « + H ( w ) d H ( w ) d d
4! + H ( w ) 16 „ d « . (C4)The lower-order moments Σ , V and σ are equivalent to those ineq. (3.41), while the higher-order moments d i ( i ≥ ) are cumu-lants of the true moments d i = i ! σ n X { l k } ( − l − ( l − i Y k =1 l k ! “ µ k k ! ” l k , (C5)so that d = ξ , d = ξ − , and d = ξ − ξ . (C6)The central moments ξ (skewness), ξ (kurtosis) and ξ are related to the true moments respectively as ( µ σ ) ξ = µ µ − µ µ µ + 2 µ , (C7) ( µ σ ) ξ = µ µ − µ µ µ + 6 µ µ µ − µ , (C8) ( µ σ ) ξ = µ µ − µ µ µ + 10 µ µ µ − µ µ µ + 4 µ . (C9)Substituting the line-of-sight true moments µ k for k = 0 , . . . , K ,we can compute L ED K ( v ) at each position on the plane of the sky.In Fig. C1, we show an example of a LOSVD (black solidline) computed directly via eq. (3.27) for the triaxial Abel modelconstructed in § c (cid:13) , 1–33 ecovery orbital structure Figure C1.
Line-of-sight velocity distribution (LOSVD) of the triaxial Abelmodel presented in § than N = 5 are needed to accurately determine the wings of theLOSVD. Nevertheless, this comparison is important to show thecorrectness of both (independent) approaches, and that the Edge-worth expansion provides a reliable and efficient way to reconstructthe LOSVD (and obtain Gauss-Hermite moments) from true mo-ments that are in general (numerically) easier to compute than thefull LOSVD.This paper has been typeset from a TEX/ L A TEX file prepared by theauthor. c (cid:13) , 1–33 G. van de Ven et al.
Figure C2.
Maps of the surface mass density ( Σ ; in M ⊙ pc − ), mean line-of-sight velocity V and dispersion σ (both in km s − ), and higher orderGauss-Hermite moments h and h , of the triaxial Abel model constructed in § (cid:13)000