Recovery of the internal orbital structure of galaxies
We construct axisymmetric and triaxial galaxy models with a phase-space distribution function that depends on linear combinations of the three exact integrals of motion for a separable potential. These Abel models, first introduced by Dejonghe & Laurent and subsequently extended by Mathieu & Dejonghe, are the axisymmetric and triaxial generalisations of the well-known spherical Osipkov–Merritt models. We show that the density and higher order velocity moments, as well as the line-of-sight velocity distribution (LOSVD) of these models can be calculated efficiently and that they capture much of the rich internal dynamics of early-type galaxies. We build a triaxial and oblate axisymmetric galaxy model with projected kinematics that mimic the two-dimensional kinematic observations that are obtained with integral-field spectrographs such as SAURON. We fit the simulated observations with axisymmetric and triaxial dynamical models constructed with our numerical implementation of Schwarzschild’s orbit-superposition method. We find that Schwarzschild’s method is able to recover the internal dynamics and three-integral distribution function of realistic models of early-type galaxies.
keywords: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 an elliptical or lenticular galaxy is completely described by its distribution function (DF) in the six-dimensional phase space of positions and velocities. The recovery of the DF from observations is difficult, as for galaxies other than our own, we can usually only measure the projected surface brightness and the line-of-sight velocity distribution (LOSVD) of the integrated light as a function of position on the plane of the sky. Moreover, we generally do not know the intrinsic shape of the galaxy, nor the viewing direction, or the contribution to the gravitational potential provided by a super massive central black hole and/or an extended halo of dark matter. By Jeans (1915) theorem, the DF is a function of the isolating integrals of motion admitted by the potential, but it is not evident how to take advantage of this property other than for the limiting case of spherical systems. Orbits in axisymmetric geometry have two exact integrals of motion, the energy and the angular momentum component parallel to the symmetry -axis, but the third effective or non-classical integral obeyed by all regular orbits is generally not known in closed form. In stationary triaxial geometry is conserved, but regular orbits now have two additional effective integrals of motion, and , which are not known explicitly.
Schwarzschild (1979, 1982) devised a numerical method which sidesteps our ignorance about the non-classical integrals of motion. It allows for an arbitrary gravitational potential, which may include contributions from dark components, integrates the equations of motion for a representative library of orbits, computes the density distribution of each orbit, and then determines the orbital weights such that the combined orbital densities reproduce the density of the system. The best-fitting orbital weights represent the 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 implementations of Schwarzschild’s method for axisymmetric geometry which fit the projected surface brightness and line-of-sight velocity distributions 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). Applications include the determination of central black hole masses (see also van der Marel et al. 1997; Cretton & van den Bosch 1999; Verolme et al. 2002; Cappellari et al. 2002; Gebhardt et al. 2003; Copin, Cretton & Emsellem 2004), accurate global dynamical mass-to-light ratios (Cappellari et al. 2006), as well as dark matter profiles as a function of radius (Cretton, Rix & de Zeeuw 2000; Thomas et al. 2005), and recovery of the DF (Krajnović et al. 2005). Van de Ven et al. (2006) and van den Bosch et al. (2006) included proper motion measurements in order to model nearby globular clusters, and determine their distance, inclination as well as mass-to-light ratio as function of radius. Finally, Verolme et al. (2003) and the companion paper van den Bosch et al. (2007, hereafter vdB07) describe an extension to triaxial geometry that includes all line-of-sight kinematics.
Although Schwarzschild models have significantly increased our understanding of the dynamical structure and evolution of early-type galaxies, questions remain about the uniqueness and the accuracy with which they are able to recover the global parameters as well as the internal dynamics of these galaxies. Many tests have been done to establish how the axisymmetric code recovers known input models, but these generally have been limited to spherical geometry or to an input axisymmetric DF that is a function of and 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ć et al. 2005).
One could construct a numerical galaxy model with Schwarzschild’s method itself, compute the observables, and then use these as input for the code and determine how well it recovers the input model. This is useful, but does not provide a fully independent test of the software. An alternative is to consider the special family of models with gravitational potential of Stäckel form, for which all three integrals of motion are exact and known explicitly. These separable potentials have a core rather than a central cusp, so the corresponding models cannot include a central black hole, and are inadequate for describing galactic nuclei. However, they can be constructed for a large range of axis ratios (Statler 1987), and their observed kinematic properties are as rich as those seen in the main body of early-type galaxies (Statler 1991, 1994; Arnold, de Zeeuw & Hunter 1994).
A small number of analytic DFs have been constructed for triaxial separable models. The ‘thin-orbit’ models (Hunter & de Zeeuw 1992) have the maximum possible streaming motions, but their DF contains delta functions, and they are therefore not particularly useful for a test of general-purpose numerical machinery. Dejonghe & Laurent (1991, hereafter DL91) constructed separable triaxial models in which the DF depends on a single parameter , which is a linear combination of the three exact integrals , and admitted by these potentials, and is quadratic in the velocity components. For a given radial density profile, the DF follows by simple inversion of an Abel integral equation. These so-called Abel models have no net mean streaming motions, and are the axisymmetric and triaxial generalisations of the 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 including two families of DF components with net internal mean motions around the long and the short axis, respectively, and compared the resulting models with observations of Centaurus A. Although the Abel character of the non-rotating components is no longer conserved, the expressions for the velocity moments in these more general models can still be evaluated in a straightforward way. When the entire DF depends on the same single variable the famous ellipsoidal hypothesis (e.g., Eddington 1915; Chandrasekhar 1940) applies, so that self-consistency is only possible in the spherical case (Eddington 1916; Camm 1941). This does not hold for Abel models with a DF that is a sum of components for which the variable has different values of the parameters and . Such multi-component Abel models can provide (nearly) self-consistent models with a large variety of shapes and dynamics.
Here, we show that for Abel models, in addition to the velocity moments, the full LOSVD can be calculated in a simple way. Next, we construct axisymmetric and triaxial Abel models to test our numerical implementation of Schwarzschild’s method. We assume a convenient form for the gravitational potential, and construct the DF that reproduces a realistic surface brightness distribution. We compute the LOSVDs of the models and derive two-dimensional maps of the resulting kinematics. We show that, despite the simple form of the DF, these models display the large variety of features observed in early-type galaxies with integral-field spectrographs such as SAURON (Emsellem et al. 2004). By fitting axisymmetric and triaxial three-integral Schwarzschild models to the simulated observables we find that Schwarzschild’s method is able to recover the internal dynamics and three-integral DF of early-type galaxies. In this paper we fix the mass-to-light ratio and viewing direction to those of the Abel models, while in our companion paper vdB07 we investigate how well these global parameters can be determined by Schwarzschild’s method, along with a full description of our numerical implementation in triaxial geometry.
This paper is organised as follows. In Section 2 we summarise the properties of the triaxial Abel models of DL91 and MD99 and present the intrinsic velocity moments in a form which facilitates their numerical implementation. We describe the conversion to observables in Section 3, including the computation of the LOSVD. In Section 4 we construct a specific triaxial galaxy model and in Section 5 we fit the simulated observables with our triaxial Schwarzschild models to investigate how well the intrinsic moments and three-integral DF are recovered. In Section 6 we consider Abel models in the axisymmetric limit and construct a three-integral oblate galaxy model to test our axisymmetric implementation of Schwarzschild’s method. We summarise our conclusions in Section 7. In Appendix A, we describe the simpler Abel models for the elliptic disc, large distance and spherical limit, and link them to the classical Osipkov–Merritt solutions for spheres. Readers who are mainly interested in the tests of the Schwarzschild method may skip Sections 2 – 4 and 6.1 – 6.3.
2 Triaxial Abel models
The triaxial Abel models introduced by DL91 have gravitational potentials of Stäckel form, for which the equations of motion separate in confocal ellipsoidal coordinates. We briefly describe these potentials, and refer for further details to de Zeeuw (1985a). We then make a specific choice for the DF, for which the velocity moments simplify.
2.1 Stäckel potentials
We define confocal ellipsoidal coordinates () as the three roots for of
with () the usual Cartesian coordinates, and with constants and such that . From the inverse relations
and similarly for and by cyclic permutation of , it follows that a combination () generally corresponds to eight different points (). In these coordinates, the Stäckel potentials have the following form (Weinacht 1924)
where is an arbitrary smooth function . The right-hand side of eq. (2.3) can be recognised as the second order divided difference of . Henceforth, we denote it with the customary expression , which is symmetric in its arguments (see Hunter & de Zeeuw 1992, eqs 2.1–2.3, 2.13 and 2.14). Addition of a linear function of to does not change .
The density that corresponds to can be found from Poisson’s equation
or alternatively by application of Kuzmin’s (1973) formula (see de Zeeuw 1985b). This formula shows that, once we have chosen the confocal coordinate system and the density along the short axis, the mass model is fixed everywhere by the requirement of separability111A third method for the calculation of the density is to use , where the fifth-order divided difference is of the function with and defines the potential as in eq. (2.3). This result was obtained by Hunter in 1989 (priv. comm.) and by Mathieu & Dejonghe (1996). Similar expressions exist for the related families of potential-density pairs introduced in de Zeeuw & Pfenniger (1988).. For centrally concentrated mass models, has the -axis as long-axis and the -axis as short-axis. In most cases this is also true for the associated density (de Zeeuw, Peletier & Franx 1986).
2.2 Orbital structure
The Hamilton-Jacobi equation separates in for the potentials (2.3), so that every orbit has three exact integrals of motion (cf. de Zeeuw & Lynden-Bell 1985)
where , and are the velocity components in the Cartesian coordinate system, and , the component of the angular momentum vector parallel to the -axis. The other two components, and , follow by cyclic permutation of and . Furthermore, is a triaxiality parameter defined as
and is the third-order divided difference of . All models for which have a similar orbital structure and support four families of regular orbits: boxes with no net rotation, inner and outer long-axis tubes with net rotation around the -axis, and short-axis tubes with net rotation around the -axis (Kuzmin 1973; de Zeeuw 1985a; Hunter & de Zeeuw 1992).
According to Jeans (1915) theorem the phase-space distribution function (DF) is a function of the isolating integrals of motion (cf. Lynden-Bell 1962; Binney 1982). The velocity moments of the DF are defined as
where , and are non-negative integers, and , and are the velocity components in the confocal ellipsoidal coordinate system. Many of the velocity moments vanish due to the symmetry of the orbits in these coordinates. The zeroth-order velocity moment is the mass density that corresponds to the DF
In self-consistent models, must equal given in eq. (2.4), the mass density that is related to the potential by Poisson’s equation.
2.3 Abel distribution function
Following DL91, we choose the DF to be a function of the three integrals of motion , and as given in eq. (2.2) through one variable
and and are two parameters222In contrast with DL91 and MD99, we choose and , consistent with e.g. de Zeeuw (1985a).. This choice for the DF is equivalent to the celebrated ellipsoidal hypothesis (e.g., Eddington 1915; Chandrasekhar 1940). Self-consistency is only possible in the spherical case (Eddington 1916; Camm 1941). On the other hand, these DFs can produce realistic (luminous) mass densities , which differ from the (total) mass density , as in galaxies with dark matter (see also § 2.4 below when we combine DFs of the form [2.9] with different values for and .)
DL91 and MD99 divided the DF into three types of components. The non-rotating (NR) type is made of box orbits and tube orbits with both senses of rotation populated equally. The two rotating types, LR and SR, consist of tube orbits, and have net rotation around either the long axis or the short axis.
2.3.1 Velocity moments
Due to the choice (2.9) of the DF, the general expression (2.7) for the velocity moments can be simplified, as shown by DL91 for the non-rotating components and by MD99 for the rotating components. We recast their expressions into a different form to facilitate the numerical implementation. The resulting velocity moments are given by
and set to zero at positions for which . The terms , and in the square root in front of the integral are defined as
with . Orbits are confined to the region of space for which all three terms are non-negative. In general, this condition will not be satisfied for all points, so that the Abel components have finite extent. From the requirement that at least the origin should be included, we find the following limits on and
The factor in the integrand as well as the upper limit of the integral are different for each of the three Abel component types NR, LR and SR, and are discussed in §§ 2.3.2–2.3.4 below. The lower limit of the integral has to be at least as large as the smallest value possible for the variable . This limiting value depends on the choice of the DF parameters and in (2.9), as is shown in Fig. 1 (cf. Fig. 7 of DL91). The boundaries follow from (2.12) and the separatrices and are given by
At a given position , orbits with different values of the integrals of motion , and , and hence different values of , can contribute to the integral (2.10). The restriction to bound orbits () together with the requirement that , and all three have to be non-negative determines the part of the integral space that is accessible by orbits that go through . An example of the resulting tetrahedron in the -space is shown in Fig. 2 (cf. Fig. 1 of MD99). The largest possible value of is given by the top of this tetrahedron
which is thus a function of the position . At the origin , which is the central value of the potential . In what follows, we normalise by setting , so that .
2.3.2 Non-rotating components (NR)
Since the non-rotating component type can exist everywhere in the accessible integral space (the tetrahedron in Fig. 2), we simply have that . Spatially the NR components are thus bounded by the surface .
The factor follows from the cross section of the -plane within the tetrahedron and can be written in compact form as
where is the beta function of three variables333The beta function of variables in terms of the complete gamma function is defined as .. Since is independent of it can be taken out of the integral (cf. eq. [3.10] of DL91), which then becomes of Abel form. Unfortunately, the inversion of eq. (2.10) for any chosen moment , including the case , is generally impossible, as the left-hand side is a function of three variables, while the DF depends on only one variable, . The density specified along any given curve will define a different . A case of particular interest is to choose the density along the short axis to be . This defines a unique , and hence gives everywhere. Kuzmin’s formula applied to similarly defines the density everywhere. For single Abel DF components these will not be the same, except in the spherical limit (see Appendix A.3).
Since the orbits have no net rotation, the velocity moments are only non-zero when , and are all three even, and vanish in all other cases.
2.3.3 Long-axis rotating components (LR)
The long-axis rotating component type only exists in the part of the integral space that is accessible by the (inner and outer) long-axis tube orbits. Within the tetrahedron for all orbits this is the region for which at . It follows that .
The term follows from the cross section of the -plane within the tetrahedron and with the above boundary plane at . Without any further constraint this results in zero net rotation, because each orbit with positive rotation around the long axis with , is balanced by an orbit with opposite direction of rotation with . Therefore, we restrict to orbits with , resulting in maximum streaming around the long axis for each LR component. This reduces the accessible integral space, and thus also the term , by a factor of two, so that the latter becomes
with , the parameters and defined as
which for are non-negative, and
The function is defined in Appendix B, where we evaluate it in terms of elementary functions (odd ) and elliptic integrals (even ).
The LR components have maximum streaming around the long axis, but the motion parallel to the intermediate axis and short axis cancels. As a result, the velocity moments vanish when or are odd444Since is even, the factor in eq. (2.16) is always real.. Multiplying with results in maximum streaming in the opposite direction. By choosing different weights for both senses of rotation, we can control the direction and the amount of long-axis streaming motion for each LR component.
2.3.4 Short-axis rotating components (SR)
The short-axis component type reaches the part of integral space accessible by the short-axis tube orbits. Within the tetrahedron for all orbits this is the region for which both at and (Fig. 2). The latter requirement is equivalent to . In this case, .
The form of the term depends on the cross section of the -plane within the tetrahedron and with the above two boundary planes. In case each SR component has maximum streaming around the short axis (), it is given by
The parameters and follow from and defined in (2.17) by interchanging , and in turn and follow from and by interchanging . For the terms we have two possibilities, I and II,
where is given in Appendix B, and and follow from
For the assignment of the labels and , we discriminate between four cases
The SR components have maximum streaming around the short axis, so that the velocity moments vanish when or are odd. Multiplying with results in SR components with maximum streaming around the short axis in the opposite direction.
2.4 Combination of multiple DF components
Until now, we have chosen the Abel DF to be a function of a single variable , and we have separated it in three component types, NR, LR and SR, but we have not made any assumption about the form of the DF (apart from the obvious requirement that it has to be non-negative everywhere and that it decreases to zero at large radii). Following MD99, we choose the DF to be a linear combination of basis functions of the form
which, like the velocity moments (2.10), are non-vanishing as long as . The exponent is a (non-negative) constant.
Once the Stäckel potential (2.3) is known by defining the function , we can use the above relations (§ 2.3) together with the expressions in Appendix B, to compute for a given basis function the velocity moments (2.10) for the NR, SR and LR components in an efficient way, where at most the integral over has to be evaluated numerically. For the NR components this integral can even be evaluated explicitly, resulting in
where (cf. eq. 2.14).
Each DF component and corresponding velocity moments thus depend on the choice of the DF parameters , and , the type of component, and for the rotating components (LR and SR), they also depend on the sense of rotation around the axis of symmetry. By summing a series of DF basis functions over , and , one might even expect to cover a large fraction of all physical DFs. Due to the different values of and , such a sum of DF components is no longer a function of the same, single variable , so that the ellipsoidal hypothesis does not apply. Consequently, it becomes possible to construct (nearly) self-consistent dynamical models, with the (combined) luminous mass density equal (or close) to the mass density associated to the potential.
We describe how the intrinsic velocity moments can be converted to projected velocity moments on the plane of the sky. Alternatively, these line-of-sight velocity moments follow as moments of the LOSVD, which we show can be calculated in a straightforward way for Abel models. Parameterising the LOSVD as a Gauss-Hermite series, we obtain the observable quantities: the surface brightness, the mean line-of-sight velocity , velocity dispersion , and higher-order Gauss-Hermite moments , , …
3.1 From intrinsic to observer’s coordinate system
In order to calculate line-of-sight velocity moments, we introduce a new Cartesian coordinate system , with and in the plane of the sky and the -axis along the line-of-sight. Choosing the -axis in the -plane of the intrinsic coordinate system (cf. de Zeeuw & Franx 1989 and their Fig. 2), the transformation between both coordinate systems is known once two viewing angles, the polar angle and azimuthal angle , are specified. The intrinsic -axis projects onto the -axis, which for an axisymmetric galaxy model aligns with the short axis of the projected surface density . However, for a triaxial galaxy model the -axis in general lies at an angle with respect to the short axis of . This misalignment can be expressed in terms of the viewing angles and and the triaxiality parameter (defined in eq. 2.6) as follows (cf. eq. B9 of Franx 1988)
with and . A rotation over transforms the coordinate system to , with the -axis and -axis aligned with respectively the major and minor axis of , whereas is along the line-of-sight.
The expressions in § 2.3 involve the velocity components in the confocal ellipsoidal coordinate system . The conversion to line-of-sight quantities can be done by four successive matrix transformations. First, we obtain the velocity components in the first octant of the intrinsic Cartesian coordinate system via , of which the first element is given by (cf. DL91)
and the other elements follow horizontally by cyclic permutation of and vertically by cyclic permutation of . The second matrix uses the symmetries of the orbits to compute the appropriate signs of the intrinsic Cartesian velocities in the other octants. The result depends on whether or not the orbit has a definite sense of rotation in one of the confocal coordinates. For the three types of Abel components this results in the following matrices
Finally, the conversion from the intrinsic to the observer’s Cartesian velocities involves the same projection and rotation as for the coordinates. We represent these two coordinate transformations respectively by the projection matrix
and the rotation matrix
In this way, we arrive at the following relation
where the full transformation matrix is thus a function of , the constants and the viewing angles .
3.2 Line-of-sight velocity moments
We can now write each velocity moment in the observer’s Cartesian coordinate system as a linear combination of the velocity moments in the confocal ellipsoidal coordinate system
with . The coefficients are products of elements of the transformation matrix in eq. (3.6). They can be obtained with the following recursive algorithm
with the first order expressions given by
and . The index is the th element of the vector , where the number of integers 3 () is equal to the value of the velocity moment index , and similarly and .
The line-of-sight velocity moments now follow from (numerical) integration of along the line-of-sight
which are thus functions of position on the sky plane.
3.3 Line-of-sight velocity distribution
where we have introduced the LOSVD
Although the integral over in general can only be evaluated numerically, we show that for the choice (2.9) of the DF, the double integral over the velocities can be simplified significantly.
Our analysis generalises the results for the well-known spherical Osipkov–Merritt models. We describe the spherical limit together with the elliptic disc and large distance limit in Appendix A, while we present axisymmetric Abel models in § 6.1.
3.3.1 Abel LOSVD
Substituting the expressions (2.2) for the integrals of motion in , we obtain
and similarly and by cyclic permutation of , we can write the expression for as
For a given position , each value of thus defines the surface of the unit sphere in the variables . In these variables, we can write the integral of the DF over velocities, i.e., the stellar mass density, as
This is the same expression as for the zeroth-order velocity moment of the DF, , in eq. (2.10), where is equal to the integral between square brackets.
The matrix in eq. (3.6) provides the conversion from the velocity components in the confocal ellipsoidal coordinate system, , to those in the observer’s Cartesian coordinate system, . Hence, for a given line-of-sight velocity , we find
The coefficients , and are defined as
and normalised with respect to given by
These coefficients are functions of the position , the constants and the viewing angles through the components of the matrix , and also depend on the DF parameters and through the terms , and . It follows that
which is a function of the variable .
We thus find that each combination of values of and results in the cross section of the surface of the unit sphere in eq. (3.15) with the plane in eq. (3.17), i.e., a circle, in the variables . We rotate the latter coordinate system such that the normal vector of the plane of the circle coincides with the -axis of the system given by
where the rotation angles and follow from
In these coordinates the circle is conveniently parameterised as
where . We can now rewrite the integral between square brackets in eq. (3.16) as
where the vector 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 and and differs for the different types of Abel components as we show below.
where after changing the order of integration in the last step, the upper limit of is given by , with
and vanishes when exceeds the ’terminal velocity’ . The expressions for and follow from eqs (3.19) and (3.26), whereas and are different for each of the three Abel component types and are considered next.
3.3.2 Non-rotating components (NR)
As for the intrinsic moments in § 2.3.2, we have for the non-rotating component type that , and, since the full integral space is accessible, , independent of and .
In the case of a basis function as defined in eq. (2.24), the integral over can be evaluated explicitly resulting in
3.3.3 Long-axis rotating components (LR)
The integral space accessible by the (inner and outer) long-axis tube orbits is given by at , so that immediately