Full text of DOI:10.1046/j.1365-8711.2000.03056.x

Mon. Not. R. Astron. Soc. 311, 297±306 (2000)
Approximate third integrals for axisymmetric potentials using local
È
Stackel fits
V. De Bruyne,w F. Leeuwin and H. Dejonghe
Astronomical Observatory, Gent University, Krijgslaan 281 S9, B-9000 Gent, Belgium
Accepted 1999 August 9. Received 1999 August 9; in original form 1999 May 4
A B S T R AC T
È
We use a set of Stackel potentials to obtain a local approximation for an effective third
integral in axisymmetric systems. We present a study on the feasibility and effectiveness of
this approach. We apply it to three trial potentials of various flattenings, corresponding to
nearly ellipsoidal, discy and boxy density isophotes. In all three cases, a good fit to the
È
È
potential requires only a small set of Stackel potentials, and the associated Stackel third
integral provides a very satisfactory, yet analytically simple, approximation to the trial
potentials effective third integral.
Key words: methods: analytical ± methods: numerical ± galaxies: elliptical and lenticular,
cD ± galaxies: kinematics and dynamics ± galaxies: structure.
them. Spectral analysis offers the possibility to simplify the
expression for the orbits and to reduce storage requirements for
1
I NTRO DUC TION
It is well known that axisymmetric potentials do not in general
allow a global and isolating third integral in addition to the energy
and the component of the angular momentum along the symmetry
axis. However, for many axisymmetric potentials, numerical
experience shows that the majority of orbits are constrained by an
effective third integral over a large number of orbital periods. As
appears from modelling observed galaxies, a third integral is often
required to describe their dynamical structure (e.g. Binney, Davies
& Illingworth 1990; Dejonghe et al. 1996). Hence it seems
unescapable that distribution functions depending a priori on three
integrals should be available for general models of elliptical
galaxies.
the orbits (Binney & Spergel 1982; Papaphilippou & Laskar 1996,
1998; Carpintero & Aguilar 1998; Valluri & Merritt 1998).
However, for this type of modelling methods, the computational
cost is relatively high. Moreover, smoothing is necessary because
the singular boundaries of the orbital densities produce a fairly
awkward sum of numerical functions, and the distribution function
is a numerical function out of which the physical content may not
be very easily extracted. The implicit reference to a third integral,
which makes Schwarzschild methods so flexible, is at the same
time a handicap for extracting information about the role played
by this integral.
Analytical techniques would be more explicit. However, no
expression for a global third integral is known in general
In a stellar dynamical context, the problem of approximating
the effective third integral is logically connected to the philosophy
behind the construction method of dynamical models for galaxies:
obviously, numerical models can do with numerical, or even
implicit third integrals, while analytical models need analytical
approximations.
È
axisymmetric potentials ± except for Stackel potentials (cf.
de Zeeuw 1985; Dejonghe & de Zeeuw 1988). Therefore
analytical approximations have been built using perturbation
methods. Systems deviating moderately from spherical symmetry
have been considered, and an approximate integral derived, which
reduces to the total angular momentum in the spherical limit (see
Petrou 1983). Gerhard & Saha (1991) report on the use of three
different perturbation techniques applied to a model that is
initially spherical. (1) Methods based on the KAM theory, which
have a validity that is often restricted to very small perturbations.
In practice, they find indeed that, whatever their order, these
methods fail to track changes in the phase-space topology that
occur when spherical symmetry is broken by a significant amount.
(2) The averaging method (see, e.g., Verhulst 1979 and de Zeeuw
& Merritt 1983) in their example also fails to track the box-orbits
ꢀLz ˆ 0† that emerge in a flattened system, and it gives at first
order only a rough approximation to the third integral. (3) Finally,
a resonant method using Lie transforms does yield a good
Amongst the numerical modelling methods, those based on
Schwarzschild's approach (Schwarzschild 1979) are undoubtedly
the most direct ones. In a given potential, a library of orbits is
computed, while the time-averaged properties of the orbits are
stored. The orbit library is supposed to sample integral space. A
reproduction of the observables is sought by combination of
orbits, populated with a non-negative number of stars (e.g. Rix
et al. 1997; van der Marel et al. 1998). This powerful and very
general technique allows us to construct general three-integral
models, without any explicit reference to the third integral, since
in principle orbits can be labelled using any phase-space point on
w E-mail: vero@izar.rug.ac.be
q 2000 RAS

298
V. De Bruyne, F. Leeuwin and H. Dejonghe
approximation to the third integral. It tracks the new orbital
families emerging around resonances, and may in principle be
carried on to high order. The analytical expressions are, however,
rather intricate, specially for orders higher than 1.
boxy and discy isophotes. These trial potentials are presented in
Section 4. The results can be found in Section 5. Section 6 gives a
discussion on the application of the method, and the conclusions
are given in Section 7.
The last technique has been used at first order by Dehnen &
Gerhard (1993) to construct three integral oblate models for the
perturbed isochrone sphere. An application to the boxy E3±E4
galaxy NGC 1600 has recently been presented by Matthias &
Gerhard (1999).
È
T H E D E S I G N O F T H E S E T O F S TAC KE L
P OT E N T I AL S
2
2.1 The principle
However, one may need to model a galaxy with a potential that
is far from any known integrable potential, so that global
perturbation techniques become ineffective. An alternative
approach is then to derive a so-called `partial integral'. The idea
is that a third integral that can be easily determined for certain
groups of orbits is probably also a good approximation for similar
orbits. de Zeeuw, Evans & Schwarzschild (1996) and Evans,
Orbits with a given set of integrals ꢀE; J†; with E the (positive)
binding energy and J the component of the angular momentum
along the rotation axis, fill a certain volume S in space (a
meridional section of it is shown in Fig. 1b). The intersection of
this volume with a meridional plane is called the zero-velocity
curve (ZVC). One proves easily that orbits with a set of integrals
(E0, J), with E0 . E, fill a volume that is completely embedded in
the volume filled by orbits with the set of integrals (E, J).
Similarly, orbits with E and J0 larger than J, have a ZVC that lies
inside the ZVC defined by (E, J). This means that all orbits in the
shaded rectangle labelled R in integral space (see Fig. 1a) will
remain interior to the shaded region indicated as S in Fig 1(b).
Any distribution function defined over R in integral space can
È
Hafner & de Zeeuw (1997) derived such a partial integral for thin
and near-thin tube orbits in scale-free potentials.
Yet another well-known approach is the use of separable
potentials. One is led to this idea because motion around the origin
can be expanded in a Taylor series and, since the lowest order,
non-zero terms are quadratic, can be described as perturbed
harmonic motion. Moreover, van de Hulst (1962) has shown that
motion around the origin can be seen as (unperturbed) motion in a
È
Stackel potential. Basically, one then uses a much enlarged set of
reference integrable potentials instead of merely quadratic ones.
È
The Stackel potential is found by fitting term-to-term the Taylor
series for the original potential (up to quartic terms) to the
È
expansion for the Stackel potential. Local fits of similar nature
have been performed by de Zeeuw & Lynden-Bell (1985) and
Kent & de Zeeuw (1991), who applied it to the solar neighbour-
hood, by expanding the effective potential in the vicinity of
circular orbits.
However, such a local fitting does not work for all orbits. For
instance, it may fail for orbits with large radial extent. In many
cases, though, an orbit can still be approximated by motion in a
Figure 1. Panel (a): A rectangle R in integral space, of which the `upper
left corner' corresponds to a gridpoint (E, J). Panel (b): All orbits with
È
Stackel potential, provided that the potential is chosen to be a
E0 > E and J0 > J will remain interior to the region S. In practice, the
good average approximation to the true potential in the region
covered by that particular orbit (see Kent & de Zeeuw 1991). A
report on the calculation of such a good average and global
approximation to the Galactic potential, together with some
indications for the numerical implementation, can be found in
Dejonghe & de Zeeuw (1988) and Batsleer & Dejonghe (1994).
In this paper we want to improve the generality of the
application of this family of potentials. We propose to obtain a set
0
È
local Stackel potential is determined in the rectangle S .
È
of local approximations to the potential using Stackel potentials,
each of which is an average approximation to the original
potential in some region. Rather than expanding around an
equilibrium point, which would make the extent of the region
where the approximation holds difficult to handle, we partition
integral space and fit the potential in the corresponding regions of
configuration space. We use the QP method (Dejonghe 1989) to
È
adjust each local Stackel potential. Section 2 describes this. Once
È
the set of local Stackel potentials is obtained, we proceed by
checking the quality of the orbits representation, and of the
approximation for an effective third integral, where it exists.
Obviously, we do not expect to describe the eventual stochastic
motion this way. Moreover, we may not be able to reproduce all
minor orbital families. These questions are also addressed in
Section 2. In Section 3 the actual fitting procedure is explained. Its
performance is illustrated with three trial potentials, which
correspond to mass densities with respectively nearly ellipsoidal,
Figure 2. For one value of
corresponding to seven different values for J and points used to determine
seven values for I3.
E in the grid, zero-velocity curves
q 2000 RAS, MNRAS 311, 297±306

Approximate third integrals for axisymmetric potentials
299
therefore be associated with a potential that needs only be defined
in S. Hence a subdivision of the (E, J)-plane into a number of
rectangles R is equivalent to subdividing space into a number of
(i) the proportion of non-regular orbits. If there is no additional
integral of motion besides E and J, the orbits should fill uniformly
the area bounded by the zero-velocity curve on their SoS
(Richstone 1982). If there is a third effective integral (effectively
meaning over the time the orbits are integrated), orbits are regular,
and appear as curves on the SoS. Experience has taught that in
axisymmetric models, stochastic orbits exist only for very small
values of J. In models with a core, as considered here, these
orbits are due to the large excursions in R and z for small J,
which allow resonances between the two degrees of freedom
(Merritt 1999).
È
bounded domains S. For each of these domains, a Stackel
potential can be determined that is locally a good approximation
to the original potential. If, in the end, the whole (E, J)-plane is
covered by a set of rectangles R(E, J), also space will be covered by
a set of domains S(E, J), and the original potential will be
È
completely fitted by a set of locally fitted Stackel potentials.
(ii) the minor families of orbits which may be present around
È
2.2 Practically
some resonances in the original potential that our Stackel potential
In practice, the specification of a set of rectangles R in integral
space is equivalent to the specification of a grid, each gridpoint of
which being the `upper left corner' of R. A grid in E can be
constructed by dividing the system into equal-mass shells and
computing the circular orbit energy for each shell radius. For
every value of E, we consider an equidistant grid in J, with seven
values in the interval ]0, Jmax(E)[, the upper bound corresponding
to the circular orbit with that energy (see also Fig. 2).
does not generate. We call them unrecovered minor families. The
thin tube orbit at fixed (E, J) determines the main orbital family.
This family encircles the origin on the (z, vz) SoS, or a point
(Rth, 0) on the (R, vR) SoS, where Rth(E, J) is the point where the
thin-tube orbit of given (E, J) intersects the equatorial plane.
Minor orbital families may be generated by other resonances than
the ones for small J, and will show up on a SoS as `islands'
surrounding other points.
È
In Stackel potentials, orbits are determined by three integrals of
È
(2) Orbital densities: Since the Stackel potential is aimed to be
motion (E, J, I3). Following van der Marel et al. (1998), the
different third integrals are parametrized with angles vi, that
correspond to the angle between the horizontal axis and the radius
to the most external of the two points where the orbit intersects
the ZVC (see Fig. 2). The point of the ZVC where I3 ˆ I3max is
found. This is the point where the thin tube orbit touches the
ZVC. The coordinates (R, z) of that point determine the angle
vmax ˆ p 2 arctan{z=‰R 2 RcꢀE†Š}. The vi are chosen linearly
between 0 and vmax, from which the values for I3(vi) are derived.
used for modelling, which is essentially assembling orbits to fit a
given density, a good reproduction of the orbital densities is
important.
These orbital densities n(R, z; E, J, I3) are functions of (R, z) for
each given orbit, defined for any set (E, J, I3). The orbital densities
for both potentials can be compared by measuring the fraction of
the mass that is located in the same place in both potentials. If the
total mass of every orbit is normalized to 1, and the cells surface is
constant over the grid, the mass fraction correctly located (MC)
can be calculated as
È
As for the fit itself, it suffices to perform the Stackel fit to the
potential associated with a rectangle R(E, J) only in the region
interior to the velocity curve S(E, J). In practice, we fit within the
X
MC ˆ 1 2 dM ˆ 1 2
j‰norꢀRk; z`† 2 nSꢀRk; z`†Šj=2:
ꢀ1†
smallest rectangle S0
which encompasses S(E, J)
.
(E, J)
k;`
The construction of a set of potentials that covers the complete
system potential is a multistep process. It is most efficient to start
In this, dM is the mass fraction that is not located in the same cell
in both potentials, nor stands for the original density, and nS stands
È
for the density associated with the Stackel potential. Double
counting is avoided by dividing the sum by 2.
È
with a global fit to the potential, i.e., fitting a Stackel potential in
the most extended domain S0, determined, e.g., by the data. This
domain defines the smallest value for E. The following steps
improve the fit to the potential by performing a number of local
fits in smaller domains.
To calculate this, each orbit determined by a given set (E, J, I3),
is integrated over a large number of periods, and its orbital density
n(R, z; E, J, I3) (normalized to 1) is computed. This is done on a
rectangular grid covering ‰RminꢀE; J†; RmaxꢀE; J†Š  ‰0; zmaxꢀE; J†Š,
by evaluating the time fraction spent in each cell.
The choice of domains where a fit will be performed is in
È
principle arbitrary, and depends only on how well the Stackel
potential has to approximate the system potential (within the
limitations of what a regular potential can achieve).
In principle, the orbit should be integrated until the orbital
density reaches quasi-stationarity. This we may define by
demanding that during a given time interval the value nkl in
each cell (Rk, zl) has varied by less than some small fraction. In
practice, it takes an extremely long time to reach stationarity for
orbits close to a resonance m : n with m, n large, or having a
moderate degree of stochasticity ± these two cases being difficult
to distinguish (cf. Binney & Tremaine 1987, p. 176).
2.3 When is a fit considered to be a good fit?
The fit is checked by comparing orbits in both potentials. Orbits
start from the ZVC defined by (E, J), at the point determined by
the value of I3. The integration of orbits uses a fourth-order
Runge±Kutta scheme with variable time-steps, proportional to the
smallest of the radial and azimuthal periods, estimated respec-
tively as TR , R=VR and TV , 2pR2=J. This ensures conserva-
Also modelling based on the Schwarzschild method has to deal
with stationarity (see, e.g., Schwarzschild 1993 and Wozniak &
Pfenniger 1997). Since modelling aims at constructing equili-
brium models, only orbits for which the orbital density reaches
stationarity in limited time (i.e., less than a Hubble time) need in
principle to be taken into account. Accordingly, the non-stationary
orbital densities would not need to be precisely reproduced in the
tion of energy over 100 TV with a precision better than 1026
.
One can consider various criteria for the comparison.
È
(1) Surfaces of section: The Stackel potential should generate
orbits that are very similar to those in the original potential. A first
inspection can be done by comparing surfaces of section (SoSs) in
both potentials. What we mainly expect to find this way is:
È
Stackel potential.
(3) Conservation of I3: If locally a Stackel potential is a good
È
q 2000 RAS, MNRAS 311, 297±306

300
V. De Bruyne, F. Leeuwin and H. Dejonghe
È
approximation, the Stackel third integral I3 should be approxi-
mately constant along the regular orbits. Since the goal is to
provide a good approximation for an effective third integral for the
largest possible number of orbits, the variation on I3 should be
made as small as possible. If possible, fits on smaller domains will
be performed if this improves the approximation for the integral.
The approximation is evaluated through the maximal variation of
I3 along the orbit with respect to its initial value, dI3/I3.
The fitting procedure is based on a Quadratic Programming
method. The aim is to find a linear combination of basis functions
that provides a good approximation to the original function. When
fitting the potential, the basis functions are
GM
Fꢀt† ˆ 2
:
ꢀ2†
s
ꢀd ‡ t p†
The basis functions and their coefficients are chosen by
minimizing the variable
On the other hand, one can argue that what really matters is a
nominal variation of I3, dI3/I3max(E, J): for a given (E, J), how
precise can our identification of orbits be when we estimate the
third integral of motion by I3 within [0, I3max(E, J)]? A related
question is how precise this identification should be or,
alternatively, how large is the tolerated variation on I3?
X
x2 ˆ
wk`‰foriginalꢀrk; z`† 2 fcalculatedꢀrk; z`†Š2;
ꢀ3†
k;`
with k and ` an index covering the grid points in the domain S0.
The weights wk` can be used to give different relative weights to
the points in the grid, and f stands for the function to fit.
È
The Stackel potentials will be used as basis for a dynamical
model. One of the aims of these models is to reveal the role of the
third integral. Hence it is important to know that all structure
caused by the dependence on this third integral can be traced.
In an ellipsoidal coordinate system 2g < n < 2a < l, with a
and g negative, and the focal distance D ˆ ja 2 gj (de Zeeuw
p
È
1985), a Stackel potential can be written as
È
Consequently, the Stackel potential should be precise enough, so
Vꢀl; n† ˆ glFꢀl† ‡ gnFꢀn†;
with F an arbitrary function [here corresponding to (2)] and
ꢀ4†
that the uncertainty on the resulting approximation for I3 is not
larger than the scale of the smallest feature in I3 of the distribution
function. For example, in the limit where the dynamical model
would turn out to be essentially 2-integral, there would be no need
to set requirements on the constancy of I3.
l ‡ a
gl ˆ
:
ꢀ5†
l 2 n
This makes determining the maximum tolerance for errors in
the fit to the potential an iterative process, which includes the
modelling of the galaxy. First, a fit to the potential has to be
obtained, and a model has to be constructed with that potential.
The results of the model will then indicate whether the error on I3
was small enough or not.
We now present the relevant formulae from Mathieu &
Dejonghe (1996), adapted for an axisymmetric system, that lead
to an expression for the density in terms of the basic functions.
The expression for the axisymmetric density is
rꢀl; n† ˆ g2lc0ꢀl† ‡ gn2c0ꢀn† ‡ 2glgnc‰n; lŠ;
where c[n, l] is the first-order divided difference of c:
ꢀ6†
(4) Topology of the orbits: As pointed out by Gerhard & Saha
(1991), approximate conservation of I3 along orbits alone does not
ensure that I3 provides a correct labelling for orbits (see also de
Zeeuw & Merritt 1983). One still needs to check that the topology
of orbital tori is similar to that of constant (E, J, I3) surfaces: for
then each orbit may be uniquely associated with one value of I3.
This is done by comparing surfaces of section for orbits in the
original potential, and I3 ˆ constant curves for given (E, J).
cꢀl† 2 cꢀn†
c‰n; lŠ ˆ
;
ꢀ7†
l 2 n
and c 0(l) is the derivative c[l, l].
The connection between c(t) and F(t) [like in (2)] is given by
t ‡ g
2pGcꢀt† ˆ 2ꢀt ‡ g†F0ꢀt† 2 Fꢀt† ‡ 2
‰Fꢀt† 2 Fꢀ2a†Š: ꢀ8†
t ‡ a
È
HOW TO OBTA IN A STACK EL POTENT IA L
3
The density can be expressed in terms of the basic functions by
means of a third-order divided difference:
The fitting procedure is derived from a deprojection method for
triaxial systems presented in a paper by Mathieu & Dejonghe
(1996). In that paper, a family of potential-density pairs is
presented that can be used as building blocks for triaxial mass
models. The spatial mass density and the potential can be both
expressed in terms of the same basis functions.
rꢀl; n† ˆ H‰l; n; l; nŠ;
with H(t) defined as
Hꢀt† ˆ ꢀt ‡ a†2cꢀt†;
with c(t) given in (8).
ꢀ9†
ꢀ10†
Using this method, it would technically be possible to obtain a
È
Stackel approximation in two ways: (1) fitting the spatial mass
density (2) fitting the potential.
The divided difference of order n 2 1 of a function G(t) is a
function of the divided difference of order n 2 2, and is given by
The first way is not ideally suited for performing local fits,
since the calculation of a local value for the potential requires an
integration of the density over the entire volume of the galaxy. In
the case of a fit in a restricted part of space, there is obviously no
control on the behaviour of the mass density at larger distances,
and this can have serious implications for the potential obtained
through integration. The correspondence between the forces
generated by the original and the fitted potential is likely to be
a more useful indication of the quality of the approximation than
the correspondence between the spatial densities. For a global fit,
G‰t1; t3; ¼; tnŠ 2 G‰t2; t3; ¼; tnŠ
G‰t1; ¼; tnŠ ˆ
:
ꢀ11†
t1 2 t2
4
T H RE E T E S T P OT E N T IA L S
We build a smooth flattened model by combining a few spherical
harmonics Y0l ꢀu† (Binney & Tremaine 1987), as follows:
jbj
rꢀr; u† ˆ C‰r0ꢀr†Y00ꢀu; f† ‡ br2ꢀr†Y20ꢀu; f† ‡
r4ꢀr†Y04ꢀu; f†Š;
5
È
on the other hand, it does make sense to construct a Stackel
potential through a fit on the spatial density.
ꢀ12†
q 2000 RAS, MNRAS 311, 297±306

Approximate third integrals for axisymmetric potentials
301
1
r0ꢀr† ˆ
r2ꢀr† ˆ
;
;
ꢀ13†
2
3
ꢀ1 ‡ r2†
r2
ꢀ14†
ꢀ15†
ꢀ1 ‡ r2†
r4ꢀr† ˆ r2ꢀr†:
b is a parameter which determines the flattening (b , 0 for oblate
p
systems). We set C ˆ 2=ꢀp p† for a total mass equal to 1. This
simple choice for the model ensures r / r24 and a fairly constant
flattening. The resulting density for b ˆ 20:3 is shown in Fig. 3.
It is nearly ellipsoidal, with axis ratio b=a . 0:8. We will in the
following refer to the corresponding potential as Cell.
The corresponding potential is
Figure 3. The spatial density, with its axis ratio, and potential
corresponding to Cell, given by the sum of harmonics given in equations
(12) and (16), with b ˆ 20:3.
cellꢀr; u† ˆ 24pGC‰c0ꢀr†Y00ꢀu; f† ‡ bc2ꢀr†Y20ꢀu; f†
jbj
‡
c4ꢀr†Y04ꢀu; f†Š;
ꢀ16†
5
where the fl are related to the rl terms by equation (2-122) in
Binney & Tremaine (1987).
As explained in Section 2, we use a grid in integral space that
defines the domains for the fit. The grid contains 16 points in E
(noted Ei), seven points in J (noted Ji) and seven angles vi.
For each of the domains, a fit can be produced, based on the x2
of equation (3), if the focal distance D is set. The choice of D of
the spheroidal coordinate system affects the quality of the fit. An
average D is obtained by computing a number of orbits in the
original potential, and evaluating the focal distance for which
sections of the l ˆ constant ellipses most closely follow the inner
and outer envelopes Rmin(z), Rmax(z) of the orbits. Of course, there
is more variation in D for the more extended domains S0.
We started with a global fit for the potential Cell, i.e., a fit in the
largest domain S0 (extending to about 180 in R and in z),
corresponding to the orbits that have (E16,0). We take J ˆ 0
instead of J1 in order to avoid the possibility that orbits with
smaller E might not be covered by the potential because their
Figure 4. The contours for the two force components generated by Cell
0
È
(full lines) and the Stackel potential for the domain S (11,0) (dashed
lines).
J , ꢀJ1†
.
E16
This combination (E16,0) is obviously the most unfavourable for
the fit, and the variations in I3 are expected to set a standard for the
worst case: the maximum (dI3)max for all seven orbits with
different vi is taken as the maximum tolerable variation. We
demand that the variation of I3 along the orbits with lesser spatial
extent (i.e., E0 . E16; J0 . 0 and the 7 vi, with (E0, J0) on the grid
in integral space) is not worse than this. In this case, it turns out
that the variation on I3 is larger for orbits with (E15, J1). Thus the
next fit is done in the domain S0 determined by (E15, 0).
In Fig. 5 the points in the (E, J)-grid for Cell are displayed as an
example of how a set of potentials can be built; different symbols
indicate the different potentials used. As can be seen from the
different symbols, some potentials of the set are used for a small
number of orbits, while others are used for a large number of
orbits. The fact that the potential fitted in the largest domain is
used only for orbits with only one value for E (indicated by the
pluses) suggests that this method can really improve the
approximation. Also the orbits that remain close to the centre
seem to require a potential fitted in a limited part of space.
As a second trial case, a boxy density, with potential Cbox, is
obtained from a combination of harmonics as described in equations
(12) and (16) with b ˆ 20:5, and has an axis ratio b=a . 0:6. The
After this fit, (dI3)max is the error on I3 for the orbits with
(E15, J1), and this new value is now used to decide on the next fit.
In this way, fits are done in the domains corresponding to (E16, 0),
(E15, 0), (E14, 0), (E11, 0), (E2, 0) and (E1, 0).
If the correspondence between the orbits in a fitted and original
potential happens to be very bad, it could be due to a poor
correspondence between the force components of both potentials.
Therefore, prior to the calculation of the orbits, it is useful to
check the behaviour of the force components in the fitted and the
original potential. Fig. 4 shows the contours for the force
contours are shown in Fig. 6. We also used a 16E Â 7J Â 7I grid
3
È
and approximated Cbox with a set of eight Stackel potentials.
Also a Miyamoto±Nagai (MN) potential with intermediate
flattening, b=a , 0:7, is taken as trial potential. MN models
exhibit very strong disciness in the density contours, as can be seen
in Fig. 7. This means that we are considering a difficult case, that is
actually on the verge of being unrealistic for elliptical galaxies.
The grid has the same dimensions as the one described in the
previous section. The MN potential is approximated with eight
È
components generated by Cell (full lines) with the Stackel
potential for the domain S0(E11, 0) (dashed lines) in overlay.
Since the force components in both potentials do not appear too
different, it is sensible to proceed with the comparison of
integrated orbits.
È
potentials, and the set of Stackel potentials is constructed
following the same strategy.
q 2000 RAS, MNRAS 311, 297±306

302
V. De Bruyne, F. Leeuwin and H. Dejonghe
Figure 5. The (E, J)-grid for Cell, with different symbols for different
È
potentials. We used a set of six Stackel potentials to approximate the
Figure 8. Surfaces of section for typical orbits in the MN potential (small
È
dots) and the Stackel potential (large dots).
original potential.
(z, vz) for vR ˆ 0, and for those orbits having vR . 0 just after
crossing.
The SoSs immediately tell us that for the three potentials, all
orbits of our phase-space grid have an effective third integral.
Stochasticity, if present, was mild and would have required very
long integrations to be visible on the SoSs. We did not carry out
such integrations, because the SoSs were only meant as a
preliminary check. In practice, anyway, mildly stochastic orbits
behave very much like regular orbits, and we can hope to
È
approximate them by an orbit in a Stackel potential.
For cell and cbox, all orbits in our library have similar SoSs in
both potentials. There is no evidence for stochasticity, nor for
unrecovered minor families. For the MN potential, all orbits also
appear as regular, but we found a few unrecovered minor orbital
families.
Figure 6. The spatial density, its axial ratio and the potential
corresponding to CBox, given by the sum of harmonics given in equations
(12) and (16), with b ˆ 20:5.
Most of the orbits trapped around resonant tube orbits ± what
we called `minor orbital families' in Section 2.3 ± were present in
È
the local Stackel potential.
in Fig. 8 SoSs of a few orbits are displayed for the MN potential
È
(small dots) and the Stackel potential (large dots). Both panels
show a `minor resonance' evidenced by small `islands'. In the left-
hand panel, showing orbits with large I3, the orbits trapped around
È
minor resonances are well reproduced within the Stackel potential.
The right-hand panel displays orbits with smaller I3. Two orbits
are trapped around a 1:1 resonance which does not exist in the
È
Stackel potential. They are examples of unrecovered minor orbital
families. They are always found for small values of I3. Therefore
they correspond to orbits remaining close to the equatorial plane,
where the MN potential is very much distorted by a strong
disciness.
Figure 7. The spatial density, its axial ratio and the potential for a
Miyamoto±Nagai potential with b=a , 0:7. The contours for the spatial
density are strongly discy.
Even for such extreme `disciness', these orbits only represent
,1 per cent of the orbit library for MN. These obviously are orbits
that we will not be able to take into account when building a
È
model that uses Stackel potentials. Although the unrecovered
5
P RESE NTAT ION OF T HE RESU LTS
minor families are only a small fraction of the entire orbital
library, we will exclude them from the statistics of the other
checks discussed hereafter. As such, we are not taking those
orbital families into account that no one expects to be
5.1 Surfaces of section
Typical surfaces of section display (R, vR) at z ˆ 0, for the orbits
having vz . 0 when they cross the equatorial plane or display
È
approximated by a Stackel potential.
q 2000 RAS, MNRAS 311, 297±306

Approximate third integrals for axisymmetric potentials
303
È
were not fitted by our Stackel potentials. Excluding those, the
average MC is 99 per cent for Cell, 98 per cent for Cbox, and 95 per
cent for MN.
5.2 Orbital weights
For every orbit, the conserved mass fraction MC is calculated
following equation (1), using a 20 Â 20 grid in (R, z). We integrate
each orbit for at least 50TV, and at most 200TV. We check every
50TV whether quasi-stationarity has been reached. Quasi-statio-
narity is assumed when for every cell the orbital density has varied
by less than 1 per cent in 50TV.
5.3 How constant is the StaÈckel I3?
For every potential, the conservation of I3 along the orbits
is estimated by computing the relative variation dr ;
dI3ꢀE; J; I3†=I3ꢀE; J; I3†: At given (E, J), the maximum of this
The results are presented in Fig. 9, which gives the cumulative
distribution for MC, the fraction of mass `correctly located' in the
quantity is derived among orbits with different I3: dmr ax
maxiˆ1;7‰dI3ꢀE; J; I3;i†=I3ꢀE; J; I3;i†Š:
;
orbital densities. A few orbital densities for MN have
MC , 60 per cent; these correspond to minor resonances that
Fig. 10 shows how logꢀdrmax† is affected as more potentials are
added to the set approximating cell
.
The black dots represent the E ꢀJ ˆ 0† of the domain S0 where
the potentials are fitted. Dark grey shades represent small
variations on I3, and light grey shades represent larger variations.
For each new local potential that is considered, a number of
previously void cells of the (E, J) grid are filled and a number of
values are replaced by new ones. It is clear that adding new
potentials to the set improves the conservation of the third integral
along the orbits, and the difference between the global fit (upper
left panel of Fig. 10) and the results of the complete set (lower
right panel) is remarkable.
In the left-hand panels of Fig. 11 we present histograms of dr
for from top to bottom: cell, cbox and MN. For the MN potential,
È
the resonances absent from the Stackel approximation are not
considered, a complete orbital library contains 784 orbits.
The average for dr is 2 per cent for cell, while 75 per cent of the
orbits have dr , 1:7 per cent. For cbox, the average for dr is
,2:6 per cent and 75 per cent of the orbits have dr , 2:7 per cent.
For MN, the average value is ,10 per cent, and 75 per cent of the
orbits have dr , 12 per cent. The variation of dr along each orbit
is of the order of the fitting error on the derivatives of the
potential.
Figure 9. Cumulative distribution for the fraction of mass `correctly
located' in the orbital densities. The full line is for Cell, the dotted line is
for Cbox, and the dashed line is for MN.
r
Figure 10. The evolution in log(d ) as more potentials are added to the set of Stackel potentials approximation Cell. The black dots indicate the E of the
È
domain where the potentials were fitted (J ˆ 0 for all six potentials). Dark (light) grey shades represent small (large) variations on I3.
q 2000 RAS, MNRAS 311, 297±306

304
V. De Bruyne, F. Leeuwin and H. Dejonghe
Figure 11. Histograms of dr in the left-hand panels, of dn [i.e., dI3
normalized to I3max(E, J)] in the middle panels, and of dI3 normalized to
I3max(E) in the right-hand panels.
In the middle panels of Fig. 11 the results in terms of a nominal
variation dn ; ‰dI3ꢀE; J; I3†=I3maxꢀE; J†Š are shown. For the nearly
ellipsoidal potential Cell, 75 per cent of the orbits have
dn , 0:3 per cent, with an average value of ,0:2 per cent. For
C
box, the average is ,0:7 per cent, while 75 per cent of the orbits
have dn , 1 per cent. For MN, the average value is ,1:7 per cent,
and 75 per cent of the orbits have a nominal error on I3 less than
2.3 per cent.
Dehnen & Gerhard (1993) used a similar indicator, but they
normalized dI3 to I3max(E), which gives smaller values for the
nominal variation, as can be seen in the right-hand panels of Fig. 11.
The averages for Cell, Cbox and MN are 0.1, 0.5 and 1.3 per cent
respectively.
5.4 Topology of constant I3 varieties
Ç
We plot on the same graph (Fig. 12): (i) the (R, R) SoS of orbits as
they cross the z ˆ 0 plane (with z_ . 0) in the original potential;
and (ii) the invariant curves defined by I3 ˆ constant (at z ˆ 0).
This is shown here for some intermediate value of E (E4 of our
grid) and J ˆ 0. Radial orbits are in principle difficult to map,
because of the transition between loop and box orbits (see Gerhard
& Saha 1991). They are, however, very well described here, both
for Cell and for Cbox. The agreement is somewhat poorer in the
case of MN for the transition region between box and loop
orbits; however, the topology is still correctly described, in the
sense that we can establish a one-to-one correspondence between
the two sets of curves. This shows that I3 may be used to label the
orbits.
Figure 12. Comparison of SoSs and I3 level curves at given ꢀE4; J ˆ 0†,
for Cell (top panel), Cbox (middle panel), and MN (bottom panel). Dots:
Ç
surfaces of section (R, R) at z ˆ 0 with z_ . 0. Lines: the intersection of
I3 ˆ constant surfaces with the z ˆ 0 plane.
6
DI S C U S S I O N
È
process. In a Stackel potential there is a function of one variable
È
that can be freely chosen. This freedom is advantageous, and can,
of course, be exploited at the fullest when performing a fit to the
galaxy potential.
For the construction of dynamical models, the use of Stackel
potentials yields a number of advantages, amongst them (1) the
property that the density and dynamical moments can be
calculated analytically, and (2) the absence of regularisation
problems. Finally, once the dynamical model has been completed,
the distribution function is known in an easy-to-use and analytical
form. Although these models require a special form of the
potential, they offer a great flexibility during the further modelling
È
Analogous to local Stackel potentials, dynamical models built
within this approximation will yield local distribution functions,
simply because we have adapted the coordinate systems in each
È
domain S where a local Stackel potential is constructed. The use
of a QP-method (Dejonghe 1989) for dynamical modelling gives a
q 2000 RAS, MNRAS 311, 297±306

Approximate third integrals for axisymmetric potentials
305
large freedom for the basis functions. In this case, we can use
basis functions for the distribution function that contribute only in
limited parts of integral space. In practice, these limited parts will
correspond to the rectangles R from the local potentials.
7
C ONC LUSION S
In this paper we present an extension of the available approxima-
tions for the effective third integral in axisymmetric systems,
È
obtained by using a set of Stackel potentials as representation for a
È
The fact that the space of Stackel potentials has measure zero in
È
galaxy potential, instead of one single Stackel potential. We have
studied the feasibility and effectiveness of this method.
the space of all axisymmetric potentials, raises the obvious
È
concern whether Stackel potentials are suitable for the representa-
tion of any galaxy potential.
È
The creation of a set of local Stackel potentials is done through
The scientific objections cluster around the following arguments:
a fit on the system potential. There is a large freedom in the choice
of the domains where the local fits are done and the composition
È
(1) `What about central cusps?' Indeed, central density cusps
of the set of Stackel potentials. The set of approximating local
È
generally cannot be generated by a Stackel potential in an
ellipsoidal coordinate system. However, Sridhar & Touma (1997)
potentials can be extended until the desired precision is obtained.
We have tested the method on three potentials: (A) a harmonical
potential (Cell) that behaves very smoothly with a nearly
ellipsoidal density; (B) a harmonical potential (Cbox) with a
boxy density; and (C) a Miyamoto±Nagai (MN) model with a
density that has an exaggerated discy structure.
È
showed that a Stackel potential in parabolic coordinates can create
a central cusp. Whether these potentials can be used in the scheme
we propose here remains to be investigated. On the other hand, it
may well be possible that there is no point in providing a good
approximation to the third integral in regions where the cuspiness
of the density becomes important. Indeed, it is likely that the
central density cusps are related to the presence of a central
massive object (see, e.g., the reviews by Kormendy & Richstone
1995 and Richstone et al. 1998). This may cause a higher fraction
of ergodic orbits over a Hubble time. Such orbits will fill regions
of phase-space where the distribution function depends essentially
on the energy, or, for less ergodic regions, also on angular
momentum (cf. Merritt 1999). The fraction of regular to ergodic
motion in a real galaxy stands as a major unknown. The most
popular scenario at the moment is that ergodicity has gradually
been introduced in the system, as orbits that pass close to the
central mass concentration were perturbed by it (Gerhard &
Binney 1985; Valluri & Merritt 1998). One may then suppose that
the original orbital tori which have not been disrupted still
underlie most of the features in phase-space, with essentially
homogeneous 6D density between them.
As one may expect, the model with smallest disciness/boxiness
È
is easiest to approximate using sets of Stackel potentials.
However, what really motivates this study is to check the estimate
made for the effective third integral, if we approximate it by the
È
third integral I3 that the Stackel potentials define.
The quality of the approximation is checked in several ways: (1)
surfaces of section, that reveal possible resonances and irregular
orbits, (2) conservation of orbital weights, which is important for
the reconstruction of the spatial density, (3) conservation of I3
along the orbits, in order to validate the labelling of orbits, and (4)
topology of the orbital space. The conservation of orbital weights
and the conservation of I3 are the criteria that can be best
expressed numerically.
According to the expectations, the approximation of Cell has
orbits that are very similar to the orbits in the original potential, as
can be judged from the surfaces of section and the topology. The
results for the conservation of orbital weights (an average of
99 per cent for Cell) and I3 ꢀdr ; dI3ꢀE; J; I3†=I3ꢀE; J; I3† , 2 per
cent and dn ; ‰dI3ꢀE; J; I3†=I3maxꢀE; J†Š , 0:2 per cent for Cell)
confirm the quality of the approximation.
È
(2) `Does a Stackel approximation allow an accurate fit to a
typical galactic potential?' As has been already noticed for some
time, the main orbit families found by numerical integration in
È
general triaxial potentials are present in a Stackel potential
The method also proves to be successful for Cbox, with an
average of 98 per cent for the conservation of orbital weights,
dr , 2:6 per cent and dn , 0:7 per cent. The SoSs and the
topology of the orbits also confirm this good result.
(Schwarzschild 1979; de Zeeuw 1985), but obviously there is no
place in an integrable potential for smaller orbital families nor
stochastic orbits. These minor orbital families appear to occupy
only a small volume fraction of phase-space, as long as figure
rotation is unimportant (Gerhard 1985; Binney 1987) In practice,
The success of the method on the MN potential is somewhat
less, because of the strong disciness. Still, using a small set of
È
È
Stackel potentials, we are able to reproduce most orbits with
satisfactory accuracy, except for a few resonances. The resonances
Stackel potentials do turn out to provide reasonably good global
fits for systems without central mass concentration (Dejonghe
et al. 1996) or for regions beyond the influence of the central mass
concentration (Emsellem, Dejonghe & Bacon 1999). It is true that
È
that are not reproduced by the set of Stackel potentials all have
small values of I3, i.e., the orbits lie in the region where the
disciness is important. Leaving these resonances out of considera-
tion, the orbital weights seem to be well conserved (on average
95 per cent for MN), and the topology of the orbits is well
È
the use of one single Stackel potential assumes a single confocal
system where the streamlines of the mean stellar motion coincide
with the coordinate lines of the ellipsoidal coordinate system.
È
È
Hence the construction of such a single Stackel potential inevit-
reproduced. Also for these potentials, the Stackel I3 can be used as
ably involves some sort of averaging (de Zeeuw & Lynden-Bell
1985; Dejonghe & de Zeeuw 1988), and is not always considered
to be a sufficiently general approach (at least in triaxial cases; see
Binney 1987 and Merritt & Fridman 1996). However, the use of
confocal streamlines seems to be a valid assumption for individual
orbits (Anderson & Statler 1998). Moreover, studying logarithmic
potentials, those authors find that the mean velocities can be well
fitted using local coordinate systems with confocal coordinate
lines, with the focal distances taken within a fairly narrow range.
Our approach allows us to carry this idea one step further, on the
level of the distribution function.
labels for the orbits (dr , 10 per cent and dn , 1:7 per cent). The
potentials of observed elliptical galaxies are generally rounder
than the strongly discy MN models.
Given the positive results found for the three rather different
models considered, it seems to be possible, using a reasonable
number of local StaÈckel potentials, to provide good approxima-
tions for a third integral, suitable for labelling orbits in dynamical
models.
With respect to existing approximations for a third integral in
axisymmetric systems, the advantage of this new one is mainly
that, while it can be envisaged for application to systems with
q 2000 RAS, MNRAS 311, 297±306

306
V. De Bruyne, F. Leeuwin and H. Dejonghe
Gerhard O., 1985, A&A, 151, 279
arbitrary flattening, it also yields a simple analytic expression for
the approximate third integral locally. This method will be fully
exploited if it is used to build, for roughly axisymmetric galaxies,
explicit distribution functions that depend on three integrals.
Gerhard O., Binney J., 1985, MNRAS, 216, 467
Gerhard O. E., Saha P., 1991, MNRAS, 251, 449
Kent S. M., de Zeeuw P. T., 1994, AJ, 102
Kormendy J., Richstone D., 1995, ARA&A, 33, 581
Mathieu A., Dejonghe H., 1996, A&A, 314, 25
Matthias M., Gerhard O., 1999, MNRAS, submitted, astro-ph/9901036
Merritt D., 1999, PASP, 111, 129
RE FE RENC ES
Merritt D., Fridman T., 1996, ApJ, 460, 136
Papaphilippou Y., Laskar J., 1996, A&A, 307, 427
Papaphilippou Y., Laskar J., 1998, A&A, 329, 451
Petrou M., 1983, MNRAS, 202, 1195
Anderson R. F., Statler T. S., 1998, ApJ, 496, 706
Batsleer P., Dejonghe H., 1994, A&A, 287, 43
Binney J. J., 1987, in de Zeeuw T., ed., Structure and Dynamics of
Elliptical Galaxies. Reidel, Dordrecht, p. 229
Richstone D. O., 1982, ApJ, 252, 496
Binney J. J., Spergel D., 1982, ApJ, 252, 308
Richstone D. O. et al., 1998, Nat, 395, 14
Rix H.-W., de Zeeuw P. T., Cretton N., van der Marel R. P., Carollo C. M.,
1997, ApJ, 488, 702
Binney J. J., Tremaine S., 1987, Galactic Dynamics. Princeton Univ. Press,
Princeton
Binney J. J., Davies R. L., Illingworth G. D., 1990, ApJ, 361, 78
Carpintero D. D., Aguilar L. A., 1998, MNRAS, 298, 1
Dehnen W., Gerhard O. E., 1993, MNRAS, 261, 311
Dejonghe H., 1989, ApJ, 343, 113
Schwarszchild M., 1979, ApJ, 232, 236
Schwarszchild M., 1993, ApJ, 409, 563
Sridhar S., Touma J., 1997, MNRAS, 292, 657
Valluri M., Merritt D., 1998, ApJ, 506, 686
van de Hulst H. C., 1962, Bull. Astron. Inst. Neth., 16, 235
van der Marel R. P., Cretton N., de Zeeuw T. P., Rix H.-W., 1998, ApJ,
493, 613
Dejonghe H., de Zeeuw P. T., 1988, ApJ, 329, 720
Dejonghe H., De Bruyne V., Vauterin P., Zeilinger W. W., 1996, A&A,
306, 363
de Zeeuw P. T., 1985, MNRAS, 216, 273
Verhulst F., 1979, Phil. Trans., Royal Society (London), 290, 435
Wozniak H., Pfenniger D., 1997, A&A, 317, 14
de Zeeuw P. T., Lynden-Bell D., 1985, MNRAS, 215, 713
de Zeeuw P. T., Merritt D., 1983, ApJ, 267, 571
de Zeeuw P. T., Evans N. W., Schwarzschild M., 1996, MNRAS, 280, 903
Emsellem E., Dejonghe H., Bacon R., 1999, MNRAS, 303, 495
A
This paper has been typeset from a T X/LT X file prepared by the author.
È
Evans N. W., Hafner R. M., de Zeeuw P. T., 1997, MNRAS, 286, 315
E
E
q 2000 RAS, MNRAS 311, 297±306