A Numerical Model For Triaxial Stellar System In Dynamical Equilibrum Martin Schwarzchild [624275]

197 9ApJ. . .232 . .2363 The Astrophysical Journal, 232:236-247, 1979 August 15
© 1979. The American Astronomical Society. All rights reserved. Printed in U.S.A.
A NUMERICAL MODEL FOR A TRIAXIAL STELLAR SYSTEM IN DYNAMICAL EQUILIBRIUM
Martin Schwarzschild
Princeton University Observatory
Received 1978 November 8; accepted 1979 February 15
ABSTRACT
A sample numerical model has been computed for a triaxial stellar system in dynamical
equilibrium. A four-step procedure was employed. (1) A density distribution with a modified
Hubble profile—to approximate elliptical galaxies—and axis ratios of 1:1.25:2 was chosen.
This figure was further chosen not to rotate. (2) The potential corresponding to the chosen
density distribution was computed. (3) About 1500 orbits were computed with this potential,
one at a time, each covering typically 100 oscillations through the system, i.e., ~ 1 billion years.
These orbits belong to two families, viz., box orbits and tube orbits around the long axis of the
system. (4) A reproduction of the chosen density distribution—in terms of its mass in 285 cells in
each octant—was sought by superposition of a subset of the available orbits, each populated by
an appropriate, nonnegative number of stars. The application of linear programming led to a
numerical solution. The main results are: First, the majority of the orbits computed in step 3
turned out to have three effective integrals, i.e., two nonclassical ones in addition to the energy
integral. Second, the existence of the numerical solution found in step 4 suggests the existence
of triaxial self-consistent systems in dynamical equilibrium with density profiles fitting elliptical
galaxies.
Subject headings: galaxies: internal motions — galaxies: structure — stars: stellar dynamics
I. INTRODUCTION
In the past, elliptical galaxies have generally been
assumed to be well represented by figures of axial
symmetry, such as oblate spheroids. Recently, how-
ever, doubt has been cast on this assumption by new
observational data such as unexpectedly low rotational
velocities (Illingworth 1977) and nonvanishing mean
motions on the projected minor axis (Schechter and
Gunn 1978) as well as by certain results from large-
scale iV-body calculations (Miller and Smith 1979;
Aarseth and Binney 1978). Hence the question arises
whether elliptical galaxies are generally better repre-
sented by triaxial figures, such as ellipsoids (for reviews
see Gott 1977; Binney 1978). If such genuinely three-
dimensional figures are to be equilibrium configura-
tions for long-lived stellar systems like elliptical
galaxies, then the corresponding potential must—
except in special cases—permit three effective integrals,
i.e., two in addition to the classical energy integral
(for the corresponding two-dimensional problem
see, e.g., the early analyses by Contopoulos 1960;
Lynden-Bell 1962; Hénon and Heiles 1964; and the
recent paper by Hunter 1977). Research regarding
the existence of two nonclassical integrals and the
existence of equilibrium configurations without axial
symmetry has become an active field in recent years
in which a variety of analytical methods (e.g., Lynden-
Bell 1963; Freeman 1965; Vandervoort 1978) as well
as A-body calculations have been employed.
The aim of this paper is to describe a numerical
procedure by which three-dimensional self-consistent
configurations can be computed for stellar systems in dynamical equilibrium for which the effects of en-
counters are negligible (§ II). The paper also contains
the application of this procedure to a specific triaxial
model for which a numerical solution for an equilib-
rium configuration was found to exist and has been
computed (§ III).
The procedure here employed, in contrast to Y-body
calculations, cannot give any insight regarding the
formation of a stellar system since it aims directly at
the final equilibrium state. Furthermore, no stability
tests are included except those regarding individual
orbits. On the other hand the present procedure, in
contrast to some of the analytic approaches, should
be applicable to figures deviating grossly from any of
the special classes of figures for which three integrals
are analytically accessible.
II. GENERAL PROCEDURE
The procedure here employed consists of the follow-
ing sequence of five steps.
a) Choice of Density Distribution
The spatial density distribution of the equilibrium
configuration to be constructed is chosen at the outset
and not altered during the procedure. Thus one may
choose a density distribution which, in projection,
simulates an observed galaxy—if the latter is in dyna-
mical equilibrium, i.e., if it is not recently formed or
recently perturbed. The specification of the density
distribution has to include the choice of a density
figure at rest in an inertial frame or of a figure rotating
with a chosen speed around a chosen axis.
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 NUMERICAL MODEL FOR TRIAXIAL SYSTEM 237
For use in step B the density distribution has to be
given in fairly detailed numerical or even analytical
form. For step E the following, generally coarser,
representation will have to suffice for computational –
economy. The space occupied by the model is divided
into N cells and the chosen density distribution is
represented by the mass, D(J) with / = 1 to iV, it
places into each cell.
b) Derivation of Gravitational Potential
Given the density distribution, the potential can be
computed from the Poisson equation by any of the
standard methods. Again for computational economy
in step C, it is advantageous if the potential and its
first three partial derivatives can be expressed at least
partly algebraically, rather than totally numerically,
such as by detailed three-dimensional tables which
require much computer core space and persistent three-
dimensional interpolations.
c) Computation of Orbits
Given the gravitational potential, any stellar orbit,
defined by its starting parameters, can be computed—
one at a time in contrast to A-body calculations—by
any of the standard schemes for numerical inte-
grations. For the present aim it is necessary that all
orbits used in the final configuration be computed for
a time interval long enough, say 100 oscillations
through the system, to represent reasonably securely
its characteristics during the entire equilibrium phase
of the galaxy.
During the computation of an orbit, data charac-
terizing it have to be stored. In future investigations a
variety of data, particularly velocities, may be wanted.
For the narrow aim of the present investigation, how-
ever, only one set of data needs to be stored, namely,
the fraction of time the orbit spends in each of the
cells referred to under step A. Thus, when the calcu-
lation of an orbit, designated by the index /, is com-
pleted a record of the time fractions, £(/, J) with the
cell index / = 1 to A, is available. These B values can
be considered as representing the density distribution
produced by one orbit when averaged over a long time.
d) Existence of Additional Effective Integrals
Before engaging in the laborious final step of this
procedure, one clearly wants to inspect the orbits to
answer the key question whether they have three
effective integrals or only the one energy integral. The
term “effective integral” is here used for any function
in phase space which is substantially constant along
an orbit for a time span as long as the age of the
universe and which is at the same time isolating. This
term therefore includes all exact integrals and formal
integrals but may not be restricted to them. The
interesting question whether orbits with two effective
integrals might exist (Froeschle and Scheidecker 1973;
Contopoulos, Galgani, and Giorgilli 1978) has not
been investigated here. If alb orbits in the potential for a given density
distribution have no additional effective integrals then
each orbit will enter all cells within its limiting poten-
tial surface. Furthermore, all orbits with the same
energy will give identical B(f J) values. Both these
points are easily recognized by inspection of the B(I, J)
tables for the computed orbits. If such inspection
indicates the absence of any additional effective
integral, all B(I9J) tables form a one-parameter
family, the one parameter being the energy. Such a
one-parameter set of orbital density distributions will
generally not suffice to reproduce the model density
distribution chosen at the outset. Hence, if these cir-
cumstances are encountered, it appears reasonable to
conclude that no dynamical equilibrium configuration
exists which has the chosen density distribution.
Accordingly, in such a case to proceed with step E
would seem worthless.
On the other hand, if the majority of orbits in a
given potential have additional effective integrals, then
each one of these orbits will not enter all the cells
within its limiting potential surface but will fill a
volume smaller than that of the energy limit. Specific-
ally, if the B(I9J) tables indicate that orbits of the
same energy provide a two-parameter family of
orbital density distributions, then orbits of all bound
energies will together form a three-parameter family of
orbital density distributions, indicating the existence
of three effective integrals for the majority of orbits.
This level of variety among the orbits is just sufficient
to make it possible, though not certain, that the model
density distribution can be reproduced by an appro-
priate combination of orbital density distributions. In
this case, therefore, proceeding to the final step E of
this procedure is in order.
e) Reproduction of Model Density Distributions
After the computation of, say, M orbits the re-
maining question is: Can the model density distri-
bution, D(J)9 be reproduced by a superposition of the
M orbital density distributions, B(I9 J)9 if each orbit is
randomly occupied by an appropriate number of
stars, here called occupation number C(7) ? Or, more
precisely, does the following set of equations,
M
2 C(I) B(I, J) = D(J) (J = l, N), (1)
7 = 1
have a physically acceptable solution, i.e., one in which
all occupation numbers are nonnegative,
C(/)>0 (/= 1, M) ? (2)
Conditions (1) and (2) represent a problem typical for
linear programming. Fortunately a systematic method
for computing a solution of a problem of this class, or
for proving the nonexistence of a solution, is well
developed by now.
If, for a particular case, linear programming gives
the answer “No solution” for conditions (1) and (2),
two alternative interpretations may be considered. It
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 238 SCHWARZSCHILD Vol. 232
could be that the M orbits included do not cover suffi-
ciently completely the types and parameter ranges of
all possible orbits. In this case additional orbits need
to be computed with the aim of increasing the variety
in the tables, and the linear programming
computation has to be repeated with the inclusion of
additional orbits. On the other hand, if the answer
“no solution” persists even after an exhaustive study
of the full variety of possible orbits, then the alter-
native interpretation would consist of the tentative
conclusion that the chosen density distribution may
not correspond to an equilibrium configuration in
spite of the existence of three effective integrals for the
majority of the orbits. Whether the second alternative
can in reality arise cannot be determined from the
present investigation. (For the particular case of
density distributions with axial symmetry, a case in
which all orbits have at least two integrals and in
which the density distribution varies only in two
dimensions, Lynden-Bell 1962 has disproved the
occurrence of the second alternative just discussed.)
In contrast, in the happy situation that linear pro-
gramming provides a solution of conditions (1) and
(2) for a particular case, one may conclude that the
chosen density distribution does indeed correspond to
a dynamical equilibrium configuration and that this
configuration can be represented by a superposition of
orbits, each occupied by the number of stars given by
the linear programming solution. Clearly, this con-
clusion should be accepted with the caution appro-
priate for any result derived by numerical methods in
which continua are replaced by finite sets.
There remains to be discussed the question of
uniqueness of a solution for a specific density distribu-
tion. Problems of the type given by conditions (1) and (2)
normally have either no solution or a finite, often large
set of “basic” solutions (each involving exactly N
orbits) which in turn can be linearly combined with
positive weights to form a multiparameter family of
solutions. Thus, the question arises: To what degree
does this apparent nonuniqueness arise from the finite
cell approximation or is it inherent in the basic physical
problem ? In this investigation this uniqueness question
has been consistently ignored, giving first preference to
deriving one solution for at least one density distri-
bution.
The preceding description of the present procedure
makes it apparent that it does not use explicitly the
collisionless Boltzmann equation, quite in contrast to
most constructions of equilibrium configurations. This
difference is, however, purely a matter of technique:
In step C the orbits are computed by explicit use of the
equations of motion which are themselves the basis of
the Boltzmann equation. The parallelism is further
emphasized if one considers the set of C(I) B(I,J)
values as a representation of the phase space distri-
bution function,/, by a set of delta functions. For each
C(/) • £(/, J) value the space coordinates are directly
given by /, i.e., the coordinates of the center of cell 7,
and the velocity components are indirectly given by /
through the following circumstance. A specific orbit
(designated I) which has three effective integrals will pass a specific point designated J (or its immediate
neighborhood) with a specific velocity vector. (More
precisely, it does so with one of a small number of
distinct velocity vectors, a number which reaches 8 for
interior points of a three-dimensional box orbit; for
simplicity of the argument we shall, however, ignore
here this finite multiplicity.) Thus / designates in-
directly a velocity vector. Altogether then the set of
C(I) B(I,J) values has exactly the character of a
phase-space distribution function.
This parallelism brings out two further points.
(1) Equation (1) can now be looked at as an approxi-
mate representation of the ordinary equation for the
total density at a point expressed as an integral of /
over the three velocity components. (2) Similarly,
condition (2) represents here (with Æ > 0 by definition)
the equivalent to the general condition of non-
negativeness for / everywhere in phase space. This
general condition is an essential feature of any con-
struction of an equilibrium configuration, but for-
tunately in the form of condition (2) it can be taken
care of by the standard methods of linear program-
ming, though not without extensive calculations.
III. CONSTRUCTION OF A TRIAXIAL MODEL
The procedure described in the preceding section has
been applied to one specific case for two purposes:
first, to ascertain that the procedure does not encounter
unforseen major obstacles and, second, to add a new
argument for the existence of triaxial self-consistent
equilibrium configurations for stellar systems. The
description of the specific sample case in this section
will follow closely the description of the general
procedure in the preceding section.
a) Choice of Density Distribution
The density distribution of the dynamical con-
figuration to be constructed is chosen to have the semi-
analytical form
P(X, r,Z) = F{R) – G(R) • (2Z2 – X2 – Y2)I(2R2)
+ H(R) 3(X2 – Y2)/R2 . (3)
The first term is spherical. The second term has the
form of the spherical harmonic with 1 = 2 and m = 0.
It shortens the Z axis and lengthens the X and Y axes
equally. The third term has the form of the spherical
harmonic with 1 = 2 and m = 2. It lengthens the X
axis and shortens the Y axis. The major axis of the
figure will lie along the X coordinate, the median one
along the Y coordinate, and the minor one along the
Z coordinate.
To simulate elliptical galaxies Hubble’s profile, with
King’s modification, is chosen for the first term:
F(R) = (1 + R2)'312 . (4)
This function is regular at the center, the correspond-
ing mass diverges with R, but only logarithmically,
and the potential well has a finite depth. In equation
(4) the central density and the “core radius,” at which
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 No. 1, 1979 NUMERICAL MODEL FOR TRIAXIAL SYSTEM 239
the surface density drops to one-half, are normalized
to 1.
The two remaining functions, G(R) and H(R), are
used to give the equidensity surfaces fixed axis ratios
(simulating, but not reproducing, similar ellipsoids).
The values chosen for the ratio of the major to the
minor axis, a, and the ratio of the median to the minor
axis, ß9 are
a = 2.00 and 0=1.25. (5)
For lack of knowledge of any more relevant source,
these values were taken from Stark’s analysis (1977)
of the bulge of the Andromeda Nebula. The conditions
to be fulfilled are then
p(aZ, 0,0) = p(0,0, Z) and P(0,0Z, 0) =/>(0,0, Z).
(6)
To apply these conditions, one may first note that
along the axes equation (3) takes the simple forms
P(R, 0, 0) = F(R) + $G(R) + 3H(R) ,
p(09 R, 0) = F(R) + ¿GCR) – 3H(R) ,
/>(0, 0, R) = F(R) – G(R) . (7)
By combining equations (6) and (7) such as to elimi-
nate H(R) one obtains of the three functions F, G, and H at small as well as
large R:
F(R) = 1 – 1.5R2 + •••,
G(R) = 0.880952R2 + •••,
H(R) = 0.54762R2 + •••,
F(R) = l/R3 – 1.5/R5 + • • •,
G(R) = 0.726106/R3 + • • •,
H(R) = 0.276034/R3 + • • •. (9)
The density distribution here chosen has equidensity
surfaces which are not strictly ellipsoidal. In fact they
are strongly dimpled on the Z-axis and mildly so on
the Y-axis. This complicates the chosen figure and
makes this first model a somewhat more severe test
for the procedure as a whole, but it also produces very
simple – equations for the potential—a circumstance
which speeds up the orbit integrations.
As a final specification for the density distribution,
zero figure rotation is chosen—not as a particularly
likely case for elliptical galaxies but as a simplifying
choice for this first model.
b) Derivation of Gravitational Potential
Since the density distribution is expressed by
equation (3) in terms of three spherical harmonics, the
potential can be expressed in similar form:
<¿(X, Y,Z) = U(R) – V(R) (2Z2 – X2 – Y2)/(2R2)
F(aR) + F(ßR) – 2F(aßR) + W(R) 3(X2 – Y2)IR2 . (10)
= G(aR) + G(ßR) + G(aßR) . (8)
Since the left side of this equation is a known function
in consequence of equation (4), the only unknown in
this equation is G(R). Generally, however, equations
of this type do not have unique solutions or may have
no solution if regularity conditions both at small R and
large R are added. In the present case it turns out that
there is only one solution of equation (8) which is
developable in powers of R and, similarly, only one
which is developable in l/R. These two solutions turn
out not to be identical but to cross each other near
R = 0.6 and again near R = 1.0, and have moderate
differences between these crossings. The final choice
for G(R) is then the interior solution for R < 0.6, the
exterior solution for R > 1.0, and a weighted average
of the two solutions for 0.6 < Æ < 1.0 with the weights
arranged such that G(R) and its first derivative are
continuous everywhere. This compromise represents
just a minor modification of the axis ratio ß from its
originally chosen value of 1.25, but only in a limited
range of R and nowhere by more than 2%.
With F(R) and G(R) decided upon, H(R) can be
computed by successive use of the third of equations
(7), the first of conditions (6), and the first of equation
(7). It may be worth noting the asymptotic behavior The Poisson equation, with 47tG normalized to 1, then
gives the standard differential equations for the
functions U9 V9 and W:
1 d2(RU) __
R dR2
1 d\RV) 6
R dR2 R2 G,
1 d2(RW)
R dR2 -r>w=h- (11)
The first of equations (11), with F given by equation
(4), can be solved explicitly but not in a form con-
venient for numerical orbit computations. All three
equations (11) can be integrated numerically by stand-
ard methods and give unique solutions when account
is taken of the physical boundary conditions at R = 0
and at R -> oo.
The asymptotic behavior of U9 V9 and W is given
for small by
£/ = —!+ 0.166667R2 – 0.075R4 + • • •,
V = — 0.061064R2 + 0.062925R4 + • •,
W = -0.015403R2 -f 0.011054R4 + • • •, (12)
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 SCHWARZSCHILD Vol. 232 240
and for large R by
£/ = -(InR)IR – 0.69360/R – 0.25/R3 + • • •,
V = -0.12102/R + 0.7641/R3 + • •,
W = -0.04601/jR + 0.8711/iR3 + • • •. (13)
These equations show that (a) the depth of the
potential of this model, with the present normaliza-
tions, is exactly 1, and (b) the leading terms of
equations (12) present for R « 1 three independent,
perpendicular harmonic oscillators, each with its
separate energy integral—characteristic for any
potential which is regular at the bottom of its well.
(c) Equations (13) give for the present model, in spite
of its midly diverging mass, a spherical potential in
the limit of very large R.
For efficiency in the calculation of the acceleration
components AX, AY, and AZ, which are needed at
every step of an orbit integration, five auxiliary
functions are defined:
UD-m- ™ =
WD^^WIR*), (14)
VT = VIR , WT = W¡R .
With the help of these functions the following ex-
pressions for the accelerations can be derived by
differentiating equation (10) for the potential:
AX = -{X¡R) [UD + VD-i(X2 + Y2 – 2Z2)
+ WD 3(X2 – Y2) + VT + WT 6],
AY = —(Y/R) [UD + VD-}(X2+ F2 – 2Z2)
+ WD 3(X2 – Y2) +VT- WT 6],
AZ = -(Z/RMt/Z) + VD-i(X2 + r2 – 2Z2)
+ WD 3(X2 – Y2) – VT 2]. (15)
Thus, once the six one-dimensional functions, i.e., U
and the five auxiliary functions, are computed and
given in tabular form for a set of R values sufficiently
dense to permit linear interpolation, then the compu-
tation of the potential and the acceleration components
can be efficiently accomplished by equations (10) and
(15).
All the preceding equations are expressed in normal-
ized units which were fixed by setting the “core
radius” rc, the central density pc, and the constant of
gravitation 4ttG equal to 1. These normalized units are
related to the defining quantities as follows:
Length:
Mass:
Time: [200 pc],
[3 x 109 Mq] ,
llV(4*Gpc) [2xl05yr],
Velocity: rc\/(47rGpc) [1000 km s' (16) The numerical values in brackets correspond to the
data given by Miller and Prendergast (1962) for the
elliptical galaxy NGC 3379.
c) Computation of Orbits and Orbital
Density Distributions
Each orbit was computed by a standard second
order scheme not requiring iterations. For the time
step a value of 0.1 was persistently used. A small-
amplitude oscillation along the Z-axis has a period of
about 10 in the present model. Hence, even for the
fastest oscillations a period was covered by about 100
integration steps. The overall accuracy may be judged
by the constancy of the energy, which was found to
fluctuate by a few units in the fifth decimal (compared
to a potential well depth of 1) for an orbit typically
carried over a time span of about 104 and by a couple
of units in the fourth decimal for the very few orbits
computed for a time span of 105. The calculation of
an orbit of typical length took about 12 sec execution
time on an IBM 360-91 computer.
Before an orbit integration could result in a B{I, J)
table, as required for the present procedure, a scheme
for dividing the model into cells had to be decided
upon. The following scheme was adopted here. In the
radial dimension, space was divided into a central cell
and 15 surrounding shells. The 16 boundary surfaces
were not chosen to be spheres, but rather equipotential
surfaces. This has the obvious advantage that the key
question of whether an orbit is ergodic or not is
explicitly answered by its resulting B{I, J) table. The
boundary potential values were picked (as shown by
Table 1) such that the average R value increases from
surface to surface by a factor of 21/3. This has the
advantage that the chosen model density distribution
places approximately equal masses into the shells and
that the density varies in the radial direction within a
shell by no more than a factor 2. The latter is also true
for the central cell, which occupies the volume within
the innermost of the boundary potential surfaces. The
portion of the model density distribution falling out-
side the outermost potential surface of Table 1 was not
included in the representation by orbital density
distributions. This neglect, introduced for obvious
practical reasons, should have no more effect on the
results than to enhance the occupation numbers for
the high-energy orbits needed to represent the model
in shell 15 to make up for the density which would
TABLE 1
Potential Values of Cell Boundaries
<R> <R>
0.7937 0.91698 1.0000 0.88146
1.2599 0.83647 1.5874 0.78269
2.0000 0.72191
2.5198 0.65671
3.1748 0.58984 4.0000 0.52376 5.040.
6.350.
8.000.
10.079.
12.699. 16.000.
20.159. 25.398. 0.46046
0.40131 0.34712
0.29829 0.25488
0.21670
0.18344
0.15469
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 No. 1, 1979 NUMERICAL MODEL FOR TRIAXIAL SYSTEM 241
Fig. 1.—Pattern of 19 facets used as cell boundaries
have been placed into shell 15—and to a lesser degree
into shell 14, etc.—by the very high energy orbits
which would have been needed to represent the
neglected portion of the model. Thus the present
representation reaches out to an average R of about
25, i.e., to 5 kpc according to tabulation (16).
In the angular dimensions each octant of each shell
was divided into 19 facets, as shown by Figure 1. This
facet pattern was chosen to fulfill the following con-
ditions: (1) The pattern should look the same when
viewed from any of the three axes. (2) The facet
boundaries should be great circles passing through
one of the axes so that each is described by a fixed value
for the ratio of a pair of coordinates. (3) All 19 facets
should subtend the same solid angle as nearly as
possible. This last condition was fulfilled by an
optimum choice of the two free angles indicated in
Figure 1 (a = 16?9860, ß = 17?2153), resulting in a
2% maximum deviation from the average. By com-
bination of the shell boundaries with the facet bound-
aries each octant was divided into 286 cells, 19 in each
of 15 shells plus one central cell. The designating J
value for each cell was arbitrarily chosen by starting
with the 19 cells of the outermost shell (/ = 1 to 19)
and ending with the central cell (/ = 286).
The buildup of a B{I,J) table during an orbit
integration was carried out as follows. At the end of
an integration step a new X, F, Z point is reached.
The potential is computed from equation (10) and
compared with the boundary values of Table 1. Thus
the shell into which the new point falls is established.
Next the coordinate ratios X/ Y, Y/Z, Z\X are formed
and compared with the boundary values for the facet
pattern and the facet containing the new point is
found. Combining shell and facet gives the cell—and
its J value—into which the new point falls. Finally the
particular entry in the B{I,J) table for which /
designates the orbit being integrated and J the cell just
established is increased by one [appropriate if a
constant time step is used and efficient in providing
integer entries in the i?(/, J) tables].
Must eight B{I,J) tables, one for each octant, be retained for each orbit or does a single table per orbit,
based on the absolute values of the coordinates,
suffice? For the present model the latter is the case.
The density distribution of the model and, conse-
quently, its potential have eightfold symmetry (re-
flection symmetry in each coordinate). The same is not
true for the density distribution produced by each
single orbit. However, inspection of the equations of
motion shows that in consequence of the symmetry of
the potential each orbit belongs to a group of eight
orbits related to each other by a change of sign in one,
two, or three coordinates. (The eight orbits of one
group are often all or in part identical except for phase
shifts.) If the eight orbital B{I, J) tables are combined
with equal weights, a group B(I, J) table results which
has the eightfold symmetry of the model density. This
group £(/, J) table can be obtained from a single orbit
by ignoring the signs of the coordinates in the buildup
of the orbital B(I, J) table—a shortcut here persistently
used.
d) Orbit Families and Effective Integrals
During an initial orientation phase the two-
dimensional orbits in the Y-Z plane of the present
model were studied with the methods now commonly
used for two-dimensional systems (for an early clear
exposition, see Hénon and Heiles 1964). The results
from this phase agree with those typical for other such
systems, including the case of the Galaxy with its
“third integral” (Contopoulos 1960).
With this background the computation of approxi-
mately 1500 three-dimensional orbits was undertaken.
No attempt was made to study the full variety of orbit
types. Except in a few experimental cases all orbits
computed belong to two major families, box orbits and
Y-tube orbits. Box orbits are here defined as orbits for
which one may choose as starting condition an
arbitrary point in space but a definite velocity
namely, zero. Such orbits were found to have a vanish-
ing average angular momentum—when averaged over
a large number of oscillations—as might be expected
from the reversal of direction when the orbit passes
through the immediate neighborhood of the initial
stationary point. Y-tube orbits are here defined as
orbits for which one may choose as starting conditions
arbitrary X and Z coordinates and an arbitrary Y
velocity component but with the other three phase-
space coordinates set to zero. To avoid overlap
between the two families one should add to the de-
finition of Y-tube orbits that their average angular
momentum should have a nonvanishing X component.
Clearly the initial values of these orbit families
contain three independent parameters. The key
question then is whether the resulting orbital density
distributions also form three-parameter families. The
computed B(I,J) tables have given a decisively
affirmative answer. Only a small fraction of the B(I, J)
tables leave open the possibility of ergodicity of the
corresponding orbits. Furthermore, the intriguing
issue of whether some of the computed orbits might
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 i 1 1 n —i 1 r
Fig. 2
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 NUMERICAL MODEL FOR TRIAXIAL SYSTEM 243
have two rather than three or one integral (Froeschle
and Scheidecker 1973; Contopoulos, Galgani, and
Giorgilli 1978) has not been studied here. In spite of
these unsolved items, one may conclude that, for the
model here chosen, the majority of orbits belonging to
the two major families of orbits covered by the present
integrations have three effective integrals, this term
being used here as defined in § lid. Thus the central
necessary—though not sufficient—condition for con-
structing a self-consistent dynamical model for the
chosen density distribution is fulfilled. (In view of the
relatively arbitrary choice of this density distribution
it would seem unlikely that it would be one of the
special cases, if any, for which this central condition is
not necessary.)
For one typical box orbit the existence of three
effective integrals has been ascertained by a special
study, additional to the inspection of its Æ(7, J) table.
For this orbit the velocity vector has been printed out
every time the orbit passed through a tiny cell surround-
ing a preselected arbitrary point. (Passages through the
seven mirror images of the cell were also exploited.)
To obtain a sufficient sample of passages (actually 18),
the integration had to be carried for a time span of
nearly 106, i.e., 2 x 1011 yr. The 18 velocity vectors
thus found did not scatter randomly over the sphere
prescribed by the energy integral of the orbit, but fell
on eight points on the sphere with two points of each
of four pairs lying diametrically opposite to each
other. This is exactly what should be expected for an
orbit with three effective integrals, with the number
eight being appropriate for an interior point of a three-
dimensional box orbit—just as for an orbit in a three-
dimensional harmonic potential with irrational period
ratios.
Figures 2 and 3 present an attempt to give a semi-
graphical picture of one selected orbit for each of the
two major families here investigated, always a
difficult task for three-dimensional orbits covering
many oscillations (Froeschle 1970). Here each orbit is
represented by three perpendicular cross sections
through its orbital density distribution, based on its
J5(7, f) table. In both figures all three cross sections
show distributions clearly not filling the total space
out to the energetically permitted potential surface.
This shows that both orbits are not ergodic. The box
orbit given in Figure 2 was started with zero velocity
at a point close to the Y — Z plane with exactly equal
Y and Z coordinates. Despite this initial symmetry in
Y and Z, the Y — Z cross section of Figure 2 indicates
a strong asymmetry with respect to Y and Z for the
orbit as a whole. This shows that even the rather
moderate difference in the Y and Z axes of the chosen
density distribution strongly affects box orbits. In this
regard the present model differs significantly from,
say, prolate spheroids. The box orbit of Figure 2 has eight corners (the
starting point and its seven mirror images) which lie
close to the Y-Z plane, and their projections onto
that plane are indicated by four crosses. The sides of
the box containing the orbit are far from plane. As
Figure 2 indicates, the two sides crossing the Z axis
are highly concave, while those crossing the X axis are
strongly convex. In fact, the X- Y cross section of this
box orbit shows a substantial orbital density near the
X axis quite far from the center, even though the zero-
velocity starting point lies close to the Y-Z plane. This
tendency of box orbits to occupy the neighborhood of
the major (X) axis has a significant effect on the final
step of the present procedure.
Turning to the Y-tube orbit shown in Figure 3, it
should be stated that it was selected to show most
strikingly the main characteristics of this orbit family;
most X-tube orbits enter more cells than the selected
one does. All three cross sections of Figure 3 show
that the orbit nowhere comes close to the X-axis. It
may bè approximately described as circling around the
X-axis while simultaneously executing oscillations in
the X coordinate—thus forming a tube around the
X-axis. The X-Z and X-Y cross sections show that
the tube has its greatest diameter at the extremes of the
oscillations in X where it also spends more time than
in its narrow throat at X = 0. This orbit is unusually
thin-walled. Typical X-tube orbits have thicker walls,
i.e., they may reach far away from the X-axis at any
X value but may range simultaneously close to the
X-axis, thus entering more cells than the example of
Figure 3. Accordingly this orbit family, just as the box
orbits, has three significant geometrical parameters,
which for this family can be taken as average diameter,
length, and wall thickness.
e) Model Construction by Linear Programming
The final step of the present procedure is to re-
produce the chosen model density distribution by an
appropriate superposition of orbital density distri-
butions, as formulated by equations (1) and (2). Even
though modern methods in linear programming
provide straightforward techniques for solving prob-
lems like the one here encountered, in practise this last
step has a trial-and-error character. The reason is that
the blind inclusion of an arbitrarily large number of
orbits would lead to prohibitively long computations.
Hence insight has to be developed by trial and error
regarding the types of orbits and ranges of parameters
required by the solution so that the number of orbits
involved can be kept within reasonable bounds. The
following description will concentrate on the main
operations required to arrive at the final solution, with
the trial method indicated at the end.
Three sets of orbits were computed: “centered”
box orbits, “shifted” box orbits, and X-tube orbits.
Fig. 2.—Three cross sections through a selected box orbit. Only the cells in the inner nine shells are shown. The potential surface
corresponding to the energy of the orbit is indicated by the small circles. The number in each cell represents the fraction of time
(scale arbitrary) which the orbit spends in the cell (one, two, or three dots stand for the corresponding numbers). Cells without
numbers are not entered by the orbit. The unentered cells located within the limiting potential surface indicate the existence of
additional integrals for the orbit.
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 NUMERICAL MODEL FOR TRIAXIAL SYSTEM 245
The starting points for each of 285 centered box orbits
(one each for the 19 cells in each of 15 shells, exclusive
of the center cell) was positioned near the center of the
outer, high-potential boundary facet of the cell. The
starting point for each of 285 shifted box orbits (again
one for each cell except the center one) was positioned
near that corner of the outer boundary facet of the cell
which produced the highest orbital density in that cell.
Some of the box orbits were excluded from further use
because their £(/, J) tables indicated that they happen
to fall near a closed resonance orbit. In such a case an
orbit drifts only slowly through its permitted space and
would require, say, 1000 oscillations, rather than the
100 oscillations normally covered here, to obtain a
reliable 2?(/, J) table (as was checked in one test case).
Early trials in which these two sets of box orbits were
used exclusively, failed. This failure turned out to be
caused by the fact that all box orbits in this model
place a relatively high density into the cells near the
Y-axis. In consequence, when occupation numbers
C(/) are determined for box orbits so as to reproduce
the model density in cells far from the Y-axis, those
near the Y-axis are found to be overpopulated (this
situation parallels the sharp density rise in the center
of a spherical system when only orbits with zero
angular momentum are considered).
Accordingly, a set of approximately 800 Y-tube
orbits were computed since this type of orbit leaves the
neighborhood of the Y-axis unoccupied. The starting
points of these orbits were placed in the Y-Z plane in
an array of 60 points, 15 each along four lines deviating
from the Y-axis by 25?6, 45?0, 64?4, and 81?5, re-
spectively. The Y component of the initial velocity was
varied over a substantial range and 131 orbits,
averaging nine in the energy range of each shell, were
selected for inclusion into the linear program. Pro-
tracted trials with this enlarged set of orbits still
failed to give a solution. This time the cause was found
to be an overpopulation in a narrow cone of cells close
to the Z-axis in the five innermost shells. This over-
population was produced by orbits required for the
intermediate shells 6-10. This difficulty may have been
caused by the dimples around the Z-axis in the chosen
model density distribution. Instead of changing the
model, which would have impaired the basic pro-
cedure that starts with the “arbitrary” choice of the
model density, a search for additional appropriate
orbits was carried out. It was found that Y-tube orbits
with starting points along a line in the Y-Z plane and
inclined 59?5 to the Y-axis might fill the need. About
100 such orbits in the energy range of shells 5-10 were
computed, and six of these were selected by the crite-
rion of minimum density in cells directly adjacent to
the Z-axis but maximum density in cells just a little
farther away from this axis. One of these orbits, which
are quite thin-walled, is shown in Figure 3. After
adding these six orbits to the linear program, a
solution was obtained.
The character of this solution will be described in
the next subsection. Here follows a discussion of non-
standard items in the linear programming.
For the basic equation (1) the necessary input, £(/, J), for the left side is provided by the orbits just
discussed. The total number of orbits finally included
was M = 684 (547 box orbits and 137 Y-tube orbits).
The number of cells considered was N = 285. This
excludes the central cell for which, being directionally
unsegmented, any orbit of sufficiently low potential
not to enter the first shell can provide the required
density (the orbits required by the 285 other cells do not
overcrowd the central cell). The data required for the
right side of equation (1), D(J), are obtained by inte-
grating the model density distribution given by
equation (3) over each cell volume.
When the linear programming problem here en-
countered was shown to Dr. Harold Kuhn, he promptly
devised the following computer-time-saving trans-
formation. The M orbits are divided into N “basic”
orbits and M — N “nonbasic” orbits. The basic
orbits can be arbitrarily selected, except for the con-
dition that they must be linearly independent so that
the inverse of their B matrix, J9-1, exists. Equation (1)
can then be written in the form
N M
y C(J). B{1 J) + V C(7) • 5(7, J) = D(J)
J = 1 I = N + 1
for 7 = l9 N. (16)
Next, the nonbasic orbits as well as the D(J) vector
are expressed as linear combinations of the basic
orbits,
= 2 A(I,K) B(K,J) K = 1
for I = N + l, M; and J=l,N; (17)
D(J) = 2 P(K) B(K,J) for J = l,N. (18)
K^l
With the help of definitions (17) and (18) for A(I, K)
and P(Y), equation (16) can be transformed by
multiplication with B'1 into
M
C(J) + 2 C(7) • A(J’ •/) = P(<J) for / = 1, iV.
I = N +1
(19)
Next, Dr. George Soules of the Institute for Defense
Analysis at Princeton devised a code including
Kuhn’s transformation as well as the pivot operation
of the Simplex method in so efficient, transparent, and
flexible a style that the inexperienced writer could use
and adapt it effectively throughout the trial-and-error
and final solution phases, though with continued
guidance by Dr. Soules.
Turning finally to the trial-and-error method here
employed, a special feature of the present problem
needs emphasis. Any orbit can enter only those cells
in which the potential does not exceed the orbit’s
energy. Since the chosen cells can be divided into 15
subsets of 19 cells with common equipotential
boundaries, the present problem can be replaced, at
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 246 SCHWARZSCHILD Vol. 232
least in principle, by a sequence of 15 smaller problems,
as follows. For the first small problem take only the
first 19 equations (1), i.e., those which refer only to the
cells in the outermost shell and use only those of all
the available orbits which enter this shell. Solve this
restricted problem by linear programming, i.e., by the
same techniques discussed for the total problem.
Assuming the first small problem to have solutions,
pick one. Compute a reduced model density for each
cell by subtracting from each of the original D(J)
values the 19 values of C(I) B(I,J) which according
to the picked solution are not zero. The reduced model
densities will be exactly zero in the cells of the outer«
most shell and somewhat smaller than their original
values in all shells further in, since the 19 orbits re-
quired to reproduce the outermost shell will also
produce some density in the shells farther in. For the
second small problem use the reduced D(J) values,
only the 19 equations (1) referring to the cells of the
second shell (counting here from the outside in) and
only those orbits which enter the second but not the
first shell, etc.
If a sequence of several such small problems has
given solutions, shell by shell, and the next small
problem gives “no solution,” this may be caused by
one of two failures: a poor choice of the solutions
picked for the outer shells, or a lack of one or more
required orbits among those included in the calcula-
tions. Which of these failures is active in a specific case
can be determined by application of the large flexible
linear programming code at all the shells thus far
solved for, as well as the one for which no solution was
found. It will either give one solution (or more) and
its new reduced model density—and thus permit con-
tinuation of the stepwise procedure—or will give “no
solution.” In the latter case, additional orbits have to
be prepared. In such a case the search for useful
additional orbits is much helped by a study of the last
available reduced model density distribution. Specific-
ally, its low points indicate the cells which the orbits
to be added should enter as rarely as possible.
This shell-by-shell procedure was actually used for
the present model and led, with minor backing and
filling, to consistent solutions for the outer 10 shells,
but to a prolonged impasse at the next shell. This in
turn led to the addition of the six thin-walled X-tube
orbits described in the preceding subsection. After
this addition the large code gave a solution for the total
system on the first trial.
/) Characteristics of the Constructed Model
The main feature of the numerical model here com-
puted for a self-consistent triaxial stellar system is that
it exists. The existence of such a numerical model does
not prove with mathematical safety the existence of a
corresponding true solution; 2280 separate cells (285
in each of eight octants) clearly permit only a marginal
representation of the distribution of, say, 1011 stars of
an elliptical galaxy. Nevertheless, the existence of such
a numerical model should, it would seem, strongly
suggest the existence of triaxial dynamical equilibrium
configurations. Next, the question of uniqueness of the solution
might be raised. Clearly the numerical solution for the
problem in the simple form of equations (1) and (2) is
not unique; this type of problem, if it has any solution,
has normally a substantial number of basic solutions.
However, this nonuniqueness of the numerical
solution might well just reflect the indefiniteness intro-
duced into the problem by the use of a very finite
number of cells rather than indicate nonuniqueness
of the exact solution. On the other hand, one might
well wonder whether nonuniqueness of the exact
solution is not suggested by the fact that a numerical
solution was found here even though orbits of only
two types were included while other types, such as
Z-tube orbits which have been found to exist, were
excluded. This basic uniqueness question is thus left
unanswered by the present investigation.
The final item to be discussed is the rotation charac-
teristics of the present model. At the outset of this
investigation, zero rotation of the density figure was
chosen for this particular model. Hence only the
internal stream motions need be considered here. The
specific numerical solution obtained turns out to
utilize 230 box orbits and 55 X-tube orbits. Box orbits
can not contribute to stream motions because they
reverse direction whenever they come close to a zero-
velocity point (i.e., one of the corners of the box) and
thereafter repeat (in the case of zero rotation of the
potential figure) the same orbit but in the opposite
direction. Accordingly a box orbit will pass any point
in its box statistically as often in a specific direction as
in the exactly opposite direction, thus nullifying its
contribution to any stream motion at any point. The
same is not true for Z-tube orbits. An orbit of this
type moves around the Z-axis persistently in the same
sense. This might be considered a consequence of the
fact that Z-tube orbits are reasonance orbits, not in
the sharp sense of closed orbits but more in the wider
sense that for each such orbit the periods of the oscil-
lations in Y and Z are exactly equal if averaged over a
large number of oscillations. This persistence of sense
of motion around the Z-axis clearly permits Z-tube
orbits to cause stream motions of the character of a
rotation around the Z-axis, the longest axis. Here the
usual ambiguity appears, e.g., Ruiz and Schwarzschild
(1976), whether all stars occupying one orbit traverse
the orbit in the same direction or whether both direc-
tions are occupied. Maximum rotation rate is obtained
if all stars on Z-tube orbits are chosen to move around
the Z-axis in the same sense, while zero rotation results
from equal occupation of both directions. The stability
problem raised by the latter choice, as well as all other
stability problems, has not been included in this
investigation.
One may then conclude that the present numerical
model can be rotating, by internal stream motions,
though not by figure rotation, around its long Z-axis
but not around the 7- and Z-axes (for a speculative
application, see Bertola and Galletta 1978). This last
restriction might be lifted by the inclusion of Z-tube
orbits—and 7-tube orbits if they exist. However, a
more substantial change might be achieved by con-
© American Astronomical Society • Provided by the NASA Astrophysics Data System

197 9ApJ. . .232 . .2363 No. 1, 1979 NUMERICAL MODEL FOR TRIAXIAL SYSTEMA 247
structing another model for which the present density
distribution may be chosen except for one major
change, namely, nonzero figure rotation. It is the plan
of the writer to extend this work in the last mentioned
direction.
Since most of the dynamical aspects of this investi-
gation and all of the required techniques of linear
programming were new to me, help from experts was
of extraordinary importance for this work. The main
parts of such help I record here with the most earnest
gratitude. Dr. Douglas Heggie has given me extensive
instructions in stellar systems dynamics as well as introducing me into the relevant literature. Dr. James
Binney has given me the full advantage of his ex-
perience with three-dimensional systems. Dr. Harold
Kuhn has introduced me to linear programming and
even devised a special transformation for the type of
problem on hand. Dr. George Soules of the Institute
for Defense Analysis has devised a code which is
simultaneously efficient and transparent, has tested the
code with the initial data of the present problem on the
computer of his institution, has translated the code to
fit the Princeton University computer, and has guided
me all the way to the final solution. This research was
supported by NSF grant AST 75-22924.
REFERENCES
Aarseth, S. J., and Binney, J. 1978, M.N.R.A.S., 185, 227.
Bertola, F., and Galletta, G. 1978, Ap. J. {Letters), 226, LI 15.
Binney, J. 1978, Comments Ap., Vol. 8, No. 2.
Contopoulos, G. 1960, Zs.f. Ap., 49, 273.
Contopoulos, G., Galgani, L., and Giorgilli, A., 1978, Phys. Rev., A18, 1183. Freeman, K. C. 1965, M.N.R.A.S., 134, 1.
Froeschle, C. 1970, Astr. Ap., 4, 115. Froeschle, C., and Scheidecker, J. P. 1973, Ap. Space Sei., 25,
373.
Gott, J. R. 1977, Ann. Rev. Astr. Ap., 15, 235. Hénon, M., and Heiles, C. 1964, A.J., 69, 73.
Hunter, C. 1977, A.J., 82, 271.
Illingworth, G. 1977, Ap. J. {Letters), 218, L43.
Lynden-Bell, D. 1962, M.N.R.A.S., 123, 447.
. 1963, M.N.R.A.S., 124, 95.
Miller, R. H., and Prendergast, K. H. 1962, Ap. J., 136, 713. Miller, R. H., and Smith, B. F. 1979, Ap. J., 227, 407.
Ruiz, M. T., and Schwarzschild, M. 1976, Ap. J., 207, 376.
Schechter, P., and Gunn, J. E. 1978, private communication. Stark, A. A. 1977, Ap. J., 213, 368.
Vandervoort, P. O. 1978, Bull. AAS, 10, 483.
Martin Schwarzschild: Department of Astrophysical Sciences Peyton Hall, Princeton University, Princeton,
NJ 08540
© American Astronomical Society • Provided by the NASA Astrophysics Data System