© O. Bienaymé 2019

1. Introduction

It is well established that orbits in galactic potentials generally admit three integrals of motion (Contopoulos 1960; Ollongren 1962) with the possible presence of ergodic orbits (Hénon & Heiles 1964). The kinematic modelling of a stationary stellar population can therefore be done simply by using a distribution function that depends on three integrals of motion. However, very few gravitational potentials allow such distribution functions to be constructed with analytic integrals. In addition to rotationally symmetric systems with two analytic integrals, meaning the energy and the angular momentum, the most general case known concerns Stäckel potentials with three analytic integrals. These integrals make it possible to model certain potentials with triaxial symmetry¹. Apart from these cases, known potentials with analytic integrals are uncommon (see for example Hietarinta 1987) and have no obvious application for galactic dynamics.

To numerically compute integrals of motion, McGill & Binney (1990) and Sanders & Binney (2014) (see Sanders & Binney 2016, for a review of different methods) showed how the torus overlaid with orbits in phase space could be modelled by numerically determining the angle-action variables. Other related methods (Binney 2012; Sanders & Binney 2015) have been proposed with simpler numerical approaches that are sufficiently precise in many practical cases, either for axisymmetric or triaxially symmetrical potentials. These simpler methods apply for potentials for which the orbits are close to those of Stäckel potentials.

In Sect. 2 of this article we detail exact integrals of motion defined by explicit analytic expressions. Our results are based on the same principles as those of Sanders (2012) and Sanders & Binney (2014), save that the search for integrals is analytic and not numerical. These integrals of motion are exact for Stäckel triaxial potentials and their expressions depend explicitly on the potential. In Sect. 3, we test the efficiency of using these analytic expressions as approximate integrals for non-Stäckel triaxial potentials. We note that the initial principle, that is the efficiency of representing each orbit by using an integral of a Stäckel potential, has been shown by Kent & de Zeeuw (1991). In Sect. 4, we detail a distribution function with a triaxial velocity ellipsoid for a triaxial potential. We test the stationarity and accuracy of this distribution function and of the approximate integrals using the Jeans equation. Section 5 summarizes the results and proposes possible extensions of this work.

2. Triaxial Stäckel potentials and their integrals of motion

Ellipsoidal coordinates and Stäckel potentials have been detailed and discussed many times (Whittaker & Watson 1902; Morse & Feshbach 1953; Ollongren 1962; Lynden-Bell 1962; de Zeeuw & Lynden-Bell 1985; de Zeeuw 1985a,b). Stäckel potentials can be written in ellipsoidal coordinates (see Appendix), depending on three arbitrary functions ζ(λ), η(μ), and θ(ν):

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\lambda ,\mu ,\nu )=& -4\frac{\zeta \left(\lambda \right)(\lambda +\alpha )(\lambda +\beta )(\lambda +\gamma )}{(\lambda -\mu )(\lambda -\nu )}-4\frac{\eta \left(\mu \right)(\mu +\alpha )(\mu +\beta )(\mu +\gamma )}{(\mu -\nu )(\mu -\lambda )}-4\frac{\theta \left(\nu \right)(\nu +\alpha )(\nu +\beta )(\nu +\gamma )}{(\nu -\lambda )(\nu -\mu )}·\hfill \end{array}$$(1)

With this coordinate system, the equations of motion separate and it can be shown that the Stäckel potentials admit three independent integrals. Here, we show that these integrals can be written with an explicit dependence on positions and velocities and on the values of the potential at three different positions. This is obtained by noting that the usual free functions defining the potential can be written as functions of the potential along the three axes of symmetry x, y, and z of the potential. This alternative formulation can be applied to obtain approximate integrals of motion for any potential that is close to a Stäckel potential and for the related orbits. Thus, possible fitting of a given potential V with a Stäckel potential Φ consists in equating both potentials along these three different axes. Other possibilities are also presented.

If Φ is a continuous function at the origin, the formulation in Eq. (1) implies that Φ(−α, −β, −γ) = 0. Evaluating Φ at positions (λ, μ = −β, ν = −γ), we deduce ζ(λ):

$$\begin{array}{c}\hfill \zeta \left(\lambda \right)=-1/4\frac{\mathrm{\Phi }(\lambda ,-\beta ,-\gamma )}{\lambda +\alpha }·\end{array}$$(2)

Similarly, η(μ) and θ(ν) are obtained evaluating Φ(−α,−β,ν) and Φ(−α,μ,−γ):

$$\begin{array}{c}\hfill \eta \left(\mu \right)=-1/4\frac{\mathrm{\Phi }(-\alpha ,\mu ,-\gamma )}{\mu +\beta },\end{array}$$(3)

$$\begin{array}{c}\hfill \theta \left(\nu \right)=-1/4\frac{\mathrm{\Phi }(-\alpha ,-\beta ,\nu )}{\nu +\gamma }·\end{array}$$(4)

Then, substituting these expressions of ζ, η, and θ within the expression of Φ (Eq. (1)), we obtain

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\lambda ,\mu ,\nu )=& \frac{(\lambda +\beta )(\lambda +\gamma )}{(\lambda -\mu )(\lambda -\nu )}\mathrm{\Phi }(\lambda ,-\beta ,-\gamma )+\frac{(\mu +\alpha )(\mu +\gamma )}{(\mu -\nu )(\mu -\lambda )}\mathrm{\Phi }(-\alpha ,\mu ,-\gamma )+\frac{(\nu +\alpha )(\nu +\beta )}{(\nu -\lambda )(\nu -\mu )}\mathrm{\Phi }(-\alpha ,-\beta ,\nu ),\hfill \end{array}$$(5)

a valid expression if Φ(−α, −β, −γ) = 0.

The expressions of the three integrals of motion, the Hamiltonian H, and two integrals J and K are given by de Zeeuw & Lynden-Bell (1985, see Eqs. (7)–(9)). They also detail the integrals I₂ and I₃ that have useful properties particularly for building stellar distribution functions for stellar discs or halos. We use the notations and definitions adopted by de Zeeuw & Lynden-Bell (1985, Eqs. (7) and (14)), and for convenience write the exact integrals in compact form as

$$\begin{array}{c}\hfill H=\mathrm{\Phi }(\lambda ,\mu ,\nu )+1/2{v_x}^2+1/2{v_y}^2+1/2{v_z}^2,\end{array}$$(6)

$$\begin{array}{c}\hfill J=\mathrm{\Psi }(\lambda ,\mu ,\nu )+f_J(L^2,v_x,v_y,v_z),\end{array}$$(7)

$$\begin{array}{c}\hfill K=\mathrm{\Xi }(\lambda ,\mu ,\nu )+f_K(L_x^2,L_y^2,L_z^2,v_x,v_y,v_z),\end{array}$$(8)

and

$$\begin{array}{cc}\hfill I_2=& \frac{{\alpha }^{2}H+\alpha J+L}{\gamma -\alpha }=\psi (\lambda ,\mu ,\nu )+1/2{(\alpha -\beta )v_x}^2+1/2\frac{(\alpha -\beta )}{(\alpha -\gamma )}{L_y}^2+1/2{L_z}^2\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill I_3=& \frac{{\gamma }^{2}H+\gamma J+L}{\alpha -\gamma }=\xi (\lambda ,\mu ,\nu )+1/2(\gamma -\beta ){v_z}^2+1/2{L_x}^2+1/2\frac{(\gamma -\beta )}{(\gamma -\alpha )}{L_y}^2.\hfill \end{array}$$(10)

According to de Zeeuw & Lynden-Bell (1985), who give a complete discussion, “I₂ and I₃ can be considered as generalizations of the energy integrals that exist in a potential separable in Cartesian coordinates, or generalization of angular momentum integrals in axisymmetric or spherical potentials”. For instance, in the limit of an axisymmetric potential (α = β) with the z axis for the rotational symmetry, $I_2=1/2L_z^2$. We note that the third integral is null, I₃ = 0 for motion within the (x, y) plane (i.e. v_z = 0).

After substitution of ζ, η, and θ, we obtain the following expressions depending on the potential at three distinct positions, one on each of the three axes of symmetry of the potential (these expressions are valid only if Φ(−α, −β, −γ) = 0):

$$\begin{array}{cc}\hfill \psi =& \frac{(\mu +\alpha )(\nu +\alpha )}{\alpha -\gamma }[\frac{(\lambda +\beta )(\lambda +\gamma )}{(\lambda -\mu )(\lambda -\nu )}\mathrm{\Phi }(\lambda ,-\beta ,-\gamma )+\frac{(\mu +\gamma )(\lambda +\alpha )}{(\mu -\nu )(\mu -\lambda )}\mathrm{\Phi }(-\alpha ,\mu ,-\gamma )+\frac{(\nu +\beta )(\lambda +\alpha )}{(\nu -\lambda )(\nu -\mu )}\mathrm{\Phi }(-\alpha ,-\beta ,\nu )],\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \xi =& \frac{(\mu +\gamma )(\lambda +\gamma )}{\gamma -\alpha }[\frac{(\lambda +\beta )(\nu +\gamma )}{(\lambda -\mu )(\lambda -\nu )}\mathrm{\Phi }(\lambda ,-\beta ,-\gamma )+\frac{(\mu +\alpha )(\nu +\gamma )}{(\mu -\nu )(\mu -\lambda )}\mathrm{\Phi }(-\alpha ,\mu ,-\gamma )+\frac{(\nu +\alpha )(\nu +\beta )}{(\nu -\lambda )(\nu -\mu )}\mathrm{\Phi }(-\alpha ,-\beta ,\nu )].\hfill \end{array}$$(12)

Other choices of the three positions are possible but require more algebra. However, a simple case can be deduced, again evaluating Φ at three positions; determining ζ, η, and θ; and substituting in Φ, ψ, and ξ expressions, now being defined with the potential at three positions on two planes and one axis of symmetry (for all the expressions below, it is no longer necessary that the potential be zero at the origin):

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\lambda ,\mu ,\nu )=& \frac{(\mu +\gamma )}{(\mu -\nu )}\mathrm{\Phi }(\lambda ,\mu ,-\gamma )+\frac{(\nu +\beta )}{(\nu -\mu )}\mathrm{\Phi }(\lambda ,-\beta ,\nu )+\frac{(\beta -\gamma )}{(\mu -\nu )}\mathrm{\Phi }(\lambda ,-\beta ,-\gamma ),\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill \psi (\lambda ,\mu ,\nu )=& \frac{(\lambda +\alpha )(\mu +\gamma )(\nu +\alpha )}{(\mu -\nu )(\alpha -\gamma )}\mathrm{\Phi }(\lambda ,\mu ,-\gamma )-\frac{(\lambda +\alpha )(\mu +\alpha )(\nu +\beta )}{(\mu -\nu )(\alpha -\gamma )}\mathrm{\Phi }(\lambda ,-\beta ,\nu )\hfill \\ \hfill & -[\frac{(\lambda +\alpha )(\lambda +\beta )(\mu +\gamma )(\nu +\alpha )}{(\lambda -\mu )(\mu -\nu )(\alpha -\gamma )}+\frac{(\lambda +\alpha )(\lambda +\gamma )(\mu +\alpha )(\nu +\beta )}{(\nu -\lambda )(\mu -\nu )(\alpha -\gamma )}+\frac{(\lambda +\beta )(\lambda +\gamma )(\mu +\alpha )(\nu +\alpha )}{(\lambda -\mu )(\nu -\lambda )(\alpha -\gamma )}]\mathrm{\Phi }(\lambda ,-\beta ,-\gamma ),\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill \xi (\lambda ,\mu ,\nu )=& \frac{(\lambda +\gamma )(\mu +\gamma )}{(\mu -\nu )(\gamma -\alpha )}[(\nu +\gamma )\mathrm{\Phi }(\lambda ,\mu ,-\gamma )-(\nu +\beta )\mathrm{\Phi }(\lambda ,-\beta ,\nu )-(\gamma -\beta )\mathrm{\Phi }(\lambda ,-\beta ,-\gamma )].\hfill \end{array}$$(15)

A simpler combination is obtained after substitution of Φ(λ, −β, ν) deduced from Eq. (13):

$$\begin{array}{cc}\hfill \psi (\lambda ,\mu ,\nu )=& \frac{(\lambda +\alpha )(\mu +\alpha )}{\alpha -\gamma }\mathrm{\Phi }(\lambda ,\mu ,\nu )-\frac{(\lambda +\alpha )(\mu +\gamma )}{\alpha -\gamma }\mathrm{\Phi }(\lambda ,\mu ,-\gamma )-(\lambda +\beta )\mathrm{\Phi }(\lambda ,-\beta ,-\gamma ),\hfill \end{array}$$(16)

$$\begin{array}{c}\hfill \xi (\lambda ,\mu ,\nu )=\frac{(\lambda +\gamma )(\mu +\gamma )}{\gamma -\alpha }[\mathrm{\Phi }(\lambda ,\mu ,\nu )-\mathrm{\Phi }(\lambda ,\mu ,-\gamma )].\end{array}$$(17)

We note that within the plane z  =  0, ξ(λ,μ,ν=−γ) is null. Moreover, along the x axis, ψ(λ,−β,−γ)  =  (−β+α)Φ(λ,−β,−γ), is null in the case of an axisymmetric potential with α = β.

Other expressions of the three integrals depending on Φ at three distinct positions can be obtained, we suppose under the necessary condition that the three chosen positions are on different sheets defined by the ellipsoidal coordinate system and including the (λ, μ, ν) position.

In the case of an axisymmetric potential with α = β, from Eq. (17) we recover the formula (exact for a Stäckel potential) used for approximate integrals for orbits within the potential of the Besançon Galaxy Model (Bienaymé et al. 2015, Eq. (A.7)).

$$\begin{array}{c}\hfill \xi (\lambda ,\nu )=(\lambda +\gamma )[\mathrm{\Phi }(\lambda ,\nu )-\mathrm{\Phi }(\lambda ,-\gamma )].\end{array}$$(18)

For an axisymmetric potential with setting α = β, we can, for instance, define another expression from Eq. (12), and express the third integral with the potential at two positions, one on the z axis and the other one on the z = 0 plane. It gives us

$$\begin{array}{c}\hfill \xi (\lambda ,\nu )=\frac{(\lambda +\gamma )}{(\lambda -\nu )}[(\nu +\gamma )\mathrm{\Phi }(\lambda ,-\gamma )-(\nu +\alpha )\mathrm{\Phi }(-\alpha ,\nu )]\end{array}$$(19)

This may be compared with the different variants of the Stäckel fudge for axisymmetric potentials proposed by Vasiliev (2019) and based on different sets of positions. Finally, we note that it is possible to write expressions depending on more than three positions.

3. Tests

We chose a Stäckel potential given by

$$\begin{array}{c}\hfill \mathrm{\Phi }=(\beta +\gamma )x^2+(\gamma +\alpha )y^2+(\alpha +\beta )z^2+{(-\alpha -\beta -\gamma +x^2+y^2+z^2)}^2\end{array}$$(20)

to verify the accuracy of our programs, as well as the validity of the various analytic expressions of the integrals of motion presented above. With α = −2, β = −1 and γ = 0, this potential has a unique minimum at x = y = z = 0 and all orbits are bound. With these values of the potential parameters, we examined the accuracy with which the integrals of motion I₂ (Eqs. (9) and (16)) and I₃ (Eqs. (10) and (17)) are numerically preserved. Over a period of 300 rotations around the centre, the energy of the orbits with an initial value of E = 1 is maintained to within 10⁻¹⁴. This accuracy is limited by rounding errors due to double-precision calculation. The other integrals I₂ and I₃ have values between about zero and one and each keep their initial value with a variation of less than 10⁻⁸, but greater than the accuracy that would result only from rounding errors.

We summarize our results and consider a potential more likely to represent a galaxy and similar to those chosen by Sanders & Binney (2014, 2015). The latter tested the conservation of quasi-integrals, quasi-actions, with an approximation also based on a Stäckel triaxial potential. This allows us to compare their results with ours.

We present our results obtained for the logarithmic potential:

$$\begin{array}{c}\hfill \mathrm{\Phi }(x,y,z)=\frac{1}{2}log(m_0^2+m^2),\end{array}$$(21)

where

$$\begin{array}{c}\hfill m^2=x^2+\frac{y^2}{q_y^2}+\frac{z^2}{q_z^2}\\ \hfill \end{array}$$

with m₀ = 0.3, q_y = 0.95, and q_z = 0.85. These flattening values q_y and q_z correspond to those used in the model of Sanders & Binney (2015). Note that these two values of the potential flattening q_Φ correspond to density flattening values of q_ρ ∼ 0.86 and 0.56 respectively. The commonly published values for flattening the dark halo of our Galaxy do not always agree with each other but still indicate that it is rather spherical. Thus, Posti & Helmi (2019) obtained q_ρ = 1.22 ± 0.23 from the analysis of globular cluster orbits while Malhan & Ibata (2019) found a flattening of the density q_ρ = 0.86 to explain the shape and velocity of the giant stellar current GD-1. In addition, Law & Majewski (2010) proposed a triaxial halo to explain the Sagittarius Galaxy stream.

From the potential of the Eq. (21), we calculated 5000 orbits of energy E = 1 with a Runge-Kutta of order 8 (Felberg 1968) over a period of about three hundred rotations around the centre. The initial position of the orbits is randomly selected on the sphere of radius r = 1 and the modulus of the initial velocity is then of the order of 1. The direction of the velocity is also randomly drawn. This choice makes it possible to cover many different types of orbits. We note that in the case of a logarithmic potential with m₀ = 0 (i.e. a scale free potential) this choice would cover all the different types of orbit. Here with m₀ = 0.3 we excluded the very inner orbits but still include many eccentric orbits, the most difficult to model accurately with our approximate integrals.

Once the orbital trajectories are determined, we adjusted the free parameters ${\mathrm{\Delta }}_1=\sqrt{\beta -\alpha }$, ${\mathrm{\Delta }}_2=\sqrt{\gamma -\beta }$ and ${\mathrm{\Delta }}_3=\sqrt{\gamma -\alpha }$ (positions of the foci of the ellipsoidal coordinate system) to obtain the best conservation of the quasi-integral I₂ and I₃ along all orbits. We obtained essentially the same adjusted focal points positions for this operation if we consider all orbits or if we simply limit ourselves to the three main closed orbits located in the three planes of symmetry.

We obtained the values Δ₁ = 0.49 and Δ₃ = 0.33, which give the quasi-integrals with the smallest fluctuations. The integrals of motion I₂ and I₃ then take values within the intervals [−0.28, 0.24], and [0.0.578], respectively. The relative variation, σ_I2 or σ_I3, of integrals (i.e. the variation divided by the amplitude of these intervals that is respectively 0.520 and 0.578 for I₂ and I₃) are shown in Fig. 1 in histogram form. The amplitude of this relative variations tells us the possibility of deciding when two orbits are distinct. For example, a relative variation of two per cent for I₂ and I₃ implies that the measurements of I₂ and I₃ make it possible to distinguish 50 × 50 orbits with the same energy. In addition, if we also have a two per cent accuracy on the energy, this would allow us to distinguish 125 000 separate orbits just by using kinematic criteria.

Fig. 1. Histograms for dispersions, σI2 and σI3, of quasi integrals for 5000 orbits with E = 1. The blue continuous line indicates dispersion for the integral I2. The red dashed line indicates dispersion for the integral I3.

Figure 1 shows that 45 per cent of the orbits have integrals with an accuracy better than two per cent, 68 per cent better than four per cent, and 91 per cent better than ten per cent. For the smallest values of σ_I2 and σ_I3, these dispersions are very strongly correlated. Figure 2 shows a correlation between the apocentre to pericentre ratio r_max/r_min of the orbit, and the accuracy of I₂ or I₃. This shows that the tube orbits, which are the least eccentric orbits, have the most accurate quasi-integrals Thus, orbits with apocenter to pericenter ratios of less than three have integrals with dispersions of less than three per cent. When we consider σ_I3, this correlation mainly concerns the tube orbits (Fig. 2).

Fig. 2. Left: σI2 uncertainties on the I2 determinations for 5000 orbits with E = 1 versus their apocentre/pericentre ratio. Right: σI3 uncertainties on the I3 determinations for 5000 orbits with E = 1 versus their apocentre/pericentre ratio.

It is interesting to compare these integrals with the angular momenta |L_z| and $L_{\perp }=\sqrt{L_x^2+L_y^2}$ which are frequently used as quasi-integrals (for example, see Helmi et al. 2017). These last quantities are easily computed and do not require the integration of orbits. These angular momenta are approximate integrals that remain tolerable in the case of an almost spherical potential. For the potential considered here, which is close to spherical, they are however significantly less accurate than Stäckel approximations. Figure 3 shows a comparison of the histograms of the relative accuracies of these different quasi-integrals. The difference is most noticeable for the least eccentric orbits, which are also the most accurate. For these orbits, we note an improvement by a factor of five (I₂) and of three (I₃) for the accuracy gain. The same gain is obtained on the position of the peak of the I₂ and I₃ histograms, also in favour of the Stäckel approximations.

Fig. 3. Comparison of accuracies achieved using Stäckel (blue continuous line) or angular momentum (black dashed line) approximations. Left: histograms of σI2 and σLz. Right: histograms of σI3 and σL⊥.

Figure 4 shows the correlation between the accuracies σ_I3 and σ_L⊥. For σ_I2 and σ_Lz the correlation also exists but mainly for tube orbits while box orbits have lower accuracies (Fig. 4). Figure 5a shows the distribution of the three main types of orbits in the (I₂, I₃) plane: box orbits and tube orbits, these latter ones classified according to their orientation along the major axis or the minor axis of symmetry of the potential. Figure 5a can be compared to Figs. 17 and 18 in de Zeeuw (1985a), which classifies the different families of orbits (see Fig. 6). His figures show a three-dimensional classification of orbits in terms of E, I₂, and I₃ for bound orbits of a triaxial Stäckel potential and a cross-section for a given energy E. Figures 5b and c show the distribution of orbits with their accuracy colour-coded.

Fig. 4. Comparison of accuracies. Left: σL⊥ versus σI3 for 5000 orbits. Right: σLz versus σI2 for 5000 orbits.

Fig. 5. Two-dimensional classification of orbits. Each point (i2, i3) corresponds to an orbit. Left: different types of orbits are the box orbits (blue), the outer long axis tube orbits (red), and the short axis tube orbits (green). Centre: the σI2 dispersion is colour-coded. Right: the σI3 dispersion is colour-coded.

Fig. 6. From top to bottom: short-axis tube orbit, long-axis tube orbit, and box orbit.

Figure 7 shows a section of the phase space for two series of orbits, each figure corresponding to orbits in the (x, y) and (y, z) planes of symmetry of the potential. Contours in sections are calculated numerically from the orbital trajectories on the one hand, and on the other hand they are obtained using our analytic expressions of the approximate integrals. The tube orbits close to the central periodic tube orbit are best represented by the analytic expression since the parameters Δ₁ and Δ₃ of the Stäckel potentials have been adjusted on these orbits.

Fig. 7. Poincaré sections. Red: orbit computation. Blue: from the analytic quasi-integrals. Left: (x, vx) section for orbits within the (x, y) plane, (z = ż = 0). Right: (y, vy) section for orbits within the (y, z) plane, x = ẋ = 0.

It must be noted that we calculated the variations of I₂ and I₃ over a long time interval, several thousand rotations around the centre. If we were interested in studying Galactic stellar streams, for example, these would have short extension along their orbit, and the dispersions of I₂ or I₃ along that extension would be significantly smaller. Figure 8 shows that abrupt variations of the quasi-integrals occur during the transit close to the pericentre where the potential variation is the most rapid. This is very sensitive for box orbits that show the greatest variations, but between two transits, away from the pericentre, the quasi-integrals are significantly more stationary.

Fig. 8. Variation of I2 (blue) and I3 (red) for three orbits. Over the short time-length of the order of an average rotation period (Δt  =  3), the variations of the quasi-integrals are significantly smaller compared to the variations over long periods. The largest variations of the quasi-integrals are observed for I3(∼ − 0.2) (bottom curve) of the box orbit (I2 ∼ 0). The sharp variations of I3 are related to the passage at the pericentre. The short-axis tube (I2  ∼  I3  ∼  0.1) and the long-axis tube (I2 ∼ −0.1 and I3  ∼  0.5) orbits have smoother variations of the quasi-integrals.

In Fig. 7 we note the lack of visible resonant orbits: either they are actually absent or they occupy a small space. To conclude, it can be expected that some simple modifications of the analytic expressions of I₂ and I₃ should lead to better agreement between the analytic and numerical phase space sections (Fig. 7). A difficulty to obtaining better analytic quasi-integrals is to extend this modification to the entire volume of the phase space and not just to a few sections of the phase space. Such a modification with additional parameters must however exist as it has been obtained in Bienaymé & Traven (2013), but in this latter case using a considerable number of free parameters to build approximate integrals.

4. Distribution functions and Jeans equations

The developments described above to model the gravitational dynamics can have several practical applications. An initial application is the search for stellar streams in the Galactic halo, which can be carried out by identifying the membership of stars in nearby orbits Malhan et al. (2018) or by identifying concentrations in the space of the integrals of motion. The quasi-integrals developed here allow us to take into consideration the probable triaxiality of the Galactic halo (Law & Majewski 2010). Conversely and simultaneously, this makes it possible to determine the Galactic potential (Malhan & Ibata 2019). Another important application is the construction of distribution functions dependent on three integrals of motion. Such distribution functions make it possible to calculate the velocity moments and density of a stationary stellar population. In addition to describing stellar populations, they can also be used to measure the gravitational potential using the Jeans equations in the case of triaxial potentials.

We test here the accuracy of a triaxial distribution function by applying the Jeans equations to find the potential used to build this distribution function. This test has already been used to validate the stellar disc distribution functions (generalisation of Shu distribution, Shu 1969) for the Besançon Galaxy model (Bienaymé et al. 2015). This last model is axisymmetric and in addition to the energy and kinetic moment, we used a third approximate integral of the Stäckel type, a particular case of the study presented here. Sanders & Binney (2014) also used the Jeans equations as tests of their approximate integrals for axisymmetric models for non-axisymmetric models.

A second possible application of the two integrals I₂ and I₃ is the generalisation of Shu’s distribution functions for stellar discs with slight ellipticity. In this case, in Shu’s three-dimensional axisymmetric version (Eqs. (8)–(9) in Bienaymé et al. 2015), it is sufficient to replace the L_z angular momentum with its generalisation I₂, and to replace the E_circ energy of the circular orbit having the L_z angular momentum with the energy of the main closed orbit having the value I₂ as a second integral.

Another application, the one we develop here, is the use of integrals to build a triaxial distribution function that models the halo stellar distribution.

4.1. Prescribed halo distribution functions

4.1.1. f(E, L²)

In order to test the relevance of the quasi-integrals through the Jeans equations, we first constructed a distribution function for the stellar halo in the particular case of an axisymmetric logarithmic potential. This distribution function for a spherical potential is inspired by the stationary distribution:

$$\begin{array}{c}\hfill f(E,L)=exp[-(E+L^2/2{\alpha }^2)].\end{array}$$(22)

Osipkov (1979) and Merritt (1985) showed that for a spherical potential, any distribution function of the form f(E + L²/2α²) has an isotropic velocity distribution at the centre and is radially anisotropic at a large distance from the centre, (σ_r does not depend on r and α is the transition radius from the inner isotropic core to the outer regions) by following the relationship

$$\begin{array}{c}\hfill \frac{{\sigma }_{\mathrm{r}}^2}{{\sigma }_{\mathrm{t}}^2}=\frac{{\sigma }_{\mathrm{r}}^2}{\frac{1}{2}<v_{\mathrm{t}}^2>}=1+\frac{r^2}{{\alpha }^2}·\end{array}$$(23)

The ratio of radial-to-tangential velocity dispersions changes continuously with distance to the centre, and the parameter α defines the transition radius. In the specific case of a logarithmic potential, this scale factor α can simply be made ineffective by replacing it in the distribution function by an amount that is both energy dependent and proportional to the r radius. Thus, with the potential

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=v_0^{2}logr,\end{array}$$(24)

the radius of the circular orbit r_circ of energy E is given by

$$\begin{array}{c}\hfill E=v_0^{2}log(r_{\mathrm{circ}})+1/2v_0^2,\end{array}$$(25)

and we deduced the following quantity related to it

$$\begin{array}{c}\hfill r_{\mathrm{c}}(E)=exp[E/v_0^2-1/2]=rexp[\frac{1}{2}\frac{v_{\mathrm{r}}^2+v_{\mathrm{t}}^2}{v_0^2}-\frac{1}{2}].\end{array}$$(26)

We then set for the distribution function:

$$\begin{array}{c}\hfill f(E,L)=exp[-\frac{n}{v_0^2}\left(E+a\frac{L^2}{r_{\mathrm{c}}{(E)}^2}\right)],\end{array}$$(27)

with only two free parameters, n and a. We can see that with the potential being logarithmic, the term on the right-hand side separates into two terms. A first term contains only the variable r, which gives for the associated density law: ν(r)∝r⁻ⁿ. The second term contains only the velocity components, and therefore the velocity moments are independent of the r radius. It should be noted that in the case of the isothermal sphere, a = 0, the velocity dispersion is isotropic and the radial and tangential velocity dispersions are ${\sigma }_{\mathrm{r}}={\sigma }_{\mathrm{t}}/\sqrt{2}=v_0/\sqrt{n}$.

Table 1 gives the velocity dispersions for some values of the a parameter when v₀ = 1. The velocity distribution is isotropic for a = 0, and with a radial or tangential bias depending on whether a >  0 or a <  0 The cases a = −∞ or a = +∞ correspond to distributions with purely tangential or purely radial orbits, respectively. The tangential velocity has two components, according to θ and ϕ. |v_t| refers to the tangential velocity modulus, while ${\sigma }_{\mathrm{t}}/\sqrt{2}$ refers to the dispersion of a single tangential component. Figure 9 represents the corresponding f(v_r, v_t) distributions.

Fig. 9. Distribution function f(vr, vt) from top-left to bottom-right, with a = −2, 1, 0, +1, +2 from tangentially to radially anisotropic velocity distributions.

4.1.2. f(E, L_z, L_⊥)

Still within the framework of a spherical and logarithmic potential, the previous distribution function is modified in order to have a triaxial velocity ellipsoid, with the three velocity dispersions σ_r, σ_ϕ, and σ_θ, different two-by-two. To do this, we introduced a dependency on the angular momentum L_z and we wrote:

$$\begin{array}{c}\hfill f=exp\{-n[E+\frac{1}{r_{\mathrm{E}}^2}(b{L}_z^2+c(L^2-L_z^2))]\}.\end{array}$$(28)

As a result, the associated density distribution is no longer necessarily spherical, but may also be oblate or prolate with the z axis as the axis of rotational symmetry. The density distribution remains proportional to r⁻ⁿ, but is now dependent on the angle θ (the angle between the vertical axis z and the vector radius r, and ϕ is the azimuth angle). The velocity distribution does not depend on the distance to the centre r but only on the angle θ.

The parameters b and c allow modification of the values of the velocity dispersion ratios along the three main axes of the distribution f(v_r, v_ϕ, v_θ). For b = c, the distribution function has a spherical density (see Eq. (28)). In these cases, the velocity dispersions σ_ϕ and σ_θ are equal and do not depend on r, ϕ, or θ (Tables 2 and 3 in the second row). For the three special cases b = 0, c = 0, or b = c, in the z = 0 plane, two of the three components of the velocity dispersions are equal (see Tables 2 and 3). The corresponding densities are shown in Fig. 10 (also the first and third columns). It should be noted that if b <  c, in the z = 0 plane the vertical velocity dispersion σ_θ being smaller than in the spherical case, the density distribution is oblate. And conversely for b >  c, the density distribution is prolate or less oblate. The corresponding values of the velocity dispersion along the three main axes of the velocity distributions are given in Table 3. More generally, this distribution function, (see Eq. (28)), allows us to model the velocity dispersions with σ_ϕ/σ_r and σ_θ/σ_r ratios arbitrarily chosen within the z = 0 plane. Outside the z = 0 plane towards the poles θ = 0 or π/2, the velocity distributions along the three main axes vary, and, for symmetry reasons, the dispersions σ_θ and σ_ϕ become equal to the poles (Table 4).

Fig. 10. Density iso-contours (logarithmically equally spaced) for distribution function (Eq. (28)) for b and c values of Table 2 (raw and column corresponding).

4.1.3. f(E, I₂, I₃)

Finally, we considered the triaxial logarithmic potential (see Eq. (21)) with m₀ = 0.3, q_y = 0.95, and q_z = 0.85 and we defined the distribution function (with n = 3):

$$\begin{array}{c}\hfill f=exp[-n(E+\frac{1}{r_{\mathrm{E}}^2}(b{I}_2+c{I}_3))].\end{array}$$(29)

For a triaxially symmetrical potential, the quasi-integrals I₂ and I₃ are workable generalisations of $L_z^2/2$ and (L² − Lz²)/2. The figures presented here are made using for I₂ and I₃ the expressions given by Eqs. (16) and (17).

In the particular case where b = 0 and c = 0, the distribution function depends only on energy and is therefore an exact solution of the collisionless Boltzmann equation. When the b or c parameters are non-zero, the potential is no longer a Stäckel potential. The integrals I₂ and I₃ are approximated and the distribution function is therefore not strictly stationary. Moreover, the presence of a non-zero core radius, m₀ = 0.3, means that the potential is no longer strictly scale-free. As a result, the tensor of velocity dispersions varies: it is approximately isotropic in the centre and anisotropic at large radii. The associated density has a core radius of ∼0.3 and at larger radii it decreases as r⁻³.

The distribution function (see Eq. (29)) is relatively accurate in an extended volume of space and only for small values in modulus of the constants b and c. We describe two cases corresponding to oblate distributions of the tracer: (b, c) = (+0.3, +0.3) and ( − 0.2, −0.2). In both cases, the dispersions are almost isotropic near the centre. Close to the galactic plane, z = 0, the tangential velocity dispersions are similar σ_ϕ ≃ σ_θ. Figure 11 shows the three components of the velocity dispersions versus the distance x (with y = z = 0) to the centre. At large radii, for the first case the distribution function is dominated by radial orbits, in the second case by tangential orbits. These two examples correspond to triaxial density distributions of the stellar tracer with axis ratios q_z = 0.90 and 0.77, respectively. Figure 12 shows the iso-contours of the latter model that is dominated with radial orbits.

Fig. 11. Velocity dispersions versus x-distance (y = z = 0) to centre for models given by Eq. (29). The blue lines model (b, c)=(+0.3, +0.3) is dominated by radial obits. The red lines models (b, c)=(−0.2, −0.2) are dominated by tangential obits. Continuous lines indicate radial dispersions σR, while short- and long-dashed lines indicate tangential dispersions σϕ and σθ, respectively.

Fig. 12. Density for radially biased model (b = c = +0.3). Shown above are equally spaced iso-contours of the logarithm of the density.

The degree of accuracy of these distribution functions is evaluated in the following paragraph by applying the Jeans equations to verify that the forces are recovered with sufficient accuracy. The main element that limits the accuracy of the distribution function is the modelling of box orbits, for which the quasi-integral orbits are the least well-preserved along the orbits (see also for example the Poincaré sections that are less accurate for box orbits than for tube orbits in Fig. 7).

4.2. Test with Jeans equations

The previous distribution function must satisfy the collisionless Boltzmann equation and therefore the Jeans equations that we rewrote as:

$$\begin{array}{c}\hfill {\displaystyle \sum _i\frac{\partial <v_i{v}_j>}{\partial x_i}+\sum _i<v_i{v}_j>\frac{\partial log\nu }{\partial x_i}+\frac{\partial \mathrm{\Phi }}{\partial x_j}=0(j=1,2,3).}\end{array}$$(30)

These equations allow us to calculate the force field associated with the gravitational potential Φ. To do this, it is sufficient to measure the density and velocity tensor of a stellar tracer represented here by the distribution function Eq. (29). As the associated density decreases rapidly according to a power law, we used the Jeans equations in the form of Eq. (30) to obtain better numerical accuracy. We applied these equations for the potential Φ (Eq. (21)) with m₀ = 0.3. This allows the accuracy of the distribution function to be tested and the usefulness of the quasi-integral I₂ and I₃ to be estimated.

The results obtained can be compared with those established in the case of axisymmetric potentials with approximate integrals Sanders (2012) and in particular with the potential of the Besançon Galaxy model (Bienaymé et al. 2015). While for axisymmetric potentials and distribution functions modelling tube orbits, we achieved accuracies of the order of a thousandth or hundredth over a large volume of the model, here we obtained accuracies of the order of a few per cent over a volume, which, although extended, is comparatively more limited. These results are similar to the work of Sanders & Binney (2014) in many respects. They use the actions as approximate integrals of motion which are also constructed by approximating orbits using triaxial Stäckel potentials. They also find a more limited accuracy than in axisymmetric cases. The main reason for this is that the modelling of box orbits passing near the centre is less accurate than that of tube orbits that remain far from the centre.

For the model with predominantly radial velocity dispersions (b = +0.3, c = +0.3), that is, one dominated with box orbits and that also best corresponds to our knowledge of the velocity distribution of galactic halo stars (Bland-Hawthorn & Gerhard 2016), we find the force field with an accuracy of a few per cent within the volume of radius r = 1. This is to be compared to the core radius m₀ = 0.3, and up to ten per cent at distances r = 2 from the centre. The panels of Fig. 13 show for the three force components and along three axes, the differences between the potential gradient ∇Φ, and the forces recovered from the Jeans equations. In addition, near the (x, y) plane, if |z|< 1, the distribution function allows us to accurately find the force field to within a few per cent and this for r radii up to ∼40.

Fig. 13. Accuracy of Jeans equations for model with b = c = −0.2. Top: three forces Fx, Fy, Fz along the three axes (x, y = 1, z = 1), (x, y = 1, z = x), and (x, y = x/2, z = x/3). Bottom: relative error for the same axes and forces.

The model with the dominant tangential velocity dispersion (b = −0.2, c = −0.2) allows us to find the force field on a much larger volume going at least to r = 50. Figure 14 shows that the relative force errors reach a plateau at r = 20. Since the distribution function is dominated by tangential orbits (essentially tube or centrophobic orbits) whose quasi-integral orbits are better preserved, the Jeans equations are therefore satisfied with greater precision.

Fig. 14. Accuracy of Jeans equations for model with b = c = +0.3. Top: three forces Fx, Fy, Fz along the three axes (x, y = 1, z = 1), (x, y = 1, z = x), and (x, y = x/2, z = x/3). Bottom: relative error for the same axes and forces.

A second and complementary test is the comparison of the total mass of the model, ρ_total, as given by the Poisson equation applied directly to the potential, and that deduced by using the forces obtained by the Jeans equations (see Figs. 15 and 16). For the radial, velocity-dominant model and for r <  1, the density ρ_total is found with systematic deviations, varying according to the position with an average of 6%, and with a dispersion of 5%. For r <  2 the density ρ_total is also found with systematic deviations with a mean systematic of 6% but with a greater dispersion of 9% depending on position. For the tangential, velocity-dominant model, the results are quite similar, with a relative error on the total mass density that depends on position and varies from about zero in the centre to 9% towards r = 2. Globally, we have better accuracy in finding the total mass density using the distribution function dominated by tangential orbits.

Fig. 15. Accuracy of the Jeans equation for model b = c = −0.2 to recover total mass density ρtotal from Poisson equation. Left: ρtotal along the axis (x, y − 0, z = 0). Right: relative error on ρtotal.

Fig. 16. Accuracy of the Jeans equation for model b = c = +0.3 to recover total mass density ρtotal from Poisson equation. Left: ρtotal along the axis (x, y − 0, z = 0). Right: relative error on ρtotal.

5. Conclusion

We model the orbits and distribution functions of stellar populations in a triaxial potential. This is done using an approximation of orbits in the frame of Stäckel potentials. Kent & de Zeeuw (1991) already noted that an accurate approximation of a quasi-invariant integral can be achieved by adjusting a Stäckel potential for each of the orbits of a galactic potential. Here, we show that quasi-integrals can be written explicitly as potential dependent, by generalising the usual formulation of Stäckel potential integrals. These analytic expressions of integrals can be used to plot sections in phase space and by a simple quadrature allow the derivation quasi-actions. They can also be used to build distribution functions. These integrals allow us to model four different types of orbits in a triaxial potential, and they are better preserved than the angular moment components, which are sometimes used as approximate integrals. We note that tube orbits are modelled with very high accuracy, while the modelling of box orbits passing near the galactic centre are less accurate. The use of quasi-integrals allows us to build distribution functions for a triaxial potential with a triaxial velocity moment tensor. The distribution functions, that are dominated by tube orbits (i.e. dispersions of velocities dominated by tangential components), are more precise than those dominated by radial velocities. We use the Jeans equations and find the components of the force of the gravity field with good accuracy over an extended volume. On the other hand, the total mass density generating the potential and the gravity field is recovered with a significantly lower accuracy. It is indeed necessary to calculate second derivatives of the model distribution, which would require an accuracy on the distribution function that is better by one to two orders of magnitude.

We have already performed a similar analysis for an axisymmetric version of the Besançon Galaxy Model (Bienaymé et al. 2015). By generalising this work to triaxially symmetrical potentials, we improve the simultaneous modelling of four families of tube and box orbits. The latter have a more complex structure and occupy a very large volume in the configuration space and are modelled with less precision than tube orbits. Our results can be compared to those of McMillan & Binney (2008), Sanders (2012), Sanders & Binney (2014), Posti et al. (2015) on axisymmetric or triaxially symmetrical potentials. Their analyses are also based on the determination of quasi-integrals, which in their case are actions that they used to define the distribution functions. The accuracies they obtain are of the same order of magnitude as ours. It should be noted that in our case the quasi-integrals are analytic and easier to calculate. As for the actions, they can be evaluated here by a simple quadrature using the expressions of quasi-integrals employed, for example, to plot the different sections of the phase space (Fig. 7). Our method and the Stäckel fudge (Sanders & Binney 2014) are based on the same principles and are less accurate compared to the torus fitting (McGill & Binney 1990, and see Sanders & Binney 2016, for a review of different methods). These methods are included in numerical tool boxes proposed by Bovy (2015) and Vasiliev (2019). It is then necessary to test their accuracy for each new potential and for each distribution function considered. Their accuracy remains limited by the accuracy of the modelled orbits, especially the box orbits. We have shown through an example that if the modelling of integrals and distribution functions made it possible to recover the forces of the gravity field with an accuracy of a few per cent over a large volume of space, the same is not true when it comes to measuring the total mass density of the model, meaning the one that generates the potential.

The advantage of the method presented here is its numerical speed. However, more precise methods such as torus fitting are necessary if high precision is desired to accurately recover the potential and the total mass density of a model. It should be noted that this last method does not make it possible to simultaneously model all the resonances, and that resonances can be numerous for box orbits, for example, in the case of a significantly flattened logarithmic potential (Bienaymé & Traven 2013). Finally, we note that quasi-integral variations are significant when orbits pass near the Galactic centre, while they remain an order of magnitude more constant between two passages at the Galactic pericentre. This makes it possible to consider using these integrals for the identification of stellar streams that have a limited extension along an orbit.

The work presented here suggests possible improvements. Figure 7 presents two sections of the phase space and directly shows the limitation of the quasi-integrals developed here. Finding more precise integrals can then be reduced to the question of finding a more general analytic expression, depending on a few additional parameters in order to reproduce more exactly the shape of the iso-contours of the sections of the phase space. Another direction for improvement is the second family of Stäckel potentials, which has been studied very little and whose equations of motion are separable in a parabolic coordinate system (Ollongren 1962; Lynden-Bell 1962; Tsiganov 2012). These potentials also depend on a generic function and therefore have a certain generality that merits further study.

Acknowledgments

The author would like to thank the International Space Science Institute, Berne, Switzerland, for providing financial support and meeting facilities.

References

Appendix A: (x,y,z) (λ, μ, ν) transformation

Stäckel potentials are defined using ellipsoidal coordinates. The (λ, μ, ν) coordinates define confocal ellipsoids and hyperboloids and are roots of the equation

$$\begin{array}{c}\hfill \frac{x^2}{\tau +\alpha }+\frac{y^2}{\tau +\beta }+\frac{z^2}{\tau +\gamma }=1,\end{array}$$(A.1)

with α, β, and γ, which are three fixed numbers with −γ ≤ ν ≤ −β ≤ μ ≤ −α ≤ λ.

The expressions for x, y, and z Cartesian coordinates are

$$\begin{array}{cc}\hfill x^2& \hfill =\frac{(\lambda +\alpha )(\mu +\alpha )(\nu +\alpha )}{(\alpha -\beta )(\alpha -\gamma )},\hfill \\ \hfill y^2& \hfill =\frac{(\lambda +\beta )(\mu +\beta )(\nu +\beta )}{(\beta -\gamma )(\beta -\alpha )},\hfill \\ \hfill z^2& \hfill =\frac{(\lambda +\gamma )(\mu +\gamma )(\nu +\gamma )}{(\gamma -\alpha )(\gamma -\beta )}·\hfill \end{array}$$(A.2)

The expressions for λ, μ, and ν coordinates are obtained from the relations

$$\begin{array}{cc}\hfill -b& \hfill =\lambda +\mu +\nu =-\alpha -\beta -\gamma +x^2+y^2+z^2,\hfill \\ \hfill c& \hfill =\lambda \mu +\mu \nu +\nu \lambda =\alpha \beta +\beta \gamma +\gamma \alpha -(\beta +\gamma )x^2-(\gamma +\alpha )y^2-(\alpha +\beta )z^2,\hfill \\ \hfill -d& \hfill =\lambda \mu \nu =-\alpha \beta \gamma +(\beta \gamma )x^2+(\gamma \alpha )y^2+(\alpha \beta )z^2,\hfill \end{array}$$(A.3)

and they are the three solutions of the third-degree polynomial equation

$$\begin{array}{c}\hfill X^3+b{X}^2+cX+d=0,\end{array}$$(A.4)

or

$$\begin{array}{c}\hfill Y^3+pY+q=Y^3+(-\frac{1}{3}{b}^2+c)Y+(\frac{2}{27}{b}^3+d-\frac{1}{3}cb)=0,\end{array}$$(A.5)

with X = Y − b/3. The equation has three solutions given that 4p³ + 27q² <  0.

The solutions are $Y_1=rsin\frac{\phi }{3}$, $Y_2=rsin(\frac{\phi }{3}+\frac{2\pi }{3})$, $Y_3=rsin(\frac{\phi }{3}+\frac{4\pi }{3})$, with $r=2\sqrt{-\frac{p}{3}}$ and $sin\phi =-\frac{3q}{\mathit{pr}}·$

For spherical potentials, we have α = β = γ and r² = x² + y² + z² = λ + α.

For prolate potentials, (λ, μ) define oblate spheroidal coordinates and β = γ. The x axis is the axis of axisymmetry and we have ${\stackrel{\sim }{z}}^2=y^2+z^2$, $\lambda +\mu =-\alpha -\gamma +x^2+{\stackrel{\sim }{z}}^2$, and $\lambda \mu =\alpha \gamma -\gamma x^2-\alpha {\stackrel{\sim }{z}}^2$.

For oblate potentials, (λ, ν) define prolate spheroidal coordinates and α = β. The z axis is the axis of axisymmetry and we have R² = x² + y², λ + ν = −α − γ + R² + z², and λν = αγ − γR² − αz².

All Tables

All Figures

	Fig. 1. Histograms for dispersions, σI2 and σI3, of quasi integrals for 5000 orbits with E = 1. The blue continuous line indicates dispersion for the integral I2. The red dashed line indicates dispersion for the integral I3.
In the text

	Fig. 2. Left: σI2 uncertainties on the I2 determinations for 5000 orbits with E = 1 versus their apocentre/pericentre ratio. Right: σI3 uncertainties on the I3 determinations for 5000 orbits with E = 1 versus their apocentre/pericentre ratio.
In the text

	Fig. 3. Comparison of accuracies achieved using Stäckel (blue continuous line) or angular momentum (black dashed line) approximations. Left: histograms of σI2 and σLz. Right: histograms of σI3 and σL⊥.
In the text

	Fig. 4. Comparison of accuracies. Left: σL⊥ versus σI3 for 5000 orbits. Right: σLz versus σI2 for 5000 orbits.
In the text

	Fig. 5. Two-dimensional classification of orbits. Each point (i2, i3) corresponds to an orbit. Left: different types of orbits are the box orbits (blue), the outer long axis tube orbits (red), and the short axis tube orbits (green). Centre: the σI2 dispersion is colour-coded. Right: the σI3 dispersion is colour-coded.
In the text

	Fig. 6. From top to bottom: short-axis tube orbit, long-axis tube orbit, and box orbit.
In the text

	Fig. 7. Poincaré sections. Red: orbit computation. Blue: from the analytic quasi-integrals. Left: (x, vx) section for orbits within the (x, y) plane, (z = ż = 0). Right: (y, vy) section for orbits within the (y, z) plane, x = ẋ = 0.
In the text

Fig. 8. Variation of I2 (blue) and I3 (red) for three orbits. Over the short time-length of the order of an average rotation period (Δt  =  3), the variations of the quasi-integrals are significantly smaller compared to the variations over long periods. The largest variations of the quasi-integrals are observed for I3(∼ − 0.2) (bottom curve) of the box orbit (I2 ∼ 0). The sharp variations of I3 are related to the passage at the pericentre. The short-axis tube (I2  ∼  I3  ∼  0.1) and the long-axis tube (I2 ∼ −0.1 and I3  ∼  0.5) orbits have smoother variations of the quasi-integrals.

In the text

	Fig. 9. Distribution function f(vr, vt) from top-left to bottom-right, with a = −2, 1, 0, +1, +2 from tangentially to radially anisotropic velocity distributions.
In the text

	Fig. 10. Density iso-contours (logarithmically equally spaced) for distribution function (Eq. (28)) for b and c values of Table 2 (raw and column corresponding).
In the text

	Fig. 11. Velocity dispersions versus x-distance (y = z = 0) to centre for models given by Eq. (29). The blue lines model (b, c)=(+0.3, +0.3) is dominated by radial obits. The red lines models (b, c)=(−0.2, −0.2) are dominated by tangential obits. Continuous lines indicate radial dispersions σR, while short- and long-dashed lines indicate tangential dispersions σϕ and σθ, respectively.
In the text

	Fig. 12. Density for radially biased model (b = c = +0.3). Shown above are equally spaced iso-contours of the logarithm of the density.
In the text

	Fig. 13. Accuracy of Jeans equations for model with b = c = −0.2. Top: three forces Fx, Fy, Fz along the three axes (x, y = 1, z = 1), (x, y = 1, z = x), and (x, y = x/2, z = x/3). Bottom: relative error for the same axes and forces.
In the text

	Fig. 14. Accuracy of Jeans equations for model with b = c = +0.3. Top: three forces Fx, Fy, Fz along the three axes (x, y = 1, z = 1), (x, y = 1, z = x), and (x, y = x/2, z = x/3). Bottom: relative error for the same axes and forces.
In the text

	Fig. 15. Accuracy of the Jeans equation for model b = c = −0.2 to recover total mass density ρtotal from Poisson equation. Left: ρtotal along the axis (x, y − 0, z = 0). Right: relative error on ρtotal.
In the text

	Fig. 16. Accuracy of the Jeans equation for model b = c = +0.3 to recover total mass density ρtotal from Poisson equation. Left: ρtotal along the axis (x, y − 0, z = 0). Right: relative error on ρtotal.
In the text