A&A 471, 419432 (2007)
DOI: 10.1051/00046361:20077672
Dynamical models with a general anisotropy profile^{}
M. Baes

E. Van Hese
Sterrenkundig Observatorium, Universiteit Gent,
Krijgslaan 281 S9, 9000 Gent, Belgium
Received 19 April 2007 / Accepted 28 May 2007
Abstract
Aims. Both numerical simulations and observational evidence indicate that the outer regions of galaxies and dark matter haloes are typically mildly to significantly radially anisotropic. The inner regions can be significantly nonisotropic, depending on the dynamical formation and evolution processes. In an attempt to break the lack of simple dynamical models that can reproduce this behaviour, we explore a technique to construct dynamical models with an arbitrary density and an arbitrary anisotropy profile.
Methods. We outline a general construction method and propose a more practical approach based on a parameterized anisotropy profile. This approach consists of fitting the density of the model with a set of dynamical components, each of which have the same anisotropy profile. Using this approach we avoid the delicate finetuning difficulties other fitting techniques typically encounter when constructing radially anisotropic models.
Results. We present a model anisotropy profile that generalizes the OsipkovMerritt profile, and that can represent any smooth monotonic anisotropy profile. Based on this model anisotropy profile, we construct a very general sevenparameter set of dynamical components for which the most important dynamical properties can be calculated analytically. We use the results to look for simple onecomponent dynamical models that generate simple potentialdensity pairs while still supporting a flexible anisotropy profile. We present families of Plummer and Hernquist models in which the anisotropy at small and large radii can be chosen as free parameters. We also generalize these two families to a threeparameter family that selfconsistently generates the set of Veltmann potentialdensity pairs. These new analytical models are an important step forward compared to isotropic or OsipkovMerritt models and can be used to generate the initial conditions for realistic simulations of galaxies or dark matter haloes.
Key words: Galaxy: kinematics and dynamics  cosmology: dark matter  methods:
analytical
1 Introduction
Numerical simulations, in particular Nbody or hydrodynamical
simulations, have become a major tool to study the structure,
dynamics, stability and evolution of stellar systems and dark matter
haloes. While at first sight this might seem to imply that analytical
studies of spherically symmetric dynamical models have become less
interesting, actually the opposite conclusion should be made. The
construction of realistic and simple dynamical models is utterly
important, as these models act as a reference frame or starting point
from which to generate the initial conditions for numerical
simulations. In this context, it is important to stress that the
initial conditions need to be generated from the correct distribution
function
which completely determines the dynamical
model. Kazantzidis et al. (2004) have recently demonstrated that the
correct details of the velocity distribution should be taken into
account. Simple adhoc Jeans dynamical models, i.e. models in which
the kinematics are assumed to be Gaussian distributions with velocity
dispersions derived from solving the Jeans equation, are not
sufficient and can lead to erroneous conclusions. The quest for
analytical dynamical models, which are simple enough on the one hand
and realistic enough on the other hand, is hence still important.
In the past few years, many dynamical models have been proposed that
are generated by various spherically symmetric potentialdensity
pairs. The earlydays models mainly represented systems with a
constant density core, such as the Plummer model or the isochrone
sphere (Plummer 1911; Hénon 1960). It now
appears that many dynamical systems such as galaxies and dark matter
haloes contain a central density cusp; also for such models a number
of representative potentialdensity pairs have been constructed and
distribution function have been derived (Hiotelis 1994; Tremaine et al. 1994; Buyle et al. 2007; Baes & Dejonghe 2004; Dehnen 1993; Hernquist 1990; Baes & Dejonghe 2002; Zhao 1996; Jaffe 1983).
Although the construction of analytical distribution functions hence
has continued to develop, most simple models still show an important
shortcoming: usually the distribution function can only be derived
under strict and unrealistic conditions on the anisotropy of the
velocity distribution. For many potentialdensity pairs, analytical
distribution functions have only been derived under the assumption of
isotropy. In this case, the phasespace distribution function that
supports any density profile can be found as a simple integration
using the famous Eddington formula. There are a few popular
generalizations of the simple isotropic case that lead to dynamical
models with an anisotropic velocity dispersion. Distribution functions
with a constant anisotropy profile can be generated with a formula
that is very similar to the Eddington formula
(An & Evans 2006a; Baes & Dejonghe 2004; Dejonghe 1986). Another special category are the
OsipkovMerritt models, engineered independently by
Osipkov (1979) and Merritt (1985). In these
essentially onedimensional models, the distribution function depends
on the energy and angular momentum integrals only through a linear
combination of both. The anisotropy profile for these models is
peculiar in the sense that the velocity dispersion tensor is
isotropic at small radii and becomes completely radial at large
radii. The OsipkovMerritt models were extended by
Cuddeford (1991) to models in which the velocity dispersion
tensor has an arbitrary anisotropy at small radii and becomes
completely radial at large radii. The transition formulae for all
these generalizations of the Eddington formula are sufficiently simple
that analytical distribution functions of these kinds can be
constructed for many simple potentialdensity pairs.
Whereas these models are a step beyond isotropic dynamical models,
they still do not correspond to the anisotropy that is observed both
in numerical simulations and real galaxies. Realistic dynamical
systems typically have a tendency towards moderate to strong radial
anisotropy at large radii. Cosmological simulations generally yield
dynamical structures that are far from isotropic. Dark matter
simulations in a cosmological CDM framework typically result in haloes
in which the anisotropy gradually increases outward to levels of
at the virial radius (Fukushige & Makino 2001; Diemand et al. 2004; Cole & Lacey 1996; Colín et al. 2000). In
general, there appears to be a connection between the logarithmic
density slope
and ,
both for pure dark matter haloes and for structures containing dark
matter and baryons (Hansen & Moore 2006). This connection was
argued to be a natural consequence of the Jeans equation
(Dehnen & McLaughlin 2005) and it was shown numerically to be an
attractor (Hansen & Stadel 2006).
In hydrodynamical simulations in which gas and stars are taken into
account, galaxies also typically have a significant radial anisotropy
at large radii (e.g. Oñorbe et al. 2007). Observational
dynamical evidence appears to support these trends. For example,
detailed stellar dynamical modelling of a large set of elliptical
galaxies by Kronawitter et al. (2000) indicates that these galaxies
generally contain an anisotropy profile that increases from near
isotropy in the central regions towards a significant radial
anisotropy at larger radii.
The anisotropy at small radii on the other hand is more subtle and
depends largely on the dynamical processes that shape the nucleus of
the system. In particular, the influence of a supermassive black hole,
now believed to be present in nearly all hot stellar systems, can have
a lasting influence on the central anisotropy. Models that contain a
supermassive black hole that grows adiabatically by slow accretion of
material are characterized by an isotropic to slightly tangential
anisotropy profile in the central region with
(Quinlan et al. 1995). On the other hand, models with a binary
supermassive black hole, formed by spiralling in due to dynamical
friction, typically have much more outspoken tangential anisotropy in
the central regions, with
(Quinlan & Hernquist 1997). This enhanced tangential anisotropy
results from the preferential ejection of stars on radial orbits
during the hardening of the black hole binary. Observationally, the
study of the anisotropy profile at very small radii is difficult and
hampered by both limited resolution effects and strong projection
effects. Gebhardt et al. (2003) found that the most massive
galaxies are strongly biased towards tangential orbits in the
innermost regions, whereas the lowermass galaxies have a range of
central anisotropies.
From this body of evidence it is clear that we need to extend our set
of simple isotropic or OsipkovMerritt dynamical models to models in
which the anisotropy profile has a more general shape, typically
rising from isotropy or moderate tangential anisotropy at small radii
to significant radial anisotropy at large radii. In the general
anisotropic case, the transition formulae between density and
distribution function are rather cumbersome (although mathematically
elegant) and involve LaplaceMellin integral transforms
(Dejonghe 1986). The goal of the present paper is to
describe a method to construct dynamical models with an arbitrary
potentialdensity pair and an arbitrary anisotropy
profile. We describe the general approach and the mathematical
formulation in Sect. 2. In the next two sections
we describe a more practical approach to construct such models: in
Sect. 3 we propose a general parameterized anisotropy
profile and in Sect. 4 we construct a general
sevenparameter family of dynamical components based on this
parameterized anisotropy profile for which the most important
dynamical properties can be calculated analytically. In
Sect. 5 we apply the results from the previous
sections to construct a family of Plummer and Hernquist models with a
smooth anisotropy profile with arbitrary values at small and large
radii. We then extend these two families to a threeparameter family
of generally anisotropic models that support the Veltmann set of
potentialdensity pairs. Finally, our results are summarized in
Sect. 6.
2 Construction of general anisotropic models
In spherical symmetry, the distribution function depends only on two
isolating integrals: the binding energy
and the magnitude of
the angular momentum L per unit mass,

(1) 

(2) 
In these expressions
is the positive binding potential of
the system, connected to the density through Poisson's equation

(3) 
and
is the transverse velocity

(4) 
When we have an expression for the distribution function, all
interesting dynamical properties can be derived. In particular, the
moments of the distribution function can be calculated as

(5) 
If we want to construct dynamical models for a given density profile
in practice, it is easier to consider an alternative characterization
of a spherical dynamical model, namely the framework of the augmented
densities. An augmented density is a bivariate function
,
which expresses the density as an explicit
function of both the radius r and the gravitational potential .
One can demonstrate that each augmented density function
completely determines a distribution function
and vice versa. There are several transition formulae
between these two equivalent formalisms, including an elegant
formalism based on combined LaplaceMellin integral transforms.
One of the advantages of the augmented density formalism is that it is
straightforward to construct selfconsistent dynamical models that
correspond to a given density profile. Indeed, we just have to solve
Poisson's equation, and construct any augmented density
that satisfies the relation

(6) 
Obviously, there are always infinitely many augmented density
functions that satisfy Eq. (6) for a given density
profile .
These different dynamical models have different
moments of the distribution function, and we often are looking for
those particular dynamical models which have a certain anisotropy
profile, defined as

(7) 
It appears that the augmented density formalism of dynamical models is
a very useful way to achieve this goal. Indeed, another advantage of
this formalism is that the moments of the distribution function can be
derived from the augmented density by subsequent derivations and
a single integration,

(8) 
where D_{x}^{m} denotes the mth differentiation with respect to x. In
particular, the radial and transverse velocity dispersions can be
found from the density through the relations



(9) 



(10) 
Now consider an augmented density that is a separable function of
and r,

(11) 
For such models, the dispersion profiles read
such that

(14) 
For models with a separable augmented density, the radial velocity
dispersion depends only on the dependent part of the augmented
density, whereas the anisotropy depends only on the rdependent part.
Equation (14) offers us in principle the opportunity to
construct dynamical models with an arbitrary density profile and an arbitrary anisotropy profile. First we solve
Eq. (14) for g(r). Next we determine the
gravitational potential
by solving Poisson's equation, we
invert this relation as
and we set
then the augmented density

(18) 
defines the desired model.
Although we have defined a way to construct an augmented density for a
dynamical model with an arbitrary density and anisotropy profile, this
is of limited use in practice. Indeed, the quantity we really want is
the distribution function for the model. In principle, the
distribution function for each
can be recovered
with the standard LaplaceMellin integral transform formulae. However,
the augmented density that results from an arbitrary choice of
and
will often be too complicated to allow an
analytical integral transform. As a result, we have to compute the
distribution function by numerical means. Unfortunately, the numerical
inversion of Eq. (5) to obtain the distribution
function from the augmented density is in general numerically
unstable. Indeed, since the density is an integration of the
distribution function over velocity space, the density is generally
much smoother than the distribution function. The inversion procedure
of determining f from
hence has the tricky job of unsmoothing
the information contained in the smoothed density. It hence comes as
no surprise that this operation is numerically unstable. For a more
precise mathematical demonstration of the unstable character of the
inversion formulae, we refer to Dejonghe (1986).
A more promising strategy is to determine a general parameterized
model for
that leads to a relatively simple function g(r). For each function ,
we can subsequently construct a
dynamical model
that has the same anisotropy
profile. Now consider a set of base functions
that are
sufficiently simple such that the augmented density
can be converted
analytically to a distribution function
.
It is easy
to see that any linear combination of such dynamical models will again
define a dynamical model with the same anisotropy profile since

(19) 
So if we can find a suitable linear combination of the base functions
such that

(20) 
we immediately obtain the required distribution function

(21) 
Finding the best linear combination of the components such that the
condition (20) is satisfied can in principle be
achieved by making a formal series expansion and determine each of the
expansion coefficients for each component. In practice however, it is
more straightforward to approximate the dynamical model by a
minimization routine. This could for example be achieved by
means of quadratic programming, in which a boundary condition is set
that forces the distribution function to remain positive in the entire
phase space. Practical realizations of this technique are already in
use by various dynamical modellers (Kuijken & Merrifield 1993; Merrifield & Kuijken 1994; Dejonghe 1989; Gerhard et al. 1998).
The remaining steps we still have to take are hence (1) create a
general parameterized anisotropy profile
that is flexible
enough to represent an array of realistic anisotropy profiles and
that still yields a simple g(r) function, (2) provide a suitable set of
base functions or components
for which, in combination
with the g(r) function derived above, the corresponding distribution
function can be computed analytically. We will look into these two
steps in the next two sections.
3 A parameterized anisotropy profile
As mentioned before, very few analytical dynamical models are known
that have a realistic anisotropy profile. Most are either completely
isotropic, have a constant anisotropy or have an anisotropy profile
according to the OsipkovMerritt type (isotropic in the centre and
completely radially anisotropic at large radii). In this paper, we aim
at the construction of spherical dynamical models with a more general
and realistic anisotropy profile. The most realistic approach would
involve an anisotropy profile that depends explicitly on the density
slope, as described in the Introduction. However, this would in
general not yield an augmented density that allows an analytical
distribution function. Instead, we will create dynamical models with
an anisotropy profile that is an explicit function of the radius
r. As long as the anisotropy profiles have enough freedom, this will
in practice enable us to construct dynamical models in which the
anisotropy and density slope are at least qualitatively connected. Our
specific goal is to create dynamical models in which the anisotropy
changes monotonically from an arbitrary value
at the central
regions to another arbitrary value
at large radii,
without a priori limitations on the values of
or
.
To reach this goal, we need to provide a functional
form that has the desired properties at small and large radii and that
can be inverted to give a reasonable simple form for g(r). Finding
such a form for
is not obvious. A straightforward candidate
that has the desired asymptotic behaviour is an exponential profile

(22) 
Solving for Eq. (14) yields, apart from an arbitrary
normalization factor,

(23) 
It comes as no surprise that such a model will not result in a simple
augmented density that can be inverted to an analytical distribution
function. In order to find a more suitable candidate function for
we inspire ourselves on the OsipkovMerritt models, where
the anisotropy profile reads

(24) 
resulting in isotropy in the central regions and complete radial
anisotropy at large radii. The resulting solution for g(r) reads

(25) 
Cuddeford (1991) generalized the OsipkovMerritt models by
considering models with an anisotropy profile

(26) 
These models also become completely radial at large radii, but at
small radii they can assume an arbitrary anisotropy .
If we
introduce this expression into Eq. (14) we obtain the
relatively simple expression

(27) 
The obvious way to generalize these results such that the anisotropy
does not have to become completely radial at large radii is to
consider the profile

(28) 
which corresponds to

(29) 
In fact, an even more general form for the anisotropy that yields a
very similar g function is

(30) 
with ,
which corresponds to

(31) 
with

(32) 
The couple (30)(31) is a suitable choice for our
goals. On the one hand this fourparameter model anisotropy contains
enough flexibility to represent any monotonically varying anisotropy
profile (with increasing or decreasing anisotropy). The parameters
have straightforward effects:
and
set the
anisotropy at small and large radii respectively,
determines
the typical transition radius between the two regimes and sets the sharpness by which this transition takes place. On the other
hand, the resulting g(r) function is sufficiently simple to lead to
analytical dynamical components for suitably chosen functions.
4 A library of components
The missing link in our approach is a set of base functions or
components with which we can construct a linear combination. We need
to look for functions
such that the augmented density

(33) 
gives rise to an analytical distribution functions
.
We will drop the subscript
in the remainder
of this section in order not to overload the notations.
The most obvious candidate components are those in which
is
a power law,

(34) 
where we limit ourselves to systems with a finite potential well
.
A finite total mass is obtained if
.
In order to calculate the distribution function
corresponding to this augmented density, we apply the general
MellinLaplace inversion formulae. After some algebra (see
Appendix A), we find that it can be expressed as a Fox
Hfunction:

(35) 
As we show in Appendix A, a sufficient (but not
necessary) condition to yield a welldefined distribution function,
i.e. continuous and nonnegative, is
or
and
.
For practical purposes, a series expansion is more useful. After some
lengthy algebra, one obtains the fairly simple expression (see
Appendix B)

(36) 
Apart from the distribution function, most other interesting dynamical
properties can be calculated analytically for this set of
components. Particularly simple are the velocity dispersions, which we
find through (12),
This set of power law components presented in the previous subsection
is very adequate to fit a broad class of models. It is however, not
fit to construct dynamical models for systems that have a central
density cusp. For such models, the isotropic distribution function
diverges in the limit
(see
e.g. Dehnen 1993; Hernquist 1990; Tremaine et al. 1994; Baes et al. 2005). This behaviour cannot be represented with the
current set of components. Therefore we extend the set of models
defined by the augmented density (34) to

(39) 
with two additional parameters
and s>0. This
sevenparameter family of dynamical components contains enough freedom
in the density profile to form a good family of base functions. If we
expand this expression as a power series in ,

(40) 
we see that we obtain a series of positive terms that all have the
form (34). Thus, our condition
is sufficient to
obtain a welldefined distribution function, which we can write down
immediately,

(41) 
or explicitly
Since the summations in j converge for every fixed value
of k, this distribution function is indeed welldefined.
Other dynamical properties, such as the moments of the distribution
function, can be derived in the same way. For the velocity dispersions
we find
with B_{x}(a,b) the incomplete Beta function.
5 Selfconsistent analytical models
In the previous section we have described a practical method to
construct analytical selfconsistent dynamical models that generate
any spherical density profile and with an anisotropy profile defined
by the general parameterized function (30). We are applying
this technique to construct dynamical models matching the density and
anisotropy profiles that have come out of cosmological simulations of
dark matter haloes, in order to gain insight into their phasespace
structure (Van Hese et al. 2007). In the remainder of this current paper we
will focus on simple analytical dynamical models, by investigating
whether the global sevenparameter family of dynamical models defined
by the augmented density (39) contains any simple
selfconsistent models (hence without the need to make a linear
combination). Thus we will look for potentialdensity pairs
that satisfy both Poisson's equation and the
relation

(45) 
with
the sevenparameter augmented density from
Eq. (39). The existence of such analytical dynamical
models would be a great step forward in our quest for simple but
realistic dynamical models that can be used as a framework in which to
initiate detailed numerical simulations.
One of the most obvious candidates is the Plummer model, as this model
has a rather straightforward potentialdensity pair
When we combine these functions with the expressions (39) and
(45), we obtain the condition

(48) 
A straightforward solution is obviously given by

(49) 

(50) 

(51) 
which yields the isotropic Plummer model, defined by
Our goal, however, is to determine the most general subspace of the
parameter space such that
the condition (48) is satisfied for all r. In particular
we aim for a subspace of the parameter space that has no restrictions
on
and
,
such that we obtain a family of
dynamical models with an arbitrary anisotropy at small and large
radii. It is obvious from Eq. (48) that, in order to
find nontrivial solutions, we have to set

(54) 

(55) 
We then obtain the expression

(56) 
These expressions are identical for all values of r if

(57) 

(58) 

(59) 
We now have constructed a general twoparameter family of
selfconsistent Plummer models with the augmented density profile

(60) 
From the condition
it follows that the central anisotropy
has to satisfy
.
This is in fact a special case of the
cusp slopecentral anisotropy theorem of An & Evans (2006b):
a model without a density cusp cannot have a radial velocity
anisotropy in the centre. The anisotropy at large radii can take any
value
.
By construction, the anisotropy profile of
this family of models reads

(61) 
which indeed tends towards
at small radii and towards
at large radii. The radial velocity dispersion profile
can be written as

(62) 
At large radii, the radial velocity dispersion profile shows a
r^{1} behaviour for all values of the parameters
and
,



(63) 
The asymptotic behaviour at small radii is



(64) 
Except for the models that are isotropic in the centre, the radial
velocity dispersions hence always tend to zero at small radii. The
asymptotic behaviour of the tangential velocity dispersions
follows immediately. At large radii
we obtain

(65) 
whereas at small radii

(66) 

Figure 1:
Top. Radial (red), tangential (green) and projected
(blue) velocity dispersion profiles for three Plummer models with
anisotropy values
and
displayed in the
figures. Bottom. The distribution function of these models,
represented by isoprobability contours in turning point space. High
values are indicated by red contours, low values by yellow
contours. In all plots, we have used normalized units with
G = M =
b = 1. 
Open with DEXTER 
In the top panels of Fig. 1 we plot the radial and
tangential velocity dispersions for a set of anisotropic Plummer
models. The three models in this figure all have the same
anisotropy
at large radii, but a
different central anisotropy. Apart from the radial and transversal
velocity dispersion, we also plot the projected velocity dispersion
on the plane of the sky, defined through

(67) 
Here,
is the projected mass density on the plane of the
sky,

(68) 
The projected velocity dispersion always reaches a finite value in the
centre and has an R^{1} dependence at large radii,

(69) 
The distribution function for our set of Plummer models can be written
as

(70) 
where the function
is a function defined as
(Dejonghe 1986)

(71) 
with
G^{m,n}_{p,q}(z) the Meijer Gfunction. This
function can conveniently be expressed as

(72) 
In the bottom panels of Fig. 1 we plot the
distribution function as a contour plot in the turning point
space for the same models as the upper panels. The change in
anisotropy from tangentially anisotropic at small radii to radially
anisotropic at large radii can easily be seen in the slope of these
contours.
An interesting subset of models in our twoparameter family of Plummer
models is the oneparameter family with .
These models are
isotropic in the inner regions and become anisotropic at large
radii. Their augmented density is given by

(73) 
and in the expression for the distribution function (70) only
the term corresponding to j=0 remains

(74) 
This subfamily of our current set of Plummer models was already
presented by Dejonghe (1987). Most of the kinematical
properties, including the projected properties such as dispersions and
higherorder moments of the line profiles, can be calculated
completely analytically.
Another extremely popular and simple potentialdensity pair is the
Hernquist model, defined by
Contrary to the Plummer model, this model has a central r^{1}density cusp and a more realistic r^{4} behaviour at large
radii. We can do the same experiment for the Hernquist model as we did
for the Plummer model. If we combine the potentialdensity pair
(75)(76) with expression (39)
and (45), we obtain

(77) 
Similarly as for the Plummer model, it is clear that we will only be
able to find a general nontrivial solution if we set

(78) 

(79) 
This yields the equation

(80) 
from which we find

(81) 

(82) 

(83) 

Figure 2:
Top. Radial (red), tangential (green) and projected
(blue) velocity dispersion profiles for three Hernquist models with
anisotropy values
and
displayed in the
figures. Bottom. The distribution function of these models,
represented by isoprobability contours in turning point space. High
values are indicated by red contours, low values by yellow
contours. In all plots, we have used normalized units with
G = M =
b = 1. 
Open with DEXTER 
We have now defined a twoparameter family of selfconsistent
Hernquist models with augmented density

(84) 
The parameter
can assume all values, whereas the
central anisotropy
is limited to
,
in agreement with the cusp slopecentral anisotropy theorem of
An & Evans (2006b). By construction, the anisotropy profile
of this family of Hernquist models reads

(85) 
which has the desired behaviour at small and large radii. The radial
dispersion profile reads

(86) 
At large radii, the radial velocity dispersion profile falls as r^{1},



(87) 
whereas the asymptotic behaviour at small radii is



(88) 
The radial velocity dispersion hence always disappears in the centre,
except for the models with the largest allowed central anisotropy
(
)
where it reaches a finite value. For the
asymptotic behaviour of the tangential velocity dispersions
at large radii we obtain

(89) 
whereas at small radii

(90) 
The projected velocity dispersion on the plane of the sky has a
R^{1} dependence at large radii,

(91) 
The distribution function can most conveniently be written as a series
of hypergeometric functions
if
,
and as
if
.
These sums only contain a finite number of terms if
is a positive integer number, i.e. when
is a positive integer or halfinteger
number. This particular subset of models, in which the outer regions
are always more tangentially anisotropic than the central regions, has
already been discussed by Baes & Dejonghe (2002).
In a similar manner as for the Plummer model, we plot in
Fig. 2 the velocity dispersions and the distribution
function for three Hernquist models with the same outer anisotropy
but different central anisotropy .
It is wellknown that the Plummer and the Hernquist potentialdensity
pairs can be generalized to a oneparameter family of models
characterized by
This potentialdensity pair was first described by
Veltmann (1979) and is a special subset (the
models) of the general set of potentialdensity pairs
considered by Zhao (1996). This potentialdensity pair
recently regained much interest because it supports dynamical models
that are hypervirial, i.e. in which the virial relation is not only
satisfied on a global but also on a local level
(Evans & An 2005; Iguchi et al. 2006; Sota et al. 2006). The parameter ,
lying in the range
,
determines the slope of the central density cusp. We
easily recognize the Plummer model with
as the only
coredensity member of the family and the Hernquist model as the model
with .
We can now repeat the same exercise as for the Plummer and Hernquist
models. After a little bit of algebra, we find that the parameters

(96) 

(97) 

(98) 

(99) 

(100) 
are the general solution for the condition of selfconsistency.
Notice that the initial condition
implies
,
which is in correspondence with the cusp slopecentral
anisotropy theorem (An & Evans 2006b). In other words, for
this family the condition
is also necessary to yield
physical models. In this manner we have constructed a threeparameter
family of dynamical models defined by the augmented density

(101) 
This augmented density selfconsistently supports the oneparameter
potentialdensity pair (94)(95) and has by
construction an anisotropy profile

(102) 
which varies smoothly from
in the centre towards
at large radii. The radial velocity dispersion can be
written as

(103) 
For general values of ,
the distribution function cannot be
simplified and should be taken as in formula (41) with the
values (96)(100). For rational values of
however, distribution function simplifies to a sum of
generalized hypergeometric functions. One such case is the model with
,
for which we obtain the potentialdensity pair
This model has a central density cups with a slope
,
which has been obtained by numerical simulations for dark
matter haloes in a CDM cosmological model (Moore et al. 1998).
6 Discussion and conclusions
In this paper, we have described a method to construct selfconsistent
spherical dynamical models that have an arbitrary density and an arbitrary
anisotropy profile. While the general formulation is not so
useful since the inversion formulae are rather complicated and
numerically unstable, we propose an approach with a parameterized
anisotropy profile. We put forward a very general fourparameter
anisotropy profile in which the four parameters control the anisotropy
at the centre and at large radii, the typical transition radius
between the two regimes and the sharpness in which this transition
takes place. This parameterized function can therefore represent a
large number of monotonically varying anisotropy profiles.
Based on these anisotropy profiles, we present a sevenparameter set
of dynamical components that all have the same anisotropy. For each of
these components the augmented density
is a
sufficiently simple function so that the corresponding distribution
function can be written analytically (as a double series). For many
reasonable choices of the parameters, the distribution function can be
written as a sum of generalized hypergeometric series. The lowerorder
moments can be expressed by simple analytical functions  for
example, the velocity dispersions of all components can be expressed
as incomplete Beta functions.
Once the anisotropy profile has been fixed, the construction of a
dynamical model with a given potentialdensity pair consists of a
determination of the best linear combination of the various
components. This can be achieved with a quadratic programming
approach, a technique already in use by various dynamical modellers
(Kuijken & Merrifield 1993; Merrifield & Kuijken 1994; Dejonghe 1989; Gerhard et al. 1998). This dynamical modelling technique is rather
advanced: dynamical models can be constructed by fitting a variety of
data, such as any moments of the distribution function, the projected
moments along the line of sight, the full line profiles or even
complete spectra. One could argue what the novelty of our approach is
compared to the general powerful techniques already presented by these
authors. The use of the quadratic programming technique is very
general and does not require all components to have the same
anisotropy profile. Any possible combination of components can in
principle be combined. The main strategy usually consists in choosing
the components in such a way that the distribution functions are
simple enough; one can for example choose only Fricke components which
are double power laws both in augmented density and distribution
function. As long as the data base of possible models is large enough,
such components allow to fit basically any dynamical model since a
linear combination of two components with a different density profile
and different constant anisotropy results in a model with a
nonconstant anisotropy. However, practice has learned that it is
generally difficult to construct models with a strong radial
anisotropy at large radii. The reason for this problem is that one
needs to populate the model with radial orbits that reach large
radii. Since such orbits also contribute to the density at small
radii, it requires a delicate finetuning of the different components
to both satisfy the density and anisotropy constraints at small radii
while still retaining the radial anisotropy at large radii. In our
present case, this problem does not surface since we only have to fit
the function ,
or equivalently the density profile. Any
combination of the base components will automatically have the correct
anisotropy profile. This was actually one of the motivations for the
construction of this set of dynamical components.
Apart from presenting an approach to construct detailed dynamical
models with a given density and anisotropy profile, our study has also
lead to a set of simple onecomponent dynamical models. The models we
have generated are selfconsistent realizations of the Plummer and
Hernquist models. We have also generalized these models to a complete
set of models that selfconsistently generate the family of Veltmann
(or generalized Plummer) potentialdensity pairs. Contrary to most
analytical dynamical models that have appeared in the literature, the
anisotropy profile of these new dynamical models can be chosen
arbitrarily in the sense that we can choose the anisotropy at both
small and large radii. Both numerical simulations and observational
evidence clearly indicates that the outer regions of galaxies and dark
matter haloes are typically mildly to significantly radially
anisotropic and that the inner regions can be significantly
nonisotropic. Realistic simulations should use this information and
not impose unrealistic constraints on the anisotropy. We therefore
believe that these new analytical models are an important step forward
compared to the isotropic or OsipkovMerritt type dynamical models
that are often used. We encourage numerical astrophysicists to use
these new distribution functions to generate the initial conditions
for their simulations of galaxies or dark matter haloes. Numerical
implementation of the most important dynamical properties of these
models can be obtained from the authors.
We are currently employing the modelling technique described in this
paper to construct dynamical models for dark matter haloes
(Van Hese et al. 2007). We are using constraints on the density and the
anisotropy profile of these haloes from detailed Nbody
simulations. This will be useful to investigate the phase space
structure of dynamical models and will also provide the community with
an advanced toy model for dark haloes from which more simulations can
be initiated. This nicely illustrates that the rise of numerical
simulations does not erase the need for analytical work, but rather
that there is a symbiosis between both approaches.
 An, J. H., &
Evans, N. W. 2006a, AJ, 131, 782 [NASA ADS] [CrossRef]
 An, J. H., &
Evans, N. W. 2006b, ApJ, 642, 752 [NASA ADS] [CrossRef] (In the text)
 Baes, M., &
Dejonghe, H. 2002, A&A, 393, 485 [NASA ADS] [CrossRef] [EDP Sciences]
 Baes, M., &
Dejonghe, H. 2004, MNRAS, 351, 18 [NASA ADS]
 Baes, M., Dejonghe,
H., & Buyle, P. 2005, A&A, 432, 411 [NASA ADS] [CrossRef] [EDP Sciences]
 Buyle, P., Hunter, C.,
& Dejonghe, H. 2007, MNRAS, 375, 773 [NASA ADS]
 Carlberg, R. G.,
et al. 1997, ApJ, 485, L13 [NASA ADS] [CrossRef]
 Cole, S., & Lacey,
C. 1996, MNRAS, 281, 716 [NASA ADS]
 Colín, P.,
Klypin, A. A., & Kravtsov, A. V. 2000, ApJ, 539,
561 [NASA ADS] [CrossRef]
 Cuddeford, P. 1991,
MNRAS, 253, 414 [NASA ADS] (In the text)
 Cuddeford, P., &
Louis, P. 1995, MNRAS, 275, 1017 [NASA ADS]
 Dehnen, W. 1993,
MNRAS, 265, 250 [NASA ADS]
 Dehnen, W., &
McLaughlin, D. E. 2005, MNRAS, 363, 1057 [NASA ADS] (In the text)
 Dejonghe, H. 1986,
Phys. Rep., 133, 217 [NASA ADS] [CrossRef]
 Dejonghe, H. 1987,
MNRAS, 224, 13 [NASA ADS] (In the text)
 Dejonghe, H. 1989,
ApJ, 343, 113 [NASA ADS] [CrossRef]
 Diemand, J., Moore,
B., & Stadel, J. 2004, MNRAS, 352, 535 [NASA ADS]
 Evans, N. W.,
& An, J. 2005, MNRAS, 360, 492 [NASA ADS]
 Fricke, W. 1952,
Astron. Nach., 280, 193 [NASA ADS]
 Gebhardt, K.,
et al. 2003, ApJ, 583, 92 [NASA ADS] [CrossRef] (In the text)
 Gerhard, O., Jeske,
G., Saglia, R. P., & Bender, R. 1998, MNRAS, 295, 197 [NASA ADS]
 Gradshteyn,
I. S., & Ryzhik, I. M. 1965, Table of integrals,
series and products, 4th edn., ed. Yu. V. Geronimus, & M. Yu.
Tseytlin (New York: Academic Press)
(In the text)
 Fukushige, T., &
Makino, J. 2001, ApJ, 557, 533 [NASA ADS] [CrossRef]
 Hansen, S. H.,
& Moore, B. 2006, New Astron., 11, 333 [NASA ADS] (In the text)
 Hansen, S. H.,
& Stadel, J. 2006, J. Cosmol. AstroPart. Phys., 5, 14 [NASA ADS] (In the text)
 Hénon, M. 1960,
Ann. Astrophys., 23, 474 [NASA ADS]
 Hénon, M. 1973,
A&A, 24, 229 [NASA ADS]
 Hernquist, L. 1990,
ApJ, 356, 359 [NASA ADS] [CrossRef]
 Hiotelis, N. 1994,
A&A, 283, 783 [NASA ADS]
 Iguchi, O., Sota, Y.,
Nakamichi, A., & Morikawa, M. 2006, Phys. Rev. E, 73,
046112 [NASA ADS] [CrossRef]
 Kuijken, K., &
Merrifield, M. R. 1993, MNRAS, 264, 712 [NASA ADS]
 Jaffe, W. 1983, MNRAS,
202, 995 [NASA ADS]
 Kazantzidis, S.,
Magorrian, J., & Moore, B. 2004, ApJ, 601, 37 [NASA ADS] [CrossRef] (In the text)
 Kronawitter, A.,
Saglia, R. P., Gerhard, O., & Bender, R. 2000, A&AS,
144, 53 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Louis, P. D.
1993, MNRAS, 261, 283 [NASA ADS]
 Mathai, A. M. 1993, A
Handbook of Generalized Special Functions for Statistical and
Physical Sciences (Oxford: Oxford University Press)
(In the text)
 Merrifield,
M. R., & Kuijken, K. 1994, ApJ, 432, 575 [NASA ADS] [CrossRef]
 Merritt, D. 1985, AJ,
90, 1027 [NASA ADS] [CrossRef] (In the text)
 Merritt, D. 1985,
MNRAS, 214, 25P [NASA ADS]
 Moore, B., Governato,
F., Quinn, T., Stadel, J., & Lake, G. 1998, ApJ, 499, L5 [NASA ADS] [CrossRef] (In the text)
 Oñorbe, J.,
DomínguezTenreiro, R., Sáiz, A., & Serna, A. 2007,
MNRAS, 376, 39 [NASA ADS] (In the text)
 Osipkov, L. P.
1979, Pis'ma Astronomicheskii Zh., 5, 77 [NASA ADS] (In the text)
 Plummer, H. C.
1911, MNRAS, 71, 460 [NASA ADS]
 Quinlan, G. D.,
& Hernquist, L. 1997, New Astron., 2, 533 [NASA ADS] (In the text)
 Quinlan, G. D.,
Hernquist, L., & Sigurdsson, S. 1995, ApJ, 440, 554 [NASA ADS] [CrossRef] (In the text)
 Sota, Y., Iguchi, O.,
Morikawa, M., & Nakamichi, A. 2006, Progr. Theor. Phys. Suppl.,
162, 62 [NASA ADS]
 Tremaine, S.,
Richstone, D. O., Byun, Y.I., et al. 1994, AJ, 107,
634 [NASA ADS] [CrossRef]
 Van Hese, E., Baes, M., &
Dejonghe, H. 2007, in preparation
(In the text)
 Veltmann,
U. I. K. 1979, AZh, 56, 976 [NASA ADS] (In the text)
 Zhao, H. 1996, MNRAS,
278, 488 [NASA ADS]
Online Material
Appendix A: Derivation of the distribution function
We now present the detailed calculation of the distribution function
corresponding with the augmented density (34). We have to
consider several cases, depending on the parameter
.
First, we consider the special case where
is zero or a
natural number. Then the augmented density is simply a finite sum of
positive Fricke components
(Hénon 1973; Fricke 1952)

(A.1) 
which leads immediately to Eq. (36).
When
we can apply the LaplaceMellin formalism: the
connection between the distribution function and the augmented density
is then (Dejonghe 1986)

(A.2) 
Since the augmented density is a separable function of
and r,
the transforms can be calculated separately:

(A.3) 
and

(A.4) 
where
lies in the convergence strip
.
Thus

(A.5) 
The inversion of the Laplace transform is straightforward,

(A.6) 
leaving us with the inversion of the Mellin transform
Taking into account the definition of the Fox Hfunction,

(A.8) 
the distribution function can be written as

(A.9) 
or equivalently

(A.10) 
The integration can be performed along three possible paths C, which are
equivalent: if the integral converges for more than one of these three
paths, then the result is the same. If the integral converges for
only one path, then that is the only one to be considered. The valid
contours are
 C is a path C_{1} from
to ,
such that the
poles of
lie on one
side and the poles of
lie on the
other side. The convergence is absolute if
or if
and
.
This condition
is found by investigating the more familiar criterion for Meijer
Gfunctions (see Appendix B). In addition, if
and
the integral is semiconvergent for
(Dejonghe 1986, Appendix).
 C is a loop C_{2}, starting and ending at ,
that
encircles the poles of
once in the positive
direction and none of the poles of
.
The integral
then converges if
.
 C is a loop C_{3}, starting and ending at ,
that
encircles the poles of
once in the
negative direction and none of the poles of
.
The integral then
converges if
.
It can easily be seen that the convergence criterion for path C_{1}ensures a welldefined distribution function, i.e. continuous and
nonnegative. Indeed, in this case the contour is a line parallel to
the imaginary axis from
to
with
.
On this path the real part of the Gamma
functions is positive (where the condition
ensures that
by choosing s_{0} sufficiently close
to
). Hence, the distribution function is nonnegative
everywhere.
In the case of
and not a natural number, the Mellin
transform does not exist. However, we can solve this problem by
rewriting the augmented density in a similar way as
Eq. (A.1): defining n as the smallest natural number such
that
,
we obtain

(A.11) 
For the individual terms the LaplaceMellin formalism does apply,
and the distribution function becomes

(A.12) 
Now, it is easily observed that we can choose the integration paths
C_{(k)} to be identical, as the contours C_{1}, C_{2} or C_{3}defined above. Therefore the summation can be performed inside the
integral. Using the properties

(A.13) 
with
the Pochhammer symbol, we find

(A.14) 
Finally, with the identity

(A.15) 
the equation for the distribution function also reduces to
Eq. (A.9).
If the convergence criterion for path C_{1} is valid, then every
contour C_{(k)} can be taken as a line
to
with
.
Again,
for each k the real part of the integrand is positive, so that the
distribution function is welldefined.
Appendix B: A practical series expansion for the distribution function
We now seek a more practical form for the distribution function. To
this aim we first consider the special case where
is a
rational number, denoting
.
Then we can write

(B.1) 
Now, using the multiplication formula for the Gamma function

(B.2) 
we can write the integral in the form of a Meijer Gfunction:

(B.3) 
with
= 
(B.4) 
Since for a general Meijer Gfunction

(B.5) 
with p=q and
real, the convergence criterion
for path C_{1} states that (Mathai 1993)

(B.6) 
or
= 0, 
(B.7) 
< 1, 
(B.8) 
we indeed obtain the conditions
or
and
,
as
used in Appendix A.
This Meijer Gfunction can be calculated as a sum of generalized
hypergeometric functions (Gradshteyn & Ryzhik 1965, Eqs. (9.303) and (9.304)). For z<1 the following equation is
valid:
where the prime by the product symbol denotes the omission
of the product when i=l, and the asterisk in the hypergeometric
function indicates the omission on the i th parameter.
Analogously, the equation for z>1 reads
With the coefficients from Eq. (B.4), we obtain for
Now, with the aid of the identities

(B.12) 
and

(B.13) 
we can simplify this expression to
so that, using again Eq. (B.2), the distribution function
reduces to
Finally, the double summation can be grouped into a single index k=i+nj,
and we obtain for

(B.16) 
Similarly, for
,
we find

(B.17) 
Although these expressions have been derived for rational values of
,
they can be generalized to any real value, since these
functions are continuous in .
Hence we indeed obtain
Eq. (36).
Copyright ESO 2007