Flattened stellar systems based on distribution functions depending on actions

TL;DR

This paper investigates the construction of flattened stellar systems using distribution functions dependent on action integrals, addressing continuity constraints and presenting an algorithm for model creation.

Contribution

It introduces a new algorithm to ensure physically consistent distribution functions for flattened stellar systems based on actions.

Findings

Developed an algorithm for constructing flattened models

Ensured velocity distribution continuity near symmetry axes

Created various self-consistent flattened stellar models

Abstract

We address an issue that arises when self-consistently flattened dynamical stellar systems are constructed by adopting a distribution function (DF) that depends on the action integrals. The velocity distribution at points on the symmetry axis is controlled by the the dependence of the DF on just one action, while at points off the symmetry axis two actions are involved. Consequently, the physical requirement that the velocity distribution evolves continuously in the neighbourhood of the symmetry axis restricts the functional forms of acceptable DFs. An algorithm for conforming to this restriction is presented and used to construct a variety of flattened models.

Equations

$${\nabla }^2\Phi (\mathbf{x})=4\pi GM\int {\mathrm{d}}^3\mathbf{v}f[\mathbf{J}(\mathbf{x},\mathbf{v})].$$

$$\frac{\partial f/\partial J_{\varphi }}{\partial f/\partial J_r}=\frac{{\Omega }_{\varphi }}{{\Omega }_r}.$$

$$\underset{|\mathbf{J}|\to 0}{\lim }\frac{\partial f/\partial J_{\varphi }}{\partial f/\partial J_r}=\frac{1}{2}.$$

$$\underset{J_{\varphi }\to 0}{\lim }\frac{\partial f/\partial J_{\varphi }}{\partial f/\partial J_r}=\frac{1}{2}at\; small{J}_z$$

$$\underset{J_{\varphi }\to 0}{\lim }\frac{\partial f/\partial J_{\varphi }}{\partial f/\partial J_z}=\underset{J_{\varphi }\to 0}{\lim }\frac{{\Omega }_{\varphi }}{{\Omega }_z}.$$

$$\underset{J_{\varphi }\to 0}{\lim }\frac{\partial f/\partial J_{\varphi }}{\partial f/\partial J_z}=1.$$

$$k(\mathbf{J})=J_r+{\eta }_{\varphi }|J_{\varphi }|+{\eta }_z{J}_z,$$

$$\frac{\partial f/\partial J_i}{\partial f/\partial J_j}=\frac{\partial k/\partial J_i}{\partial k/\partial J_j}=\frac{{\eta }_i}{{\eta }_j},$$

$$k(\mathbf{J})=J_r+\xi (\mathbf{J},{\eta }_{\varphi },{\eta }_z)|J_{\varphi }|+{\eta }_z{J}_z.$$

$$\xi (\mathbf{J},{\eta }_{\varphi },{\eta }_z)\to \{\begin{array}{ll}{\eta }_z& J_{\varphi }/|\mathbf{J}|\to 0\\ \frac{1}{2}& |\mathbf{J}|\to 0\\ {\eta }_{\varphi }& \text{otherwise}.\end{array}$$

$$\xi (\mathbf{J},\eta ,\eta )=\eta \forall \mathbf{J}$$

$$\xi({\bf J},\eta_{\phi},\eta_{z})=m_{1}\eta_{\phi}+(1-m_{1})\bigl{[}m_{2}\eta_{z}+(1-m_{2})\eta_{0}\bigr{]},$$

$$m_1(J_{\varphi })\to \{\begin{array}{ll}0& |J_{\varphi }|\to 0\\ 1& \text{otherwise},\end{array}$$

$$m_2(J_{\varphi },J_z)\to \{\begin{array}{ll}0& (J_{\varphi },J_z)\to (0,0)\\ 1& \text{otherwise}.\end{array}$$

$${\eta }_0({\eta }_{\varphi },{\eta }_z)\to \{\begin{array}{ll}\eta & \text{if }{\eta }_{\varphi }={\eta }_z\equiv \eta \\ \frac{1}{2}& \text{if }{\eta }_z>{\eta }_{\varphi }.\end{array}$$

$$m_1(J_{\varphi })$$

$$m_2(J_{\varphi },J_z)$$

$${\eta }_0({\eta }_{\varphi },{\eta }_z)$$

$$\xi=\frac{1}{1+\widetilde{J}_{\phi}^{2}}\Biggl{(}\eta_{\phi}\widetilde{J}_{\phi}^{2}+\frac{\eta_{0}+\eta_{z}(\widetilde{J}_{\phi}^{2}+\widetilde{J}_{z}^{2})}{1+\widetilde{J}_{\phi}^{2}+\widetilde{J}_{z}^{2}}\Biggr{)},$$

$$f({\bf J})=f_{0}\exp\biggl{[}-\biggl{(}\frac{k({\bf J})}{J_{0}}\biggr{)}^{\alpha}\biggr{]},$$

$$k(\mathbf{J})=J_r+k_{\varphi }|J_{\varphi }|+{\eta }_z{J}_z,$$

$$k_{\varphi }=\{\begin{array}{ll}{\eta }_{\varphi }& \text{old models (equation reffor:kj)}\\ \xi (\mathbf{J},{\eta }_{\varphi },{\eta }_z)& \text{new models (equation reffor:newB14)}.\end{array}$$

$$\frac{\partial\ln\rho}{\partial\ln R}\biggr{|}_{(R=r_{\rm c},z=0)}=-\frac{1}{2}.$$

$$\gamma (R,z)\equiv {\left(\frac{\partial z}{\partial R}\right)}_{\rho =\mathrm{const}},$$

$$\delta (u)\equiv 1-\frac{\rho (u,0)}{\rho (0,u)}.$$

$$\begin{gathered} \underset{v_{\varphi }\to 0^{+}}{\lim }\frac{\partial }{\partial v_{\varphi }}\int \mathrm{d}{v}_R\mathrm{d}{v}_{z}f= \\ -\underset{v_{\varphi }\to 0^{-}}{\lim }\frac{\partial }{\partial v_{\varphi }}\int \mathrm{d}{v}_R\mathrm{d}{v}_{z}f<0. \end{gathered}$$

$$\chi^{2}=\sum_{i=0}^{N}\biggl{[}1-\biggl{(}\frac{R_{i}^{2}}{a^{2}}+\frac{z_{i}^{2}}{c^{2}}\biggr{)}\biggr{]}^{2},$$

$$\frac{\partial {\chi }^2}{\partial a^2}=\frac{\partial {\chi }^2}{\partial c^2}=0,$$

$$a^2=\frac{(\sum R_i^2{z}_i^2{)}^2-\sum R_i^4\sum z_i^4}{\sum R_i^2{z}_i^2\sum z_i^2-\sum R_i^2\sum z_i^4},$$

$$c^2=\frac{(\sum R_i^2{z}_i^2{)}^2-\sum R_i^4\sum z_i^4}{\sum R_i^2{z}_i^2\sum R_i^2-\sum R_i^4\sum z_i^2}.$$

Full text

• Galaxy: kinematics and dynamics - Galaxy: structure

1 Introduction

Model stellar systems with known distribution functions (DFs) $f({\bf x},{\bf v})$ are powerful tools for the interpretation of observations because such a model predicts the outcome of any observation. For example, the system’s surface brightness can be obtained by integrating $f$ times the luminosity per star over velocities and the line of sight; a map of any velocity moment can be obtained by including an appropriate power of ${\bf v}$ in the integral. If the system’s stars are resolved, the likelihood of the data given the model can be computed regardless of how many phase-space coordinates have been measured and with what precision.

In general it is essential to consider a stellar system to comprise several components. For example, a globular cluster will contain a range of stellar masses and include significant numbers of massive, dark, remnants. In addition to these components, galaxies are thought to contain a population of dark-matter particles that are distributed more extensively in phase space than the stars, and they generally contain several stellar populations that differ by age and/or chemical composition.

A small number of model systems are known that have analytic expressions for density $\rho({\bf x})$ and potential $\Phi({\bf x})$ and also analytic expressions for the DF. All such systems are spherical and the only multi-component models are the rather specialised and complex two-component models described by Ciotti (1999). A wider range of models can be obtained by considering systems with DFs that are analytic functions of integrals of stellar motion but have density and potential distributions that have to be obtained numerically. Traditionally the integrals of stellar motion used as arguments of the DF have been the energy $E=\frac{1}{2}v^{2}+\Phi({\bf x})$ and the magnitude of the angular momentum ${\bf L}={\bf x}\times{\bf v}$. However, the key to producing multi-component and aspherical systems is to take the DF to be a specified function $f({\bf J})$ of the action integrals $J_{i}$. Advantages of using actions as arguments of the DF include

• •

The mass of any component is specified by the DF before the system’s density and potential have been determined.

• •

The addition of an extra component of mass $M_{i}$ changes the density distributions of other components in an intuitive way.

• •

The self-consistently generated potential can be solved for by a stable and rapidly convergent iteration.

Early uses of DFs that depended exclusively on actions examined the structure predicted by DF in assumed potentials (Binney, 2010, 2012). The potential self-consistently generated by a DF $f({\bf J})$ was first obtained by Binney (2014), but only for one-component systems. Piffl et al. (2015) solved for the potential self-consistently generated by a multi-component DF designed to model our Galaxy. Posti et al. (2015), Williams & Evans (2015) and Pascale et al. (2018) explored spherical models of early-type galaxies that are defined by DFs $f({\bf J})$ designed to model both dark haloes and the stellar content of elliptical and dwarf spheroidal galaxies.

In the case of a spherical potential, the actions comprise two components $J_{\phi}$ and $J_{z}$ of the angular momentum ${\bf L}$, and the radial action $J_{r}$. $J_{\phi}$ is the component of ${\bf L}$ about some chosen axis, and $J_{z}=L-|J_{\phi}|$ quantifies the inclination of the orbit with the respect to the chosen axis. The radial action $J_{r}=(2\pi)^{-1}\oint\text{d}r\,p_{r}$ quantifies the amplitude of a star’s radial oscillations.

If the part of $f$ that is an even function of $J_{\phi}$ depends on $J_{\phi}$ and $J_{z}$ only through the combination $L=J_{z}+|J_{\phi}|$, the system’s real-space structure will be spherical. It may, however, have net rotation around the axis that defines $J_{\phi}$: rotation around this axis is encoded in the part of $f$ that is a odd function of $J_{\phi}$. If the part of $f$ that is even in $J_{\phi}$ depends on $J_{\phi}$ and $J_{z}$ other than through the combination $L$, the model will be aspherical. If $f$ decreases with increasing $J_{z}$ faster than it does with increasing $|J_{\phi}|$, the model will be oblate.

In an oblate potential, the existence of actions is not guaranteed, but numerical orbit integrations in plausible oblate potentials reveal that the great majority of orbits are quasiperiodic (Binney & Tremaine, 2008), which implies the existence of three action integrals. For the majority of orbits these actions prove to be minor modifications of the actions $J_{i}$ familiar from the spherical case (Binney & McMillan, 2016). A minority of ‘resonantly trapped’ orbits have actions that are not simply related to the spherical actions, although they can be computed using first-order perturbation theory formulated in terms of the usual $J_{i}$ (Binney, 2016). To a good approximation, the existence of trapped orbits can be neglected when model building since an ensemble of trapped orbits generate very similar predictions for most observables to an ensemble of untrapped orbits (Binney, 2018).

Subject to this proviso regarding the treatment of resonantly trapped orbits, any non-negative function of three variables for which the integral $\int{\rm d}^{3}{\bf J}\,f({\bf J})$ through the positive octant of Cartesian space is finite defines a stellar model of mass $M$, because, given such a function, one can normalise it such that $(2\pi)^{-3}=\int{\rm d}^{3}{\bf J}\,f({\bf J})$ and can then solve for the potential $\Phi$ that satisfies

[TABLE]

This computation is rendered feasible by the Stäckel Fudge (Binney, 2012), which, given any plausible axisymmetric potential $\Phi({\bf x})$, provides approximate formulae for ${\bf J}({\bf x},{\bf v})$. Equation (1) can be solved in $\sim 5$ iterations, starting from any plausible initial potential $\Phi_{0}$, by taking the potential $\Phi_{n+1}$ to be that on the left of the equation with $\Phi_{n}$ used in the Stäckel Fudge on the right (Binney, 2014, hereafter B14).

While any non-negative, normalisable function $f({\bf J})$ defines a logically possible model, B14 already noted that unless candidate functions $f({\bf J})$ are subjected to restrictions, the final model is liable to display physically implausible structure near the origin and/or symmetry axis. Moreover, Piffl et al. (2015) discovered that the simplest DFs for our Galaxy’s dark halo predicted implausibly cusped velocity distributions. The goal of this paper is to elucidate conditions on $f({\bf J})$ that ensure that it will generate a physically plausible model.

In Section 2 we explain the physical origins of the restrictions on $f({\bf J})$ and suggest a way of satisfying them. In Section 3 we illustrate the effectiveness of our proposal by presenting a variety models with and without implementation of our proposal. In Section 4 we ask why galaxies are restricted in the distribution of stars in action space, and Section 5 sums up.

2 Restricting the DF

2.1 Physical motivation

Along the symmetry axis of an axisymmetric galaxy, the two directions that run parallel to the equatorial plane are physically equivalent. If we call these the $x$ and $y$ directions, it follows that at any point on the symmetry axis, the distributions of $v_{x}$ and $v_{y}$ must be identical.

Only orbits with $J_{\phi}=0$ can reach the symmetry axis, so when considering the velocity distribution at points on that axis we can confine attention to the plane $J_{\phi}=0$ of action space, which has axes $J_{r}$ and $J_{z}$. The orbits in this plane fall into two families: boxes (with small values of $J_{z}$) and loops (with $J_{z}$ above a threshold), as illustrated by Fig. 1. The loop orbits do not reach a central section of the symmetry axis; this section is reached by box orbits, which do not visit the part of the symmetry axis that is visited by loops.

In the central section of the axis, $J_{z}$, which quantifies the height of a box orbit, largely quantifies $v_{z}$, while $J_{r}$ quantifies $v_{x}$ and $v_{y}$. Outside the central section of the symmetry axis, $J_{r}$, which quantifies the radial excursions of a loop orbit, largely quantifies $v_{z}$ while $J_{z}$ quantifies $v_{x}$ and $v_{y}$. In both sections of the symmetry axis, the distributions of $v_{x}$ and $v_{y}$ are guaranteed to be identical because they are set by the way $f({\bf J})$ depends on the same argument: $J_{r}$ in the central section and $J_{z}$ further out. If we wish to avoid rapidly changing velocity distributions as we move between the central and outer sections of the symmetry axis, we should relate the way $f$ depends on $J_{r}$ and $J_{z}$ along the line in the plane $J_{\phi}=0$ that has boxes on one side and loops on the other.

Since the velocity distribution should be a continuous function of position, the distributions of $v_{x}$ and $v_{y}$ should be nearly the same if we move a small distance along the $x$ axis from the symmetry axis. Once we are off the symmetry axis, orbits with non-zero $J_{\phi}$ contribute to the kinematics. Fig. 2 shows two orbits with the same small value of $J_{\phi}$. The orbit on the left approaches the symmetry axis away from its central section, while the orbit on the right approaches just this central section. If we start from the central section of the symmetry axis and move parallel to the $x$ axis, it will be orbits like that on the right of Fig. 2 that contribute to the kinematics. A change to $v_{x}$ of such an orbit changes only $J_{r}$, while a change to $v_{y}$ changes both $J_{r}$ and $J_{\phi}$. Unless we restrict the way $f$ depends on $J_{r}$ and $J_{\phi}$, there is no guarantee that the distributions of $v_{x}$ and $v_{y}$ will be similar at our new location and on the symmetry axis.

If we move parallel to the $x$ axis from a point on the symmetry axis that lies outside the central section, it will be orbits like that shown on the left of Fig. 2 that contribute to the kinematics. A change to $v_{x}$ on an orbit of this type only varies $J_{z}$ (which controls the amplitude of oscillations perpendicular to the symmetry axis), while a change to $v_{y}$ mainly changes $J_{\phi}$. Hence the dependence of $f$ on $J_{\phi}$ and $J_{z}$ must be restricted if the condition of approximate isotropy just off the symmetry axis is to be satisfied.

In summary, these arguments show that

• •

At $J_{\phi}=0$ the derivatives of $f$ with respect to $J_{r}$ and $J_{z}$ should be related along a line in the plane $J_{\phi}=0$.

• •

At small $J_{z}$, the derivatives of $f$ with respect to $J_{r}$ and $J_{\phi}$ should be related in the limit $J_{\phi}\to 0$.

• •

At larger $J_{z}$, the derivatives of $f$ with respect to $J_{z}$ and $J_{\phi}$ should be related in the limit $J_{\phi}\to 0$.

In Section 3 we will show that a DF that does not guarantee near $xy$ isotropy all along the minor axis, generates unphysical density distributions.

2.2 Essential restrictions

Velocity isotropy would be guaranteed if the DF were a function $f(H)$ of the Hamiltonian. Then the derivatives of $f$ with respect to $J_{r}$ and $J_{\phi}$ would be in the ratio of orbital frequencies

[TABLE]

For definiteness, we restrict ourselves to cored models. Orbits that are confined to the core will be essentially harmonic, with the consequence that for these orbits $\Omega_{\phi}/\Omega_{r}\simeq\frac{1}{2}$. We conclude that we can ensure that the velocity distribution parallel to the equatorial plane tends smoothly to the mandatory central isotropy by requiring that

[TABLE]

It is not hard to see that satisfaction of the very similar condition

[TABLE]

ensures that the velocity distribution in the $xy$ plane tends smoothly to isotropy as one approaches any point on the central section of the symmetry axis.

Points on the symmetry axis and outside the central section are reached by orbits with $J_{\phi}=0$ but significantly non-zero $J_{z}$. In this region, $J_{r}$ quantifies the vertical velocity component, which is unrestricted, while $J_{\phi}$ and $J_{z}$ quantify the two tangential components of velocity, which should have nearly identical distributions. By the same chain of argument we deployed above, we infer that the condition of approximate isotropy in $v_{x}$ and $v_{y}$ will be satisfied if

[TABLE]

The limiting frequency ratio required here is unity, as one may convince oneself in two ways: (i) use the Torus Mapper (Binney & McMillan, 2016) to compute the frequencies of orbits for diminishing $|J_{\phi}|$; or (ii) recall that $\Omega_{\phi}-\Omega_{z}$ is the frequency at which the orbital plane precesses, and that symmetry requires this frequency to be zero for an orbit that passes right over the pole of the potential. Thus we conclude that velocity isotropy near the symmetry axis requires

[TABLE]

2.3 Implementing the restrictions

Posti et al. (2015) describe a general procedure for constructing DFs that depend on the actions through a function

[TABLE]

that is linear and homogeneous in ${\bf J}$. Without loss of generality, the coefficient of $J_{r}$ can be taken to be one (B14). Posti et al. (2015) confined their attention to the case $\eta_{\phi}=\eta_{z}$ in which $k$ and therefore $f$ become functions of $(J_{r},L)$. The model that $f$ then generates is spherical. B14 had earlier shown in the special case of the isochrone (Hénon, 1960) that if one chooses $\eta_{z}>\eta_{\phi}$, the model generated by $f$ is flattened. This idea was later exploited by Das & Binney (2016) to model the flattened inner stellar halo of our Galaxy.

Since for a DF $f[k({\bf J})]$

[TABLE]

adopting constant values of $\eta_{\phi}$ and $\eta_{z}$ is not consistent with the restrictions derived in Section 2.2. Hence we replace the coefficient of $J_{\phi}$ in $k({\bf J})$ by a function $\xi$:

[TABLE]

We require $\xi$ to be a continuous function of ${\bf J}$ such that

[TABLE]

We further require

[TABLE]

so that, when $\eta\equiv\eta_{\phi}=\eta_{z}$, $k({\bf J})$ depends on $J_{\phi}$ and $J_{z}$ through the total angular momentum $L=J_{z}+|J_{\phi}|$ and the generated model is spherical.

We satisfy the conditions (10) and (11) by writing

[TABLE]

where $m_{1}$ is a smooth function such that

[TABLE]

and $m_{2}$ is a smooth function such that

[TABLE]

Finally, we require the function $\eta_{0}(\eta_{\phi},\eta_{z})$ in equation (12), to which $\xi$ tends as $(J_{\phi},J_{z})\to(0,0)$, to satisfy

[TABLE]

Functions that satisfy conditions (13) to (15) are

[TABLE]

where $J_{0}$ is a scale action characteristic of the system’s core. The form of $\eta_{0}$ specified by equations (16) satisfies the conditions (15): (i) it ensures a continuous transition between flattened models ($\eta_{z}>\eta_{\phi}$) and spherical models ($\eta_{z}=\eta_{\phi}\equiv\eta$); (ii) in the case of even small flattening, it quickly tends to $\frac{1}{2}$. Fig. 3 plots $\eta_{0}$ as a function of $\eta_{z}$ for given values of $\eta_{\phi}$. Substituting for $m_{1}$ and $m_{2}$ from equations (16) allows one to rearrange equation (12) to

[TABLE]

where $\widetilde{J}_{\phi}\equiv|J_{\phi}|/J_{0}$ and $\widetilde{J}_{z}\equiv J_{z}/J_{0}$ are dimensionless actions.

3 Worked examples

We illustrate the benefit of using $\xi$ (equation 19) rather than $\eta_{\phi}$ as the coefficient of $J_{\phi}$ in the function $k({\bf J})$ by computing some flattened models derived from the DF that Pascale et al. (2018) introduced to model the Fornax dwarf spheroidal galaxy – Pascale et al. (2019) already explored the range of spherical models that this DF generates. The DF is

[TABLE]

where $\alpha$ is a positive, dimensionless constant that primarily controls the model’s density profile, $J_{0}$ is a scale action that sets the size of the model’s core, and $f_{0}$ is a normalising constant that ensures that $(2\pi)^{-3}=\int{\rm d}^{3}{\bf J}\,f({\bf J})$. Here we assume $\alpha=1$. In the interests of generality, in equation (20) we define

[TABLE]

where

[TABLE]

We will refer to the old models, with constant $\eta_{\phi}$ as coefficient of $J_{\phi}$ in equation (7), as $k_{\phi}=\eta_{\phi}$, and to the new models, with $\xi$ as coefficient of $J_{\phi}$ (equation 19), as $k_{\phi}=\xi$.

Fig. 4 plots various quantities for a model computed with $k_{\phi}=\eta_{\phi}=0.9$. The top panel shows the density profiles along the model’s major (red) and minor (black) axes, with distances scaled to the core radius $r_{\rm c}$, defined as the distance down the major axis at which

[TABLE]

The middle panel of Fig. 4 shows isodensity contours in the meridional plane, and the bottom panel shows the axis ratios of these contours that one obtains by fitting ellipses to the contours as explained in Appendix A. These plots reveal unphysical features that derive from the use of $k_{\phi}=\eta_{\phi}$ in equation (21).

At the left edge of the middle panel, the contours are sloping down to the left, reflecting a depression in the density along the minor axis, and, worse still, a discontinuity in the direction of the normals to isodensity surfaces where they cut the symmetry axis. The bottom panel shows a rapid decrease in the model’s flattening as one approaches the centre, which sends the axis ratio $q\equiv c/a$ through unity to values indicative of prolateness before the centre is reached. This feature can also be seen in the top panel, where, when approaching the center, the density along the $z$-axis becomes larger than the density computed along the $R$-axis of the meridional plane. These unphysical features arise from the failure of the DF (20) with $k_{\phi}=\eta_{\phi}$ rather than $k_{\phi}=\xi$ (equation 21) to respect restriction on the velocity distribution along the minor axis that we derived in Section 2.2.

Fig. 5 compares six models computed with different $k_{\phi}$: in the models shown in the upper row $k_{\phi}=\eta_{\phi}$, while in the models of the lower row $k_{\phi}=\xi$. The models have similar flattenings but their radial bias increases from left to right: their parameters are $(\eta_{\phi},\eta_{z})$ = (0.5,1), (0.7,1.4) and (1,2) and flattening increases with the ratio $\eta_{z}/\eta_{\phi}$, while increasing both $\eta_{\phi}$ and $\eta_{z}$ increases the radial bias of the velocity distribution. All models have $\alpha=1$.

The isodensity contours plotted in the upper panels of Fig. 5 have clearly discontinuous slopes across the minor axis. The contours plotted in the lower panels do not show this unphysical discontinuity, so using $k_{\phi}=\xi$ rather than $k_{\phi}=\eta_{\phi}$ banishes cuspy isodensity surfaces. The left and middle panels of Fig. 6 show the iso-density contours along the minor axis on an enlarged scale for models with $(\eta_{\phi},\eta_{z})=(0.9,1.8)$. The right column in this figure shows the variation with $z$ of the gradients of contours

[TABLE]

at three values of $\rho$: $\rho=0.01\rho_{\rm c},0.002\rho_{\rm c},0.001\rho_{\rm c}$, with $\rho_{\rm c}$ equal to the model central density. Positive values of $\gamma$ indicate a depression along the reference isodensity contour when approaching the minor axis. A depression is not unphysical – ‘peanut’ bulges of disc galaxies have such depressions – but a non-zero value of $\gamma$ as $R\to 0$ is unphysical. Blue curves show $\gamma(R)$ when $k_{\phi}=\eta_{\phi}$, while red curves show $\gamma(R)$ when $k_{\phi}=\xi$. We see that use of $\xi$ rather than $\eta_{\phi}$ ensures that the slopes of isodensity contours are only slightly positive near the axis and vanish on the minor axis.

Fig 7 shows density profiles and axis ratios for models with flattenings that increase from left to right. Panels in the top row show density profiles along the major (red) and minor (black) axes for models computed with $k_{\phi}=\eta_{\phi}$, while below them we show the corresponding models with $k_{\phi}=\xi$. In the top row, the red curves fall below the black curves at $R\la r_{\rm c}$, implying that these generally oblate models have prolate cores. In the panels of the middle row, the red and black curves approach $R=0$ together. The bottom row shows plots of the axis ratios as a function of semi-major axis. The models with $k_{\phi}=\eta_{\phi}$ have central axis ratios significantly greater than unity, whereas the more flattened models with $k_{\phi}=\xi$ have axis ratios that are always less than unity. In the least flattened of the models that uses $k_{\phi}=\xi$ (extreme left panel), $c/a$ does exceed unity at $R<0.3r_{\rm c}$, but this model is so nearly spherical that this tendency to central prolateness is of little significance. The model just to its right becomes marginally prolate at $R<0.2r_{\rm c}$, where the density is extremely close to the central density and even a tiny angular variation in density generates a significant value of $1-|c/a|$.

Further exploration of any tendency to prolateness in the core is facilitated by defining the diagnostic

[TABLE]

The small panels of Fig. 7 are plots of $\delta(u)$ in models in which $k_{\phi}=\eta_{\phi}$. Clearly, in an oblate model $\delta$ should be negative, but the figure shows that $\delta>0$ in the core to an extent that varies with ($\eta_{\phi},\eta_{z}$). Fig. 8 plots for several values of $\eta_{\phi}$ the peak value of $\delta$ as a function of the ratio $\eta_{z}/\eta_{\phi}$ that controls a model’s flattening. The top panel is for models that use $k_{\phi}=\eta_{\phi}$ and the bottom panel is for models that use $k_{\phi}=\xi$. In the top panel, the largest values of $\delta$, and therefore the most prolate cores, occur in models with $\eta_{z}/\eta_{\phi}\simeq 1.1$ i.e., nearly spherical models, as is to be expected, and in models with the largest values of $\eta_{\phi}$ and therefore the greatest radial velocity bias. In the most radially biased model ($\eta_{\phi}=1.1$), the peak in $\delta$ reaches $0.08$ and this model is prolate throughout the core. In the bottom panel, $\delta<0.01$ even in the most radially biased and least flattened model. In the vast majority of models it is much smaller. Thus replacing $\eta_{\phi}$ with $\xi$ essentially resolves the issue of prolate cores in addition to banishing cusped isodensity contours.

3.1 Distribution of $v_{\phi}$

Piffl et al. (2015) encountered a problem with the distribution of azimuthal velocity components that occurs already in spherical models of the B14 type: when the model has a radial bias, a cusp around $v_{\phi}=0$ appears in the $v_{\phi}$ distribution because

[TABLE]

When the model has tangential bias, the cusp is replaced by a dimple, so the above limit is positive. To eliminate this problem, they forced the dependence of $f$ on the actions in the limit $J_{\phi}\to 0$ to mirror that of the Hamiltonian. To achieve this goal with $k$ defined by equation (9), we would have to make $\eta_{z}$ depend on the actions such that it tends to $\Omega_{\phi}/\Omega_{r}$ as $J_{\phi}\to 0$. Fig. 9 explores the consequences of our failure to take this step by showing $v_{\phi}$ distributions at several locations in nine models of varying flattening and radial bias. The approximately isotropic and radially biased models have unexceptionable $v_{\phi}$ distributions. At $R\ga 10r_{\rm c}$ the $v_{\phi}$ distributions of the tangentially biased models do have anomalous central shapes. The double humped nature of these profiles simply reflects the absence of a part of the DF that is odd in $J_{\phi}$; such a component would reinforce one hump at the expense of the other, yielding the skew $v_{\phi}$ distribution of a rotating model. The anomalous features that one would ideally eliminate by the method of Piffl et al. (2015) are the central spikes displayed by two of the tangentially biased models.

The difference between the anomalous features in these $v_{\phi}$ distributions and those presented by Piffl et al. (2015) probably arises because the latter were adapting the DFs of Posti et al. (2015), which generate models with central density cusps, rather than building cored models.

4 Discussion

Action space is simply a way of cataloguing orbits, and one may ask why one is not at liberty to populate those orbits as one pleases by setting the DF to an arbitrary, normalisable, non-negative function of ${\bf J}$. The proper response to this claim is to say ‘yes you can populate the orbits as you please, but there are two reasons why the schemes of population occurring in real systems are restricted.’

The first reason is the requirement for self-consistency. It is possible that certain DFs are not consistent with an essentially integrable self-consistent potential. For example, some of the potentials encountered during an attempt to iterate from $\Phi_{0}$ to a self-consistent $\Phi_{n}$ might have significant chaotic zones. Stars would diffuse through these zones, irreversibly changing $f$.

The second reason why the DFs of real systems will be restricted relates to the manner in which stars are distributed in phase space. Imagine stars being scattered like confetti, or shrapnel from a shell burst at a large number of locations. Then the initial conditions will be smooth functions of ${\bf x}$ and of ${\bf v}$ and the resulting DF $f({\bf J})$ will not be one that induces discontinuities in velocity distributions. Star formation will likewise create stars with a density that is a continuous function of $({\bf x},{\bf v})$.

5 Conclusions

While any non-negative, normalisable function $f({\bf J})$ specifies a self-consistent stellar system, the system it generates will have unphysical features unless the DF satisfies the constraints (4) and (6) on its derivatives. These constraints arise from the requirement that the velocity distribution at points near the symmetry axis should differ little from the (axially symmetric) velocity distributions on the symmetry axis. If the DF is a function only of the Hamiltonian $H$, these constraints are automatically satisfied.

We have experimented with an algorithm that generates DFs for flattened and possibly radially biased models that are consistent with the DF tending to a function of the Hamiltonian as $J_{\phi}$ approaches zero. We have shown that the resulting models are free of the unphysical features near the symmetry axis that disfigure models based on the simpler DFs proposed by B14. The new DFs provide a promising basis for modelling different components of globular clusters, dwarf spheroidal galaxies and galactic bulges.

Appendix A Quantifying flattening

Each iso-density contour is characterised by a collection of $\{R_{i},z_{i}\}$ points, with $i=1,...,N$, such that $\rho(R_{i},z_{i})=\hbox{constant}$. For that contour we define

[TABLE]

where $a$ and $c$ are, respectively, the semi-major and the semi-minor axes for an oblate model, viceversa for a prolate model. We find the best ($a,c$) as the solution of

[TABLE]

which is

[TABLE]

and

[TABLE]