Search
Close
Search
Advanced Search
Search Menu
AI Discovery Assistant

Abstract

We present an approach to the design of distribution functions that depend on the phase-space coordinates through the action integrals. The approach makes it easy to construct a dynamical model of a given stellar component. We illustrate the approach by deriving distribution functions that self-consistently generate several popular stellar systems, including the Hernquist, Jaffe, and Navarro, Frenk and White models. We focus on non-rotating spherical systems, but extension to flattened and rotating systems is trivial. Our distribution functions are easily added to each other and to previously published distribution functions for discs to create self-consistent multicomponent galaxies. The models this approach makes possible should prove valuable both for the interpretation of observational data and for exploring the non-equilibrium dynamics of galaxies via N-body simulations.

galaxies: kinematics and dynamics, galaxies: structure, dark matter

1 INTRODUCTION

Axisymmetric equilibrium models are extremely useful tools for the study of galaxies. A real galaxy will never be in perfect dynamical equilibrium – it might be accreting dwarf satellites, or being tidally disturbed by the gravitational field of the group or cluster to which it belongs, or displaying spiral structure – but an axisymmetric equilibrium model will usually provide a useful basis from which a more realistic model can be constructed by perturbation theory.

By Jeans (1915) theorem, every equilibrium model can be described by a distribution function (DF) that depends on the phase-space coordinates |$({\boldsymbol x},{\boldsymbol v})$| only through isolating integrals of motion. In an axisymmetric potential, most orbits prove to be quasi-periodic, with the consequence that they admit three isolating integrals (Arnold 1978). Consequently, a generic DF for an axisymmetric equilibrium galaxy is a function of three variables.

The major obstacle to exploiting this insight is that we have analytic expressions for only two isolating integrals of motion in a general axisymmetric potential, namely the energy |$E={\textstyle {1\over 2}}v^2+\Phi ({\boldsymbol x})$| and the component of the angular momentum about the symmetry axis, |$J_\phi =({\boldsymbol x}\times {\boldsymbol v})_z$|⁠. Several authors have examined model galaxies with DFs of the two-integral form f(E, Jϕ) (Prendergast & Tomer 1970; Wilson 1975; Rowley 1988; Evans 1994), but in such models the velocity dispersions σR and σz in the radial and vertical directions are inevitably equal. This condition is seriously violated in our Galaxy and we have no reason to suppose that the condition is better satisfied in any external galaxy. Hence, it is mandatory to extend the DF's argument list to include a ‘non-classical’ integral, I3, for which we do not have a convenient expression.

Since any function J(E, Jϕ, I3) of three isolating integrals is itself an isolating integral, we actually have an enormous amount of freedom as to what integrals to use as arguments of the DF. Given that we must use at least one integral for which we lack an expression for its dependence on |$({\boldsymbol x},{\boldsymbol v})$|⁠, there is a powerful case for making the DF's arguments action integrals. These integrals are alone capable as serving as the three momenta Ji of a canonical coordinate system – this property makes them the bedrock of perturbation theory. Their canonically conjugate variables, the angles θi, have two remarkable properties: (i) along any orbit they increase linearly with time at rates |$\Omega _i({\boldsymbol J})$|⁠, so
(1)
and (ii) they make the ordinary phase-space coordinates periodic functions
(2)
The actions Ji also have nice properties. In particular, (i) any triple of finite numbers (Jr, Jϕ, Jz) with Jr, Jz ≥ 0 corresponds to a bound orbit with the orbit |${\boldsymbol J}=0$| being that on which a star is stationary at the middle of the galaxy, and (ii) the volume of phase space occupied by orbits with actions in |${\rm d}^3{\boldsymbol J}$| is |$(2\pi )^3{\rm d}^3{\boldsymbol J}$|⁠. Consequently, any non-negative function |$f({\boldsymbol J})$| that tends to zero as |$|{\boldsymbol J}|\rightarrow \infty$| and has a finite integral |$\int {\rm d}^3{\boldsymbol J}\,f({\boldsymbol J})$| specifies a valid galaxy model of mass
(3)
The actions are defined by integrals
(4)
where γi is a closed path in phase space. If we require that the first action Jr quantifies the extent of a star's radial excursions and the third action Jz quantifies the extent of its excursions either side of the potential's equatorial plane, then the actions are unambiguously defined. What we here call Jr is sometimes called JR or Ju, and what we call Jz is sometimes called Jϑ or Jv, but no significance attaches to these different notations. In a spherical potential Jz = L − |Jϕ|, where L is the magnitude of the angular momentum vector.

To obtain the observable properties of a model defined by |$f({\boldsymbol J})$|⁠, for example its density distribution |$\rho ({\boldsymbol x})=\int {\rm d}^3{\boldsymbol v}\,f({\boldsymbol J})$| and its velocity dispersion tensor |$\sigma ^2_{ij}({\boldsymbol x})$|⁠, one has to be able to evaluate |${\boldsymbol J}({\boldsymbol x},{\boldsymbol v})$| in an arbitrary gravitational potential. Recently, a number of techniques have been developed for doing this (Binney 2012a; Sanders & Binney 2014, 2015). Consequently, while the last word on action evaluation has likely not yet been written, we now have algorithms that enable one to extract the observables from a DF |$f({\boldsymbol J})$| with reasonable accuracy.

DFs |$f({\boldsymbol J})$| that depend on the phase-space coordinates only through the actions were first used to model the disc of our Galaxy in an assumed gravitational potential (Binney 2010, 2012b). Recently, Binney (2014, hereafter B14) showed how to derive the self-consistent gravitational potential that is implied by a given |$f({\boldsymbol J})$| by exploring a family of flattened, rotating models that he derived from the ‘ergodic’ DF of the isochrone model: that is the DF f(H) that depends on the phase-space coordinates only through the Hamiltonian |$H=\frac{1}{2} v^2+\Phi ({\boldsymbol x})$|⁠. Hénon (1960) derived the isochrone's ergodic DF, and in the case of the isochrone potential explicit expressions are available for |${\boldsymbol J}({\boldsymbol x},{\boldsymbol v})$| and |$H({\boldsymbol J})$| (Gerhard & Saha 1991). Substituting |$H({\boldsymbol J})$| in f(H) B14 obtained the DF |$f({\boldsymbol J})$| of the isotropic isochrone model. In this paper, we present simple analytic functions |$f({\boldsymbol J})$| that generate nearly isotropic models of other widely used models, such as the Hernquist (1990), Jaffe (1983), and Navarro, Frenk & White (1996, hereafter NFW) models.

Once a DF of the form |$f({\boldsymbol J})$| is available for a spherical, non-rotating model, the procedure B14 used to flatten the isochrone sphere and to set it rotating can be used to flatten and/or set rotating one's chosen model. So DFs for spherical models in the form |$f({\boldsymbol J})$| are valuable starting points from which quite general axisymmetric models are readily constructed.

Galaxies are generally considered to consist of a number of components, such as a disc, a bulge, and a dark halo, that cohabit a single gravitational potential. If we represent each component by a DF of the form |$f({\boldsymbol J})$|⁠, it is straightforward to find the gravitational potential in which they are all in equilibrium (e.g. Piffl et al. 2014; Piffl, Penoyre & Binney in preparation). An analogous composition using DFs of the form f(E, Jϕ, I3) has never been achieved and may be impossible, because when components are added, their potentials must be added, and the energies of physically similar orbits in a given component are quite different before and after we add in the potential of another component. For example, the orbit on which a star sits at the centre of the galaxy will have different energies before and after addition. If E is used as an argument of the DF, the change in E will change the density of stars on the given orbit, which is contrary to the fundamental idea of building up the galaxy by adding components. By contrast, the actions of the orbit on which a star sits at the galactic centre vanish in any potential, and if a component is defined by |$f({\boldsymbol J})$|⁠, it contributes the same density of stars to this orbit regardless of the external potential in which that component finds itself. This fact is a major motivation for discovering what DF of the form |$f({\boldsymbol J})$| is required to generate each component of a galaxy.

The DF of an isotropic spherical model must depend on the actions only via the Hamiltonian |$H({\boldsymbol J})$|⁠. The dependence of f on H is readily obtained from the inversion formula of Eddington (1916), but an exact expression for |$H({\boldsymbol J})$| is only available for the isochrone potential and its limiting cases, the harmonic oscillator and Kepler potentials. Our ignorance of |$H({\boldsymbol J})$| for potentials other than the isochrone amounts to a barrier to the extension of B14's approach to model building. One way to break through this barrier is to devise numerical approximations to |$H({\boldsymbol J})$| and some success has been had in this direction by Fermani (2013) and Williams, Evans & Bowden (2014). In this paper, we pursue a slightly different strategy, which is to develop simple algebraic expressions for DFs |$f({\boldsymbol J})$| that generate self-consistent models that closely resemble popular spherical systems. We also show that a very simple form of |$f({\boldsymbol J})$| generates a model that is almost identical to the isochrone sphere and we give a useful analytic expression for the radial action as a function of energy and angular momentum for a Hernquist sphere.

The paper is organized as follows. In Section 2, we use analytic arguments to infer |$f({\boldsymbol J})$| for scale-free models. These models are not physically realizable as they stand, so in Section 3 we consider models that consist of two power-law sections joined at a break radius. In Section 4, we extract realizable models from scale-free models by the alternative strategy of adding a core to the system and/or tidally truncating the model. Section 5 sums up.

2 POWER-LAW MODELS

Consider a gravitational potential that scales as a power of the distance from the galactic centre, i.e. |$\Phi (\xi {\boldsymbol x})\propto \xi ^a\Phi ({\boldsymbol x})$| with |$a\not= 0$|⁠: in the limit a → 0 the gravitational potential tends to a logarithmic potential, which is an interesting special case that we will treat in Section 2.1.

An orbit in a power-law potential has time-averaged kinetic and potential energies, K and W, respectively, that are related by the virial theorem: 2K = aW. The instantaneous total energy, given by the sum of the instantaneous kinetic and potential energies, is conserved along the orbit and consequently is given by
(5)
In any power-law potential, we need only to study orbits of one arbitrarily chosen energy E because each of these orbits can be rescaled to a similar orbit at any given energy E′. Indeed, if an orbit is rescaled by a spatial factor, i.e. |${\boldsymbol x}\rightarrow {\boldsymbol x}^{\prime }=\xi {\boldsymbol x}$|⁠, then the orbit's total energy scales as
(6)
since obviously W → W′ = ξaW. Further |$v^2\propto K={\textstyle {1\over 2}}aW$|⁠, so under rescaling |${\boldsymbol v} \rightarrow {\boldsymbol v}^{\prime }=\xi ^{a/2}{\boldsymbol v}$|⁠.
Given the scalings derived above for |${\boldsymbol x}$| and |${\boldsymbol v}$| it follows that
(7)
Thus, both the energy and the actions of an orbit that is rescaled by the spatial factor ξ are rescaled by powers of this factor.
From equations (6) and (7) we deduce that the Hamiltonian is of the form
(8)
where |$h({\boldsymbol J})$| is a homogeneous function of degree one, i.e. |$h(\zeta {\boldsymbol J})=\zeta h({\boldsymbol J})$| for every constant ζ. In particular, H is itself a homogeneous function of the three actions of degree a/(1 + a/2). It is easy to check that equation (8) gives the correct scalings |$H\propto |{\boldsymbol J}|$| and |$H\propto |{\boldsymbol J}|^{-2}$| for the harmonic oscillator (a = 2) and Kepler (a = −1) potentials. Williams et al. (2014) derive a closely related result in which a specific form is proposed for |$h({\boldsymbol J})$|⁠.
The homogeneous function h is strongly constrained by the orbital frequencies. Indeed
(9)
In a scale-free model, the frequency ratio on the left is a homogeneous function of degree zero, i.e. scale independent, in agreement with the right-hand side. A natural choice for h that we will use extensively is
(10)
In a scale-free model, this is homogeneous of degree one, as required. Moreover, so long as the frequency ratios do not change rapidly within a surface of constant energy in action space, the derivatives of h satisfy equation (9) to good precision.

In the definition (10) of |$h({\boldsymbol J})$|⁠, the modulus of the angular momentum Jϕ appears because we are concerned with the construction of the part of the DF that is even in Jϕ. If we wish to set the model rotating, we will add to this even part an odd part as discussed by B14.

Consider now the density distribution that generates a power-law potential. In the spherical case,1 we have
(12)
Hence
(13)
If a = −1, we recover the expected result ρ = 0, but for a ≠ 0 we obtain the polytropic relation for index n = 1 − 2/a (e.g. Binney & Tremaine 2008, section 4.3.3a):
(14)
From this relation, it is easy to derive the ergodic DF
(15)
from Eddington's formula (e.g. Evans 1994). From equations (8) and (15) it follows that the DF of a power-law model is
(16)
The DF of a power-law model is itself a power law of the three actions and the exponent is completely determined by that of |$\Phi ({\boldsymbol x})$|⁠.

2.1 Logarithmic potentials

Now consider the limit a → 0 when the scaling of Φ becomes additive
(17)
where vc is a constant that one can easily show is the circular speed. Since galaxies have quite flat circular speed curves, potentials of this form are very useful.
The kinetic energy K does not change on rescaling, while the potential energy |$W \rightarrow W^{\prime }=W+v_{\rm c}^2\log (\xi )$|⁠, so
(18)
The invariance of K implies invariance of |${\boldsymbol v}$| under orbit rescaling, so the scaling of the actions is
(19)
We now use each side of this equation as the argument of a homogeneous function of degree one, |$h({\boldsymbol J})$|⁠, and obtain
(20)
or on rearrangement
(21)
Here E′ and E are the energies of any two orbits whose actions |${\boldsymbol J}^{\prime }$| and |${\boldsymbol J}$| are proportional to each other. We can choose to make |${\boldsymbol J}$| an orbit with vanishing energy, and we can choose h to be the homogeneous function that satisfies |$h({\boldsymbol J})=1$| as |${\boldsymbol J}$| moves over the surface E = 0 in action space. With these choices, we have
(22)
The ergodic DF that self-consistently generates the spherical logarithmic potential is well known to be
(23)
where |$\sigma ^2=v_{\rm c}^2/2$| and E0 is a constant (e.g. Binney & Tremaine 2008, section 4.3.3b). Using equation (22) it follows that the ergodic DF is
(24)
This result is consistent with the limit a → 0 of equation (16) for a power-law model.

Note that equation (24) implies that the phase-space density diverges as |${\boldsymbol J}\rightarrow 0$|⁠. It follows that this DF unambiguously specifies the singular isothermal sphere, in contrast to the DF (23), from which one can derive both cored and singular isothermal spheres (e.g. Binney & Tremaine 2008, section 4.3.3b). It is characteristic of DFs of the form |$f({\boldsymbol J})$| that they uniquely and transparently specify the phase-space density both at the centre of the model (⁠|${\boldsymbol J}=0$|⁠) and for marginally bound orbits |$({\boldsymbol J}\rightarrow \infty )$|⁠. From a DF that depends on energy, by contrast, the phase-space density at the centre of the model is implicitly specified by the boundary condition adopted at r = 0 when solving Poisson's equation for the self-consistent potential.

The considerations of the last paragraph apply equally to the power-law DFs (16): although we used the standard form (15) of the energy-based DF of the polytropes to derive this DF, it implies infinite phase-space density at the system's centre, so it is inconsistent with familiar cored polytropes, such as the Plummer model.

3 TWO-POWER MODELS

Any power-law model is problematic in the sense that the mass interior to radius r diverges as r → ∞ if the density declines as r−b with b ≤ 3, and the mass outside radius r diverges as r → 0 when b ≥ 3. Hence, there is no value of b for which the model is physically reasonable at both large and small r. One way we can address this problem is to assume that ρ scales as different powers of radius at small and large radii. A widely used family of models of this type is given by the density profile
(25)
where rb is the break radius (e.g. Binney & Tremaine 2008). Three particular cases of importance are the Jaffe (1983) model (α, β) = (2, 4), the Hernquist (1990) model (α, β) = (1, 4), which belong to the family of Dehnen (1993) models (β = 4), and the NFW model (α, β) = (1, 3) (Navarro et al. 1996). The ergodic DFs of the Jaffe and Hernquist models are known analytic function, but that of the NFW model is not. Our goal in this section is to find analytic functions |$f({\boldsymbol J})$| that generate models that closely resemble these three classic models.
In the regime r ≪ rb the mass M(r) enclosed by the sphere of radius r is M ∝ r3−α, so the gravitational acceleration is dΦ/dr ∝ r1−α and thus the potential drop between radius r and the centre is
(26)
Setting a = 2 − α we can now employ the results we derived above for power-law potentials to conclude that
(27)
The Hernquist and NFW models both have α = 1 so we expect their DFs to have asymptotic behaviour
(28)
A Jaffe model has α = 2, so the asymptotic behaviour of the Jaffe model's DF as |${\boldsymbol J}\rightarrow 0$| is given by equation (24).
Consider now the asymptotic behaviour of a two-power model as r → ∞. If the model has finite mass, the potential will asymptote to the Kepler potential, Φ ∝ r−1, so ρ ∝ |Φ|β. In the Kepler regime, the dependence of the Hamiltonian on the actions is (e.g. Binney & Tremaine 2008, equation 3.226a)
(29)
where |$g({\boldsymbol J})$| is a homogeneous function of degree one. Although ρ is a simple power of |Φ| we cannot employ the polytropic formula (15), because that rests on Poisson's equation, which does not apply in this case: the model's envelope is a collection of test particles that move in the Kepler potential generated by its core. We instead go back to Eddington's formula
(30)
where |${\cal E}=-E$| and Ψ = −Φ. From this formula, it is easy to show that ρ ∝ Ψβ implies
(31)
Combining this with equation (29) we conclude that for β > 3 the asymptotic behaviour of a double-power DF is
(32)
For the Jaffe and Hernquist models β = 4, so for these models
(33)
Now that we have the asymptotic behaviour of f in the limits of both small and large |${\boldsymbol J}$|⁠, it is straightforward to devise a suitable form of the DF
(34)
Here M0 is a constant that has the dimensions of a mass and J0 is a characteristic action. If the two homogeneous functions are normalized such that |$h({\boldsymbol J})\simeq g({\boldsymbol J})\simeq |{\boldsymbol J}|$|⁠, orbits that linger near the break radius rb have |$|{\boldsymbol J}|\simeq J_0$|⁠. These conditions ensure that f tends to the required powers of h and g when |$|{\boldsymbol J}|\ll J_0$| and |$|{\boldsymbol J}|\gg J_0$|⁠, respectively.
We use different homogeneous functions for the regimes of small and large |${\boldsymbol J}$| because the frequency ratios in these two regimes will differ. In the Kepler regime, which is handled by g, all frequencies are equal, so if we require an isotropic model we choose
(35)
In the regime of small |${\boldsymbol J}$|⁠, Ωr > Ωϕ = Ωz, and we take h to be of the form (10) with a frequency ratio that is less than unity. Unfortunately, in this regime the frequency ratio does vary over a surface of constant energy and an exactly isotropic model cannot be constructed using constant ratios. We simply use Ωϕ/Ωr = Ωz/Ωr = 1/2, which are the frequency ratios of a harmonic oscillator.

The DF (34) is infinite on the orbit |${\boldsymbol J} = 0$| of a star that is stationary at the model's centre. Cuspy models such as the Hernquist, Jaffe and NFW models do have such centrally divergent DFs, while in other cored systems the phase-space density reaches a finite maximum. Cored systems will be treated in Section 4.

3.1 Technicalities

Here, we touch on some technical issues that arise when one sets out to recover the observable properties of a model from the DF that defines it. The first step is to normalize the DF to the desired total mass by evaluating the integral (3). When the DF depends only on the function |$h({\boldsymbol J})$| defined by equation (10) [i.e. the case |$g({\boldsymbol J})=h({\boldsymbol J})$|] with the frequency ratios ω ≡ Ωϕ/Ωr = Ωz/Ωr taken to be constant, it is convenient to change coordinates from (Jr, Jϕ, Jz) to (Jr, L, Jz) and integrate out Jz, and then to change coordinates to (h, L) and integrate out L. Then one finds
(36)
In the more general case, when |$h({\boldsymbol J})\ne g({\boldsymbol J})$|⁠, the integral (3) cannot be reduced to one dimension. Equation (36) can be written as
(37)
where |${\boldsymbol y}\equiv {\boldsymbol J}/J_0$|⁠. The integral in equation (37) is dimensionless and depends only on the model's parameters α, β and on the forms of the homogeneous functions h and g. It can therefore be computed at the outset. Then the value of M0 can be set that ensures that the model has whatever mass is required.
The physical scales of the models are determined by the action scale J0 and by the mass scale M0, so the natural length-scale is
(38)
In following sections, we will present |$f({\boldsymbol J})$| analogues of three classic models that have a finite mass: the Hernquist, Jaffe and isochrone models. For our analogue models the top row of Table 1 gives the ratio rh/r0 of half-mass radius to the scale radius defined by equation (38). The second row gives for the classical models the ratio of rh to the break radius, and we see that for the Hernquist model r0 = rb to good precision, while in the other two cases the difference between r0 and rb is less than 25 per cent.
Table 1.
Open in new tab

The ratio of the half-mass radius rh to the scale radius r0, defined by equation (38), for the |$f({\boldsymbol J})$| isochrone, |$f({\boldsymbol J})$| Hernquist and |$f({\boldsymbol J})$| Jaffe models. For comparison, we list also the ratio rh/rb, where rb is the break radius, of the corresponding classical models.

IsochroneHernquistJaffe
rh/r03.42.420.76
rh/rb3.062.411
IsochroneHernquistJaffe
rh/r03.42.420.76
rh/rb3.062.411
Table 1.
Open in new tab

The ratio of the half-mass radius rh to the scale radius r0, defined by equation (38), for the |$f({\boldsymbol J})$| isochrone, |$f({\boldsymbol J})$| Hernquist and |$f({\boldsymbol J})$| Jaffe models. For comparison, we list also the ratio rh/rb, where rb is the break radius, of the corresponding classical models.

IsochroneHernquistJaffe
rh/r03.42.420.76
rh/rb3.062.411
IsochroneHernquistJaffe
rh/r03.42.420.76
rh/rb3.062.411

Once |$f({\boldsymbol J})$| has been normalized, we are able to determine the potential |$\Phi ({\boldsymbol x})$| that the model self-consistently generates by the iterative procedure described by B14.

3.2 Worked examples

3.2.1 The Hernquist model

The Hernquist (1990) model is an interesting example both because it is a widely used model and because we can derive its ergodic DF as a function of the actions for comparison with the |$f({\boldsymbol J})$| model given by equation (34) with (α, β) = (1, 4), which hereafter we refer to as |$f({\boldsymbol J})$| Hernquist model.

In Appendix A, we derive an analytic expression for Jr = Jr(H, L) in the spherical Hernquist potential. By numerically inverting this expression, we arrive at H = H(Jr, L) for the Hernquist sphere. Combining this with the sphere's ergodic DF, which was given already by Hernquist (1990), we have the exact |$f=f[H({\boldsymbol J})]$|⁠. In Fig. 1, we show surfaces in action space on which this DF is constant together with surfaces on which DF of the |$f({\boldsymbol J})$| Hernquist model is constant. The differences are small but apparent and arise because the surfaces of constant energy are not exactly planar.

The red full curves show surfaces on which the DF of the classical isotropic Hernquist sphere is constant in the (Jr, L) plane of action space, while the black dashed curves show surfaces on which the corresponding $f({\boldsymbol J})$ DF is constant.
Figure 1.

The red full curves show surfaces on which the DF of the classical isotropic Hernquist sphere is constant in the (Jr, L) plane of action space, while the black dashed curves show surfaces on which the corresponding |$f({\boldsymbol J})$| DF is constant.

Open in new tabDownload slide
Fig. 2 compares the radial profiles of density, circular speed and radial component of velocity dispersion in the exact isotropic model and in the |$f({\boldsymbol J})$| Hernquist model. The largest discrepancy is in the velocity dispersion and reflects the fact that the model is significantly radially biased around r0. The long-dashed curve in Fig. 3 shows that the |$f({\boldsymbol J})$| Hernquist model has a slight radial bias at all radii by plotting the anisotropy parameter
(39)
By virtue of the adopted form of g (equation 35), βa → 0 in the Keplerian regime. Even though the potential is not harmonic at the centre, still the model tends to isotropy also at small radii, which justifies our simple choice for |$h({\boldsymbol J})$|⁠.
Density (left-hand panel), circular velocity (central panel), and radial velocity dispersion (right-hand panel) profiles for the classical isotropic Hernquist sphere (normalized to rb) and for the $f({\boldsymbol J})$ Hernquist model (normalized to r0).
Figure 2.

Density (left-hand panel), circular velocity (central panel), and radial velocity dispersion (right-hand panel) profiles for the classical isotropic Hernquist sphere (normalized to rb) and for the |$f({\boldsymbol J})$| Hernquist model (normalized to r0).

Open in new tabDownload slide
Anisotropy profiles for $f({\boldsymbol J})$ Hernquist, $f({\boldsymbol J})$ Jaffe, $f({\boldsymbol J})$NFW, $f({\boldsymbol J})$ isochrone and $f({\boldsymbol J})$ isothermal models. The profiles are normalized to r0 (equation 38).
Figure 3.

Anisotropy profiles for |$f({\boldsymbol J})$| Hernquist, |$f({\boldsymbol J})$| Jaffe, |$f({\boldsymbol J})$|NFW, |$f({\boldsymbol J})$| isochrone and |$f({\boldsymbol J})$| isothermal models. The profiles are normalized to r0 (equation 38).

Open in new tabDownload slide

3.2.2 The Jaffe model

The Jaffe (1983) model behaves as Hernquist's at large radii, while tending to ρ ∝ r−2 close to the centre. Fig. 4 shows the radial profiles of the |$f({\boldsymbol J})$| Jaffe model defined by setting (α, β) = (2, 4) in the DF (34), and compares them with the classical isotropic model. The discrepancies in σr are due to the slight radial bias of the |$f({\boldsymbol J})$| model around r0. The full curve in Fig. 3 shows that this bias is actually quite mild – |βa| < 0.1.

Same as Fig. 2, but for the classical isotropic Jaffe sphere and for the $f({\boldsymbol J})$ Jaffe model.
Figure 4.

Same as Fig. 2, but for the classical isotropic Jaffe sphere and for the |$f({\boldsymbol J})$| Jaffe model.

Open in new tabDownload slide

3.2.3 NFW halo

The NFW model has β = 3 with the consequence that its mass diverges logarithmically as r → ∞ and its potential is never Keplerian. Consequently, the reasoning used to construct a DF above equation (34) does not apply. If we nevertheless adopt equation (34) with (α, β) = (1, 3), we obtain a DF that implies that as J → ∞ the mass with actions less than J diverges like log J. Asymptotically, the circular speed of the standard NFW model is
(40)
so in this model the action of a circular orbit is |$J_\phi \sim \sqrt{r\log r}$|⁠. This shows that mass diverging like log J in action space corresponds, to leading order, to divergence of the mass in real space like log r. Hence, it is plausible that the DF (34) with (α, β) = (1, 3) generates a model similar to the NFW model.
Computation of ρ(r) for the |$f({\boldsymbol J})$| model with (α, β) = (1, 3) bears out this expectation. However, the slope of the model's density profile at large r is slightly steeper than desired, and a better fit to the classical NFW profile is obtained by adopting
(41)
Fig. 5 shows the radial profiles of the classical NFW model and those of the model generated by the DF (41), which we shall call |$f({\boldsymbol J})$|NFW model. The dotted curve in Fig. 3 shows that this model is mildly radially biased at radii larger than r0 and it becomes very slightly tangentially biased for r < r0. These anisotropies account for the difference between the σr profiles of the |$f({\boldsymbol J})$| and classical NFW models.
Same as Fig. 2, but for the classical isotropic NFW sphere and for the $f({\boldsymbol J})$NFW model defined by equation (41).
Figure 5.

Same as Fig. 2, but for the classical isotropic NFW sphere and for the |$f({\boldsymbol J})$|NFW model defined by equation (41).

Open in new tabDownload slide

4 CORES AND CUTS

In the last section, we addressed the problematic nature of power-law models – that their mass diverges at either small or large radii – by introducing separate slopes of the dependence of f on |${\boldsymbol J}$| at small and large |${\boldsymbol J}$|⁠. The recovered models had central density cusps similar to those of the Hernquist, Jaffe and NFW models. If a homogeneous core is required, the natural DF to adopt is
(42)
for then the phase-space density has the finite value |$M_0/J_0^3$| at the centre of the model, and the asymptotic density profile is expected to be ρ ∝ r−β. For β ≤ 3 the system has infinite mass, so for these models we taper the DF by subtracting a constant from the value given by equation (42)
(43)
where |${\boldsymbol J}_{\rm t}$| is some large action, which defines a truncation radius
(44)

4.1 Isochrone model

Fig. 6 compares the density profiles of the model equation (42) generates for β = 4 (black curves) with those of the isochrone (Hénon 1960). The two models are extremely similar, so we shall refer to the model generated by the DF (42) when β = 4 as the |$f({\boldsymbol J})$| isochrone model. The density profiles of the two models are essentially identical, but at r ≃ r0 σr is slightly smaller in the |$f({\boldsymbol J})$| isochrone than in the classical isochrone because the |$f({\boldsymbol J})$| isochrone is mildly radially biased near r0 – the thin full curve in Fig. 3 shows βa(r) for this model. It is non-zero because action space surfaces of |$f({\boldsymbol J})$| do not quite coincide with surfaces of constant |$H({\boldsymbol J})$|⁠, as the upper panel of Fig. 7 shows by plotting contours of f and H. For the isochrone potential, we have an analytic expression for the frequency ratio Ωϕ/Ωr as a function of L. The lower panel of Fig. 7 shows that the constant-energy and constant-DF contours are more closely aligned when the argument of the homogeneous function uses the exact frequency ratio.

Same as Fig. 2, but for the classical isotropic isochrone sphere and for the $f({\boldsymbol J})$ isochrone model.
Figure 6.

Same as Fig. 2, but for the classical isotropic isochrone sphere and for the |$f({\boldsymbol J})$| isochrone model.

Open in new tabDownload slide
Surfaces of constant $H({\boldsymbol J})$ (red) and of constant $f({\boldsymbol J})$ black in action space for two $f({\boldsymbol J})$ isochrone models with different choice of the function g appearing in the DF (42). In the upper panel, g is Jr + L whereas in the lower panel it is Jr + (Ωϕ/Ωr)L where the frequency ratio is a function of L.
Figure 7.

Surfaces of constant |$H({\boldsymbol J})$| (red) and of constant |$f({\boldsymbol J})$| black in action space for two |$f({\boldsymbol J})$| isochrone models with different choice of the function g appearing in the DF (42). In the upper panel, g is Jr + L whereas in the lower panel it is Jr + (Ωϕ/Ωr)L where the frequency ratio is a function of L.

Open in new tabDownload slide

Given that the exact DF of the isochrone is a complicated function of |${\boldsymbol J}$|⁠, it is astonishing that the trivial DF (42) provides such a good approximation to it.

4.2 Cored isothermal sphere

In Section 2.1, we derived an approximation (24) to the DF of the singular isothermal sphere. Here, we modify this model into one that is numerically tractable by (i) adding a core and (ii) tapering its density at large radii so the model's mass becomes finite. Then the DF is
(45)
where |$h({\boldsymbol J})$| is given by equation (10) with both frequency ratios set to 1/√2 and |${\boldsymbol J}_{\rm t}=(0,v_{\rm c}r_{\rm t},0)$|⁠. As full curves in Fig. 8 show this DF generates a model that has a core that extends to r0 and a density profile that plunges to zero near the truncation radius rt. The short-dashed curve in Fig. 3 shows the model's anisotropy parameter βa, which is always small (|βa| < 0.04).
Same as Fig. 2, but for the truncated isotropic cored isothermal sphere (equation 46) and for the $f({\boldsymbol J})$ truncated isothermal model. We show the location of the truncation radius defined by equation (44).
Figure 8.

Same as Fig. 2, but for the truncated isotropic cored isothermal sphere (equation 46) and for the |$f({\boldsymbol J})$| truncated isothermal model. We show the location of the truncation radius defined by equation (44).

Open in new tabDownload slide
An ergodic model with a simple functional form of ρ(r) to which we can compare our |$f({\boldsymbol J})$| model has
(46)
The dashed curve in the left-hand panel of Fig. 8 shows that the model defined by the DF (46) provides an excellent fit to the density profile of our |$f({\boldsymbol J})$| model. Curiously, in the |$f({\boldsymbol J})$| model σr(r) is more nearly constant within rt than in either of the models with analytic density profiles. The dashed curve in the right-hand panel of Fig. 8 shows that the model defined by the DF (46) has a significantly deeper central depression in σr than the |$f({\boldsymbol J})$| model.

5 CONCLUSIONS

Studies of both our own and external galaxies will benefit from the availability of a flexible array of dynamical models of galactic components such as disc, bulge and dark halo. The construction of general models of this type is rather straightforward when one decides to start from an expression for the component's DF as a function of the action integrals Ji. In this paper, we have illustrated this fact by deriving simple analytic forms for DFs that self-consistently generate models that closely resemble the isochrone, Hernquist, Jaffe, NFW and truncated isothermal models. In previous papers, Binney (2010, 2012b) has given simple analytic DFs that provide excellent fits to the structure of the Galactic disc, so now DFs are available for all commonly occurring galactic components.

Our models are tailored to minimize velocity anisotropy at both small and large radii. In all of them, the anisotropy parameter βa peaks at intermediate radii. The peak is by far sharpest in the |$f({\boldsymbol J})$| isochrone, but even in this model βa stays below 0.25.

Our presentation has been elementary in the sense that we have confined ourselves to spherical, almost isotropic components that live in isolation. However, B14 showed that given a near-ergodic DF |$f({\boldsymbol J})$| of a component such as those presented here, it is trivial to modify it so it generates a system that is flattened by velocity anisotropy, or by rotation, or by a combination of the two. Equally important, when the DF of an individual component is given as |$f({\boldsymbol J})$|⁠, it is straightforward to add components. Such addition was exploited by Piffl et al. (2014) in a study of the contribution of dark matter to the gravitational force on the Sun: in that study the models fitted to data comprised a sum of DFs |$f({\boldsymbol J})$| for the disc and the stellar halo. The dark halo was assigned a density distribution rather than a DF, but Piffl et al. (in preparation) represent the dark halo by the |$f({\boldsymbol J})$|NFW model, making the Galaxy a completely self-consistent object. A key point for such work is that the mass of each component can be specified at the outset.

Our approach has several points of contact with that of Williams et al. (2014) and Evans & Williams (2014), who derive approximations to |$H({\boldsymbol J})$| for models that are defined by DFs of the form f(E, L). In particular, they show that for their models better approximations to the iso-energy surfaces in action space can be obtained if one's homogeneous function has as its argument the sum of a linear function of the actions, as used here, and a small term |$\epsilon \sqrt{LJ_r}$|⁠. We expect that the anisotropy of our models could be enhanced by adding such a term.

In addition to assisting in the dynamical interpretation of observations of galaxies, the models that this work makes possible could provide useful initial conditions for N-body simulations. The first step would be the construction of a self-consistent galaxy model from a judiciously chosen DF. Then one could Monte Carlo sample the action space using the DF as the sampling density, and torus mapping (e.g. Binney & McMillan 2011) could be used to generate an orbital torus at each of the selected actions. Finally, some number n of initial conditions |$({\boldsymbol x},{\boldsymbol v})$| would be selected on each torus, uniformly space in the angles θi. The resulting simulation would be in equilibrium to whatever precision had been used in the solution of Poisson's equation, and it would experience a ‘cold start’ (Sellwood 1987). Moreover, given that it would be possible to evaluate the original DF at any phase-space point, the model would lend itself to the method of perturbation particles (Leeuwin, Combes & Binney 1993) in which the simulation particles represent the difference between a dynamically evolving model and an underlying equilibrium rather than the whole model. This method has been little used in the past on account of the lack of interesting models with known DFs, which is precisely the need that we have here supplied.

LP is pleased to thank the Rudolf Peierls Centre for Theoretical Physics in Oxford for the warm hospitality during an early phase of this work. JB is supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 321067, and by the UK Science Technology through grant ST/K00106X/1. LC and CN are partly supported by PRIN MIUR 2010-2011, project ‘The Chemical and Dynamical Evolution of the Milky Way and Local Group Galaxies’, prot. 2010LY5N2T.

1

In the non-spherical case,

(11)
Consequently, ρ and Φ have simple scalings with r but they are not necessarily functions of each other.

REFERENCES

Arnold
V. I.
, 
The Mathematical Methods of Classical Mechanics
,
1978
New York
Springer

Google Scholar

OpenURL Placeholder Text

 
Binney
J.
, 
MNRAS
,
2010
, vol.
401
pg.
231
Crossref
Search ADS
 
Binney
J.
, 
MNRAS
,
2012a
, vol.
426
pg.
1324
Crossref
Search ADS
 
Binney
J.
, 
MNRAS
,
2012b
, vol.
426
pg.
1328
Crossref
Search ADS
 
Binney
J.
, 
MNRAS
,
2014
, vol.
440
pg.
787
 
(B14)
Crossref
Search ADS
 
Binney
J.
McMillan
P. J.
, 
MNRAS
,
2011
, vol.
413
pg.
1889
Crossref
Search ADS
 
Binney
J.
Tremaine
S.
, 
Galactic Dynamics
,
2008
2nd edn
Princeton, NJ
Princeton Univ. Press

Google Scholar

OpenURL Placeholder Text

 
Byrd
P. F.
Friedman
M. D.
, 
Handbook of Elliptic Integrals for Engineers and Scientists
,
1971
Berlin
Springer-Verlag

Google Scholar

OpenURL Placeholder Text

 
Ciotti
L.
, 
ApJ
,
1996
, vol.
471
pg.
68
Crossref
Search ADS
 
Dehnen
W.
, 
MNRAS
,
1993
, vol.
265
pg.
250
Crossref
Search ADS
 
Dickson
L. E.
, 
Elementary Theory of Equations
,
1914
Hoboken, NJ
Wiley Inc.

Google Scholar

OpenURL Placeholder Text

 
Eddington
A. S.
, 
MNRAS
,
1916
, vol.
76
pg.
572
Crossref
Search ADS
 
Evans
N. W.
, 
MNRAS
,
1994
, vol.
267
pg.
333
Crossref
Search ADS
 
Evans
N. W.
Williams
A. A.
, 
MNRAS
,
2014
, vol.
443
pg.
791
Crossref
Search ADS
 
Fermani
F.
, 
DPhil thesis
,
2013
Univ. Oxford

Google Scholar

OpenURL Placeholder Text

 
Gantmacher
F. R.
, 
The Theory of Matrices Vol. 2
,
1959
Providence, Rhode Island
Am. Math. Soc.

Google Scholar

OpenURL Placeholder Text

 
Gerhard
O. E.
Saha
P.
, 
MNRAS
,
1991
, vol.
251
pg.
449
Crossref
Search ADS
 
Hénon
M.
, 
Ann. Astrophys.
,
1960
, vol.
23
pg.
474
 
Hernquist
L.
, 
ApJ
,
1990
, vol.
356
pg.
359
Crossref
Search ADS
 
Jaffe
W.
, 
MNRAS
,
1983
, vol.
202
pg.
995
Crossref
Search ADS
 
Jeans
J. H.
, 
MNRAS
,
1915
, vol.
76
pg.
70
Crossref
Search ADS
 
Leeuwin
F.
Combes
F.
Binney
J.
, 
MNRAS
,
1993
, vol.
262
pg.
1013
Crossref
Search ADS
 
Navarro
J. F.
Frenk
C. S.
White
S. D. M.
, 
ApJ
,
1996
, vol.
462
pg.
563
Crossref
Search ADS
 
Piffl
T.
et al. , 
MNRAS
,
2014
, vol.
445
pg.
3133
Crossref
Search ADS
 
Prendergast
K. H.
Tomer
E.
, 
AJ
,
1970
, vol.
75
pg.
647
Crossref
Search ADS
 
Rowley
G.
, 
ApJ
,
1988
, vol.
331
pg.
124
Crossref
Search ADS
 
Sanders
J. L.
Binney
J.
, 
MNRAS
,
2014
, vol.
441
pg.
3284
Crossref
Search ADS
 
Sanders
J. L.
Binney
J.
, 
MNRAS
,
2015
 
in press
 
Sellwood
J. A.
, 
ARA&A
, vol.
25
pg.
151
Crossref
Search ADS
 
Williams
A. A.
Evans
N. W.
Bowden
A. D.
, 
MNRAS
,
2014
, vol.
442
pg.
1405
Crossref
Search ADS
 
Wilson
C. P.
, 
ApJ
,
1975
, vol.
80
pg.
175
Crossref
Search ADS
 

APPENDIX A: ANALYTICAL EXPRESSION FOR THE RADIAL ACTION IN THE HERNQUIST SPHERE

The radial action is defined as
(A1)
where Φ(r) = −GM/(r + rb) and r1, r2 are the pericentric and apocentric radii for the given energy E and angular momentum L, i.e. the two roots of the integrand in equation (A1). Introducing the dimensionless quantities s ≡ r/rb, |$\mathcal {E}\equiv -Er_{\rm b}/GM$| and |$l = L/\sqrt{2GMr_{\rm b}}$|⁠, equation (A1) can be rewritten as
(A2)
where Ψ(s) ≡ 1/(1 + s) is the relative dimensionless potential. We now change the integration variable from s = (1 − Ψ)/Ψ to Ψ (Ciotti 1996), and have
(A3)
where Ψ1 ≡ Ψ(s1), Ψ2 ≡ Ψ(s2) and
(A4)
is a cubic in Ψ, the roots of which can be found by standard methods (e.g. Dickson 1914). Ψ1, Ψ2 are two roots in the physical range 0 ≤ Ψ ≤ 1. Let A be the third real root, so
(A5)
While it is physically obvious that two of the three real solutions of equation (A4) are in the range (0, 1) and the remaining one is outside (A > 1), we remark that the same conclusion can be reached by purely algebraic arguments by using the Routh–Hurwitz theorem (see e.g. Gantmacher 1959). By evaluating equations (A4) and (A5) at Ψ = 0 one gets |$A=\mathcal {E}/\Psi _1\Psi _2 > 0$|⁠. By splitting into its partial fractions the integrand in equation (A3), it is possible to express the integral for Jr in terms of complete elliptic integrals (see Byrd & Friedman 1971):
(A6)
where K, E, Π are, respectively, the complete elliptic integral of the first, second and third kind,
(A7)
and finally
(A8)
We have tested the formula (A6) for consistency by numerically integrating equation (A1) for a large set of orbits at different (E, L) and the numerical and analytical results agree within the error of the employed routine.
© 2015 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
Issue Section:
Article
Download all slides