© ESO 2021

1. Introduction

In stellar dynamics, the most fundamental quantity of a self-gravitating system is the phase space distribution function f(t, x, υ). It describes the density of stars in the six-dimensional phase space and contains all the dynamical information. As the mass density ρ(t, x) is obtained by integrating f(t, x, υ) over velocity space, there are an infinite number of dynamical models with different orbital structures that correspond to a given density. In stationary spherical symmetry, various techniques to generate dynamical models with a different anisotropy profile for a given density profile have been presented, for example by Osipkov (1979), Richstone & Tremaine (1984), Dejonghe (1984, 1986, 1989), Merritt (1985b), Cuddeford (1991), Gerhard (1991), and Cuddeford & Louis (1995). Well-known examples of spherical models for which different analytical distribution functions have been derived include the Plummer, Hernquist, and Jaffe models (Jaffe 1983; Merritt 1985b,a; Dejonghe 1987; Hernquist 1990; Cuddeford 1991; Baes & van Hese 2007).

It is customary to ignore possible contributions from matter that is not bound to the stellar system under consideration, hence it is convenient to adopt the symbol ℰ = −E ≥ 0 for the binding energy per unit mass. A characteristic of dynamical models that has drawn some attention is the differential energy distribution N(ℰ), that is, the distribution of mass as a function of ℰ. It is a natural diagnostic for dynamical models that is easily calculated from N-body simulations (e.g., van Albada 1982; Hanyu & Habe 2001; Di Cintio et al. 2013; Errani & Peñarrubia 2020). It has been argued that the differential energy distribution is a fundamental partitioning of a stellar system, as equilibrium stellar dynamical systems are collisionless systems in which all particles retain their energies (Binney 1982; Efthymiopoulos et al. 2007; Hjorth & Williams 2010).

Given that for a fixed mass density dynamical models with a widely varying orbital structure can be generated, one would expect that these different models would also have strongly different differential energy distributions. The opposite turns out to be the case, however, in different case studies of spherical models. Binney & Tremaine (1987) considered two different models with the same Jaffe density profile, an isotropic model, and a model consisting of only radial orbits, and they showed that they have a very similar differential energy distribution. Dejonghe (1987) presented a family of Plummer models with an isotropic anisotropy profile in the central regions, but strong anisotropy in the outer regions. He concluded that the differential energy distribution of the models in this family did not differ drastically, even when comparing models that range from completely radial to strongly tangentially anisotropic in the outer regions. A similar conclusion was drawn by Cuddeford & Louis (1995) based on a similar family of anisotropic Plummer models. Finally, Binney & Tremaine (2008) presented the differential energy distribution for three spherical models with a Hernquist density profile, but with different anisotropy profiles: a radially anisotropic, an isotropic, and a tangentially anisotropic model. While the distribution functions of these models are very different, their differential energy distributions are, again, very similar.

All of these studies came to a similar qualitative conclusion, namely that the differential energy distribution generally does not seem to vary strongly with the anisotropy profile and hence mainly depends on the mass distribution. However, no satisfactory explanation has been provided for this systematic trend. The goal of the present paper is to quantify the similarity of the differential energy distribution of steady-state dynamical models in more detail. In Sect. 2 we analyse the differential energy distribution for a family of Hernquist models that covers an extreme range in the anisotropy profile (ranging from completely radial to completely circular orbits). We show that all these models have the same mean binding energy per unit mass, which contributes to the similarity of their differential energy distributions. In Sect. 3 we generalise this result for all steady-state dynamical models corresponding to a given mass density, and we derive a general relation between the total integrated binding energy and other measures for the total energy content. In Sect. 4 we discuss our results, and Sect. 5 contains a summary.

2. The Hernquist model

We start our analysis with a simple but illustrative example: the Hernquist model. This model, introduced by Hernquist (1990), has a very simple density profile

$$\begin{array}{c}\hfill \rho (r)=\frac{M}{2\pi }\frac{b}{r{(b+r)}^3},\end{array}$$(1)

where M represents the total mass and b a scale radius. The (positive) gravitational potential is

$$\begin{array}{c}\hfill \mathrm{\Psi }(r)=\frac{\mathit{GM}}{b+r}.\end{array}$$(2)

The total potential energy of the Hernquist model is (Hernquist 1990; Baes & Ciotti 2019),

$$\begin{array}{c}\hfill W_{\text{tot}}=-\frac{G{M}^2}{6b}.\end{array}$$(3)

2.1. Constant anisotropy models

Baes & Dejonghe (2002) discussed different dynamical models that all generate the Hernquist model. They considered, amongst others, a one-parameter family of models with constant anisotropy β. In Appendix A we derive explicit expressions for the differential energy distribution of this family of models. We note that this family covers a very wide range of orbital structure, ranging from a model with all stars on purely radial orbits (β = 1), over an isotropic model (β = 0), to a model with all stars on circular orbits (β = −∞). We note that in accordance with the density slope–anisotropy relation (An & Evans 2006; Ciotti & Morganti 2010; Van Hese et al. 2011), only models with $\beta ⩽\frac{1}{2}$ have positive distribution functions, which means that models with $\beta >\frac{1}{2}$ are nonphysical. We can, however, still formally derive the differential energy distribution for models with $\beta >\frac{1}{2}$, but these are all negative at large binding energies.

In the left panel of Fig. 1 we plot the differential energy distributions for different members of this family of constant anisotropy Hernquist models. The models corresponding to $\beta =\frac{1}{2}$, β = 0, and $\beta =-\frac{1}{2}$ are the same models as shown in Fig. 4.5 of Binney & Tremaine (2008), Fig. 1 covers a wider range in the anisotropy profile.

Fig. 1. Differential energy distributions for Hernquist models with different anisotropy profiles. Left panel: models with a constant anisotropy for different values of the anisotropy parameter β. Right panel: models with an Osipkov-Merritt anisotropy profile for different values of the anisotropy radius ra. We note that the discontinuity in the Osipkov-Merritt distribution function of these models leaves its mark.

The first obvious conclusion is that the differential energy distributions are very similar, in spite of the very different anisotropy of the models. Inspecting the different curves more closely, we note that the shape of the differential energy distribution changes in a systematic way as the anisotropy of the model changes from purely radial to purely circular orbits. N(ℰ) is a decreasing function of ℰ at both low and high binding energies, and an increasing function of β at intermediate binding energies. In other words, tangentially anisotropic models seem to have a slight preference for stars on orbits with more extreme binding energies, whereas radial models slightly prefer average binding energies. On average, the different differential energy distributions seem to have the same mean value.

This qualitative impression can be quantified by looking at the moments of the differential energy distribution. As the differential energy distribution is just the distribution of mass as a function of binding energy, we obviously must have

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{{\mathrm{\Psi }}_0}N(\mathcal{E})\text{d}\mathcal{E}=M,}\end{array}$$(4)

with Ψ₀ the depth of the potential well and M the total mass of the system. All of the differential energy distributions presented in Appendix A and shown in the left panel of Fig. 1 satisfy this normalisation. Even the non-physical differential energy distribution (A.10) of the hypothetical purely radial model satisfies this normalisation.

Next, we look at the mean binding energy per unit mass,

$$\begin{array}{c}\hfill ⟨\mathcal{E}⟩=\frac{1}{M}{\displaystyle {\int }_0^{{\mathrm{\Psi }}_0}N(\mathcal{E})\mathcal{E}\text{d}\mathcal{E}.}\end{array}$$(5)

For all the constant anisotropy Hernquist models, ⟨ℰ⟩ turns out to have the same value, irrespective of the value of β,

$$\begin{array}{c}\hfill ⟨\mathcal{E}⟩=\frac{\mathit{GM}}{4b}.\end{array}$$(6)

This is an interesting finding: in spite of the fact that the anisotropy can be extremely different, ranging from purely radial to purely circular orbits, the mean binding energy is exactly the same. We can express this result in an equivalent way by introducing the total integrated binding energy B_tot as

$$\begin{array}{c}\hfill B_{\text{tot}}\equiv M⟨\mathcal{E}⟩={\displaystyle {\int }_0^{{\mathrm{\Psi }}_0}N(\mathcal{E})\mathcal{E}\text{d}\mathcal{E}.}\end{array}$$(7)

Taking into account expression (3) we find for the family of constant anisotropy Hernquist models the relation

$$\begin{array}{c}\hfill B_{\text{tot}}=-\frac{3}{2}{W}_{\text{tot}}.\end{array}$$(8)

The total integrated binding energy is hence independent of the value of β.

Since the normalisation and the mean value of all differential energy distributions are identical, deviations between them must be found in the higher-order moments only. The standard deviation of the differential energy distribution is calculated as

$$\begin{array}{c}\hfill {\sigma }_{\mathcal{E}}={[\frac{1}{M}{\displaystyle {\int }_0^{{\mathrm{\Psi }}_0}}N(\mathcal{E})(\mathcal{E}-⟨\mathcal{E}⟩{)}^2\text{d}\mathcal{E}]}^{1/2}.\end{array}$$(9)

For the family of Hernquist models with a constant anisotropy, the standard deviation decreases slightly but systematically with increasing β, ranging from 0.0458 for the model with purely circular orbits, over 0.0375 for the boundary model with $\beta =\frac{1}{2}$, to 0.0299 for the nonphysical model with only radial orbits (all in units GM/b). This confirms the systematic shift from a slight excess for extreme binding energies for tangential models to a preference for average binding energies for radially anisotropic models.

2.2. Osipkov-Merritt models

In order to exclude the possibility that the property (6) is a special property of models with a constant anisotropy, we also calculate the differential energy distribution corresponding to Hernquist models with an Osipkov-Merritt anisotropy profile (Osipkov 1979; Merritt 1985b). These models are characterised by an anisotropy parameter β(r) that increases smoothly from 0 at small radii to 1 at large radii, that is the models are isotropic in the centre and completely radially anisotropic in the outskirts. The transition from isotropic to radial is characterised by an anisotropy radius r_a. The distribution function for the Osipkov-Merritt Hernquist model can be calculated analytically (Hernquist 1990; Baes & Dejonghe 2002), but the differential energy distribution needs to be determined numerically.

The result is shown in the right panel of Fig. 1 for different values of the anisotropy parameter. This plot shows a systematic trend when the anisotropy radius increases from r_a = 0, corresponding to the hypothetical model with only radial orbits, to r_a = ∞, corresponding to the isotropic model. The differential energy distributions all asymptotically behave as the pure radial model at low binding energies, and they asymptotically approach the isotropic model at high binding energies. The larger the anisotropy radius, the more the differential energy distribution approximates the one corresponding to the isotropic model.

We have also numerically calculated the moments of the differential energy distribution for these Osipkov-Merritt models. First of all, all the differential energy distributions obviously satisfy the normalisation (4), which is good validation test for the accuracy of our calculations. Secondly, it turns out that they also all satisfy the relation (6), or equivalently relation (8), that is, the mean binding energy per unit mas is independent of the anisotropy radius. The fact that we find the same relation as for the constant anisotropy models strongly suggests that the mean binding energy per unit mass, or equivalently, the total integrated binding energy, of all Hernquist models is exactly the same, independent of the anisotropy profile.

Finally, we can also quantify the difference between the different differential energy distributions through the standard deviation. We find that, for the Osipkov-Merritt Hernquist models, σ_ℰ increases systematically from 0.0299 to 0.0403, in units of GM/b, when the anisotropy radius increases from r_a = 0 (model with only radial orbits) to r_a = ∞ (isotropic model). This is a similar trend as for the constant anisotropy models: more radially anisotropic models tend to have smaller values of σ_ℰ, or in other words, to prefer more average binding energies.

3. The total integrated binding energy of equilibrium dynamical models

The analysis in the previous section suggests that the mean binding energy per unit mass ⟨ℰ⟩, or equivalently the total integrated binding energy B_tot, of spherical dynamical models with a fixed density profile is independent of the details of the anisotropy profile. In this section we show that this hypothesis is true, not only for spherical models but for any steady-state dynamical model. This turns out to be a simple exercise.

To prove this conjecture, we start with the general expression for the differential energy distribution provided by Binney & Tremaine (2008),

$$\begin{array}{c}\hfill N(\mathcal{E})={\displaystyle \int }\text{d}\mathit{x}{\displaystyle \int }\text{d}\mathit{\upsilon }f(\mathit{x},\mathit{\upsilon })\delta (\mathrm{\Psi }(\mathit{x})-\frac{1}{2}{|\mathit{\upsilon }|}^2-\mathcal{E}).\end{array}$$(10)

The total integrated binding energy is thus given by

$$\begin{array}{c}\hfill B_{\text{tot}}={\displaystyle \int }\mathcal{E}\text{d}\mathcal{E}\int \text{d}\mathit{x}{\displaystyle \int }f(\mathit{x},\mathit{\upsilon })\delta (\mathrm{\Psi }(\mathit{x})-\frac{1}{2}{|\mathit{\upsilon }|}^2-\mathcal{E})\text{d}\mathit{\upsilon }.\end{array}$$(11)

Permuting the order of integration yields

$$\begin{array}{c}\hfill B_{\text{tot}}={\displaystyle \int }\text{d}\mathit{x}{\displaystyle \int }f(\mathit{x},\mathit{\upsilon })\text{d}\mathit{\upsilon }{\displaystyle \int }\mathcal{E}\delta (\mathrm{\Psi }(\mathit{x})-\frac{1}{2}{|\mathit{\upsilon }|}^2-\mathcal{E})\text{d}\mathcal{E}\end{array}$$(12)

or

$$\begin{array}{c}\hfill B_{\text{tot}}={\displaystyle \int }\mathrm{\Psi }(\mathit{x})\text{d}\mathit{x}{\displaystyle \int }f(\mathit{x},\mathit{\upsilon })\text{d}\mathit{\upsilon }-\frac{1}{2}{\displaystyle \int }\text{d}\mathit{x}{\displaystyle \int }f(\mathit{x},\mathit{\upsilon }){|\mathit{\upsilon }|}^2\text{d}\mathit{\upsilon }.\end{array}$$(13)

The inner integral in the first term on the right-hand side of this equation equals the mass density, whereas the integral in the second term is the total kinetic energy:

$$\begin{array}{c}\hfill B_{\text{tot}}={\displaystyle \int }\rho (\mathit{x})\mathrm{\Psi }(\mathit{x})\text{d}\mathit{x}-\frac{1}{2}{\displaystyle \int }\rho (\mathit{x})⟨{\upsilon }^2⟩(\mathit{x})\text{d}\mathit{x}.\end{array}$$(14)

The first term is the total potential energy of an equilibrium dynamical system (Binney & Tremaine 2008). With the notations

$$\begin{array}{cc}\hfill W_{\text{tot}}& =-\frac{1}{2}{\displaystyle \int }\rho (\mathit{x})\mathrm{\Psi }(\mathit{x})\text{d}\mathit{x},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill T_{\text{tot}}& =\frac{1}{2}{\displaystyle \int }\rho (\mathit{x})⟨{\upsilon }^2⟩(\mathit{x})\text{d}\mathit{x},\hfill \end{array}$$(16)

we obtain for the total integrated (binding) energy

$$\begin{array}{c}\hfill B_{\text{tot}}=-2W_{\text{tot}}-T_{\text{tot}}.\end{array}$$(17)

Invoking the virial theorem, which states that W_tot + 2 T_tot = 0 or E_tot = W_tot + T_tot = −T_tot < 0, with E_tot the total energy, we ultimately obtain the expression

$$\begin{array}{c}\hfill B_{\text{tot}}=3T_{\text{tot}}=-\frac{3}{2}{W}_{\text{tot}}=-3E_{\text{tot}}.\end{array}$$(18)

This proves our hypothesis that the total integrated binding energy of a steady-state dynamical model is independent of the details of the orbital structure. We also note that the casual argument B_tot = −W_tot − T_tot = −E_tot is erroneous.

4. Discussion

4.1. The similarity of the differential energy distributions

The original starting point of this research was the observation in previous work that the differential energy distribution of steady-state spherical dynamical models generally does not seem to vary strongly with the anisotropy profile and hence mainly depends on the mass distribution. While Binney & Tremaine (1987, 2008), Dejonghe (1987), and Cuddeford & Louis (1995) all concluded that models with fixed density profile all have a similar differential energy distribution, the level of similarity has not been quantified.

The first main conclusion of this work is a quantification of this statement: we demonstrate that all dynamical models corresponding a fixed density distribution have the same mean binding energy per unit mass. Since the normalisation and the mean value of the differential energy distribution of any dynamical model consistent with a given density distribution is identical, deviations between them must be found in the higher-order moments only. This explains, at least to some degree, why the distributions are so similar, in spite of the fact that the orbital structure can vary from completely circular to completely radial orbits.

To the best of our knowledge, this simple fact has never been demonstrated. Dejonghe (1987) does discuss μ_ℰ(r), the mean binding energy per unit mass as a function of radius, in his analysis of a set of anisotropic Plummer models. Fig. 2 from that paper shows how μ_ℰ(r) changes as a function of the anisotropy of the models (parameterised by a parameter q). At large radii, models with a tangential anisotropy have a slightly smaller mean binding energy per unit mass than models with a radial anisotropy, and vice versa at small radii. This anisotropy dependence is also clear from the explicit formula for μ_ℰ(r), as written in Eq. (28) of Dejonghe (1987). Integrating this quantity over the entire physical space results in what we denote as ⟨ℰ⟩. This integral results in ⟨ℰ⟩ = 9π/64, independent of q as required.

4.2. Systematic differences in the differential energy distributions

The second general conclusion from our work relates to these second-order differences between the differential energy distributions corresponding to spherical models with different anisotropy profiles. Both the constant anisotropy and the Osipkov-Merritt Hernquist models presented in Sect. 2 showed a systematic behaviour between the anisotropy profile and the shape of the differential energy, as parametrised by the standard deviation. The bottomline is that more radially anisotropic models tend to prefer more average binding energies, whereas models with a more tangential orbital distribution slightly favour more extreme binding energies. We note that these differences always need to be symmetric, because the mean binding energy is independent of the anisotropy profile.

Inspecting previous work that discussed the differential energy distribution as a function of the anisotropy profile, the same trend is seen. It is most clearly present in the two sets of Plummer models discussed by Dejonghe (1987) and Cuddeford & Louis (1995). As shown in their Fig. 3 and Fig. 3a, respectively, the models with a radial anisotropy have a clearly peaked differential energy distribution with a maximum around ⟨ℰ⟩ = 9π/64, whereas the differential energy distribution of tangential models is much broader and less strongly peaked. For the most strongly tangential models shown by Cuddeford & Louis (1995), the distribution is even decidedly boxy.

The reason for this systematic shift needs to be sought in the different contribution of potential and kinetic energy to the binding energy. In a dynamical model mainly consisting of radial orbits, many particles, especially at large radii, spend a considerable fraction of their time near the apocentre. Their kinetic energy is modest, so the binding energy corresponding to these orbits is roughly equal to the (positive) potential energy, with almost no kinetic energy subtracted. On the other hand, since all particles on strongly radial orbits also cross the inner regions, there is little additional room for highly bound particles on tangential orbits in the central regions. The final result is that radially anisotropic models contain few particles with either very small or very large binding energies, and thus display a preference for average binding energies. Strongly tangential models on the other hand show the opposite behaviour: they contain both tightly bound particles that remain in the very inner regions, whereas the particles in the outer regions still have a significant kinetic energy, almost at the same magnitude of the potential energy. The result is that tangentially anisotropic models prefer more extreme binding energies. Interestingly, this excess at both small and large binding energies is nicely balanced such that the average binding energy is independent of the anisotropy profile.

4.3. The total integrated binding energy

The last result from our study concerns the actual value of the ⟨ℰ⟩, or equivalently of the total integrated binding energy B_tot = M⟨ℰ⟩. It turns out to be quite straightforward to demonstrate that, for any steady-state dynamical model, B_tot can be written as expression (17). We can write B_tot in this form due to the simple fact that the binding energy per unit mass is the sum of the binding potential and the kinetic energy per unit mass, that is, it is a linear combination of these two contributions. Any higher-order moment of the differential energy distributions contains mixed combinations of these contributions. Interestingly, we now have three different linear combinations of the total potential and total kinetic energy of a steady-state dynamical model:

$$\begin{array}{cc}\hfill 2W_{\text{tot}}+T_{\text{tot}}& =-B_{\text{tot}},\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill W_{\text{tot}}+T_{\text{tot}}& =E_{\text{tot}},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill W_{\text{tot}}+2T_{\text{tot}}& =0,\hfill \end{array}$$(21)

which leads to

$$\begin{array}{c}\hfill B_{\text{tot}}=3T_{\text{tot}}=-\frac{3}{2}{W}_{\text{tot}}=-3E_{\text{tot}}.\end{array}$$(22)

The total integrated binding energy thus supplements the well-known trio consisting of total kinetic energy, total potential energy, and total energy on equal footing. The four energies (B_tot, T_tot, W_tot, E_tot) form a quadruple, in which knowledge of any one out of these four suffices to calculate the other three. This result is the more remarkable in view of the fact that knowledge of N(ℰ), say from an N-body calculation, would not seem to suffice in order to calculate the total kinetic and/or potential energies separately. Clearly it does.

These relations can be extended to the case of alternative theories of gravity where the force differs from the Newtonian r⁻¹ potential. The main result, expression (17), does not depend on the nature of the potential. In the case of an interparticle power-law potential of the form r^−p, as studied by many previous authors (e.g., Ispolatov & Cohen 2001; Iguchi 2002; Di Cintio et al. 2013, 2015, 2017; Marcos et al. 2017), the scalar virial theorem can be written as (Binney & Tremaine 2008, Problem 7.1)

$$\begin{array}{c}\hfill p{W}_{\text{tot}}+2K_{\text{tot}}=0.\end{array}$$(23)

Combining this with the expressions (17) and (20) leads to a generalisation of Eq. (22),

$$\begin{array}{c}\hfill B_{\text{tot}}=\left(\frac{4-p}{p}\right)T_{\text{tot}}=-\left(\frac{4-p}{2}\right)W_{\text{tot}}=-\left(\frac{4-p}{2-p}\right)E_{\text{tot}}.\end{array}$$(24)

5. Summary

The starting point for this paper is the qualitative observation that the differential energy distribution of steady-state spherical models does not vary strongly with the anisotropy profile. Our investigation, based on a detailed study of a family of Hernquist models and a rigorous general proof, has yielded the following conclusions.

Firstly we demonstrate that all dynamical models corresponding a fixed density distribution have the same mean binding energy per unit mass. Differences between differential energy distributions are only visible in the higher-order moments. This helps to explain why the distributions are so similar, even for orbital structures ranging from completely circular to completely radial orbits.

Secondly we show that the shape of the differential energy distribution changes systematically with the anisotropy profile for spherical dynamical models with a given density distribution, radially anisotropic models tend to prefer more average binding energies, whereas models with a more tangential orbital distribution slightly favour more extreme binding energies. These differences are always symmetric because the mean binding energy is independent of the anisotropy profile. This systematic effect can be understood as the result of the different contribution of potential and kinetic energy to the binding energy.

Finally, we demonstrate that the total integrated binding energy supplements the well-known triplet consisting of total kinetic energy, total potential energy, and total energy on an equal footing. The four energies form a quadruple, in which knowledge of any one out of these four suffices to calculate the other three.

References

Appendix A: Differential energy distributions for the constant anisotropy Hernquist models

For a spherical dynamical model characterised by a general distribution function f(ℰ,L), expression the differential energy distribution can be written as

$$\begin{array}{c}\hfill N(\mathcal{E})=16{\pi }^2{\displaystyle {\int }_0^{r_{\text{m}}(\mathcal{E})}\text{d}r{\int }_0^{\sqrt{2r^2(\mathrm{\Psi }-\mathcal{E})}}\frac{f(\mathcal{E},L)L\text{d}L}{\sqrt{2(\mathrm{\Psi }-\mathcal{E})-L^2/r^2}},}\end{array}$$(A.1)

with r_m(ℰ) the maximum radius that can be reached by a star with binding energy ℰ, defined through Ψ(r_m(ℰ)) = ℰ. In this Section we derive explicit expressions for the differential energy distribution for Hernquist models with a constant anisotropy, including the limiting cases of models consisting of purely radial or purely circular orbits.

A.1. Isotropic model

For isotropic models, the distribution function is function of ℰ only and can be calculated from the density by means of Eddington’s formula (Eddington 1916). If the distribution function only depends on binding energy, f(ℰ) can be taken outside the integrals in Eq. (A.1) and the inner integral can be evaluated analytically. The result is

$$\begin{array}{c}\hfill N(\mathcal{E})=16\sqrt{2}{\pi }^{2}f(\mathcal{E}){\displaystyle {\int }_0^{r_{\text{m}}(\mathcal{E})}\sqrt{\mathrm{\Psi }-\mathcal{E}}{r}^2\text{d}r.}\end{array}$$(A.2)

For the isotropic Hernquist model, the distribution function can be written in terms of elementary functions, and also the remaining integral in expression (A.2) can be evaluated analytically. The final result is (Hernquist 1990),

$$\begin{array}{cc}\hfill N(\mathcal{E})& =\frac{b}{12\pi G}\frac{1}{{\epsilon }^4{(1-{\epsilon }^2)}^2}\hfill \\ \hfill & \times [\frac{3arcsin\epsilon }{\sqrt{1-{\epsilon }^2}}-\epsilon (1-2{\epsilon }^2)(3-8{\epsilon }^2-8{\epsilon }^4)]\hfill \\ \hfill & \times [\frac{3(1-4{\epsilon }^2+8{\epsilon }^4)arccos\epsilon }{\epsilon }\hfill \\ \hfill & +(1-4{\epsilon }^2)(3+2{\epsilon }^2)\sqrt{1-{\epsilon }^2}].\hfill \end{array}$$(A.3)

where ε is the dimensionless binding energy,

$$\begin{array}{c}\hfill \epsilon =\frac{\mathcal{E}}{\mathit{GM}/b}.\end{array}$$(A.4)

A.2. Models with a constant anisotropy

The isotropic case can be considered as a special case of the more general class of models with distribution function of the form f(ℰ,L) = f_A(ℰ) L^−2β. These models are characterised by an anisotropy parameter β(r) = β that is independent of radius. The isotropic case corresponds to β = 0, models with a radial anisotropy have 0 < β < 1, and models with a preference for tangential orbits have β < 0. Inserting this distribution function in Eq. (A.1) and solving the inner integral, one obtains (Cuddeford 1991; Evans & An 2006)

$$\begin{array}{c}\hfill N(\mathcal{E})=\frac{{(2\pi )}^{5/2}}{2^{\beta -1}}\frac{\mathrm{\Gamma }(1-\beta )}{\mathrm{\Gamma }(\frac{3}{2}-\beta )}{f}_A(\mathcal{E}){\displaystyle {\int }_0^{r_{\text{m}}(\mathcal{E})}\frac{r^{2-2\beta }\text{d}r}{{(\mathrm{\Psi }-\mathcal{E})}^{\beta -1/2}}.}\end{array}$$(A.5)

Baes & Dejonghe (2002) presented an explicit expression for the distribution function of constant anisotropy Hernquist models in terms of hypergeometric functions, and showed that this distribution function is positive, and hence physical, as long as $\beta ⩽\frac{1}{2}$. Inserting this expression into equation (A.5) and evaluating the resulting integral, we find a general expression for the differential energy distribution of the Hernquist model with constant anisotropy β,

$$\begin{array}{cc}\hfill N(\mathcal{E})& =\frac{2b}{G}\frac{\mathrm{\Gamma }(3-2\beta )\mathrm{\Gamma }(5-2\beta )}{\mathrm{\Gamma }(\frac{9}{2}-3\beta )\mathrm{\Gamma }(\frac{7}{2}-\beta )}(1-\epsilon )\hfill \\ \hfill & {\times }_2{F}_1(\frac{1}{2}-\beta ,\frac{3}{2}-\beta ,\frac{9}{2}-3\beta ,1-\epsilon )\hfill \\ \hfill & {\times }_2{F}_1(-\frac{3}{2}+\beta ,\frac{5}{2}+\beta ,\frac{7}{2}-\beta ,\epsilon ).\hfill \end{array}$$(A.6)

For all integer and half-integer values of β, this daunting expression can be written in terms of elementary functions. Setting β = 0 recovers the expression (A.3). A particularly simple case is $\beta =\frac{1}{2}$, the model with the most radial anisotropy that still corresponds to the positive distribution function. The result is (see also Binney & Tremaine 2008, Problem 4.7)

$$\begin{array}{c}\hfill N(\mathcal{E})=\frac{3b}{G}{(1-\epsilon )}^2.\end{array}$$(A.7)

The expressions for other integer and half-integer values of β are more complex. For example, for $\beta =-\frac{1}{2}$ we have

$$\begin{array}{cc}\hfill N(\mathcal{E})& =\frac{b}{3G}\frac{(10-10\epsilon +3{\epsilon }^2)}{{(1-\epsilon )}^4}\hfill \\ \hfill & \times (1-6\epsilon +18{\epsilon }^2-10{\epsilon }^3-3{\epsilon }^4+12{\epsilon }^{3}log\epsilon ).\hfill \end{array}$$(A.8)

A.3. Radial orbit model

In the limit β → 1 we have a dynamical model with only radial orbits. As a radial orbit has zero angular momentum, such models are characterised by a distribution function of the form f(ℰ,L) = f_R(ℰ) δ(L²). Inserting this expression into equation (A.1), we immediately find

$$\begin{array}{c}\hfill N(\mathcal{E})=4\sqrt{2}{\pi }^2{f}_{\text{R}}(\mathcal{E}){\displaystyle {\int }_0^{r_{\text{m}}(\mathcal{E})}\frac{\text{d}r}{\sqrt{\mathrm{\Psi }-\mathcal{E}}}.}\end{array}$$(A.9)

For the case of the Hernquist model, this case is a bit cumbersome, as the Hernquist model cannot be supported by a model consisting of only radial orbits. Indeed, it has only a r⁻¹ cusp, whereas a r⁻² cusp is required for a radial orbit model to be consistent (Richstone & Tremaine 1984). However, we can still formally derive the differential energy distribution for this hypothetical orbital configuration by evaluating the integral (A.9). The result is

$$\begin{array}{c}\hfill N(\mathcal{E})=\frac{16b}{15\pi G}(5-6\epsilon )[arccos\sqrt{\epsilon }+\sqrt{\epsilon (1-\epsilon )}].\end{array}$$(A.10)

This expression can also be obtained by setting β = 1 in the general expression (A.6). This differential energy distribution is negative for $\frac{5}{6}<\epsilon <1$.

A.4. Circular orbit model

Finally, in the limit β → −∞, all stars are on purely circular orbits. The distribution function of such models can be written as (Richstone & Tremaine 1984)

$$\begin{array}{c}\hfill f(\mathit{x},\mathit{\upsilon })=\frac{1}{\pi }\rho (r)\delta (v_r)\delta (v_t+r{\mathrm{\Psi }}^{\prime }).\end{array}$$(A.11)

Inserting this expression in Eq. (10), we find after some calculation

$$\begin{array}{c}\hfill N(\mathcal{E})=4\pi {\displaystyle {\int }_0^{\infty }\rho (r)\delta (\mathrm{\Psi }+\frac{1}{2}r{\mathrm{\Psi }}^{\prime }-\mathcal{E})r^2\text{d}r.}\end{array}$$(A.12)

If the density and the potential are known, this integral can be evaluated and written as a function of ℰ. For the Hernquist model, this yields the following simple expression

$$\begin{array}{c}\hfill N(\mathcal{E})=\frac{2b}{G}(\frac{3}{\sqrt{1+8\epsilon }}-1).\end{array}$$(A.13)

All Figures

Fig. 1. Differential energy distributions for Hernquist models with different anisotropy profiles. Left panel: models with a constant anisotropy for different values of the anisotropy parameter β. Right panel: models with an Osipkov-Merritt anisotropy profile for different values of the anisotropy radius ra. We note that the discontinuity in the Osipkov-Merritt distribution function of these models leaves its mark.

In the text