Issue 
A&A
Volume 627, July 2019



Article Number  A123  
Number of page(s)  15  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201935573  
Published online  11 July 2019 
Integrals of motion for nonaxisymmetric potentials
Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, 11 rue de l’Université, 67000 Strasbourg, France
email: olivier.bienayme@unistra.fr
Received:
29
March
2019
Accepted:
22
May
2019
Context. The modelling of stationary galactic stellar populations can be performed using distribution functions.
Aims. This paper aims to write explicit integrals of motion and distribution functions.
Methods. We propose an analytic formulation of the integrals of motion with an explicit dependence on potential. This formulation applies to potentials with rotational symmetry or triaxial symmetry. It is exact for Stäckel potentials and approximate for other potentials.
Results. Modelling a stationary stellar population using these integrals of motion allows the force field to be found with satisfactory accuracy. On the other hand, the mass density distribution that generates the force field and the gravitational potential is recovered with less accuracy due to lower precision in modelling boxtype orbits.
Key words: Galaxy: kinematics and dynamics / methods: numerical
© O. Bienaymé 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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^{1}. 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 angleaction 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 nonStä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; LyndenBell 1962; de Zeeuw & LyndenBell 1985; de Zeeuw 1985a,b). Stäckel potentials can be written in ellipsoidal coordinates (see Appendix), depending on three arbitrary functions ζ(λ), η(μ), and θ(ν):
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 ζ(λ):
Similarly, η(μ) and θ(ν) are obtained evaluating Φ(−α,−β,ν) and Φ(−α,μ,−γ):
Then, substituting these expressions of ζ, η, and θ within the expression of Φ (Eq. (1)), we obtain
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 & LyndenBell (1985, see Eqs. (7)–(9)). They also detail the integrals I_{2} and I_{3} 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 & LyndenBell (1985, Eqs. (7) and (14)), and for convenience write the exact integrals in compact form as
and
According to de Zeeuw & LyndenBell (1985), who give a complete discussion, “I_{2} and I_{3} 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, . We note that the third integral is null, I_{3} = 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):
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):
A simpler combination is obtained after substitution of Φ(λ, −β, ν) deduced from Eq. (13):
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)).
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
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
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_{2} (Eqs. (9) and (16)) and I_{3} (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^{−14}. This accuracy is limited by rounding errors due to doubleprecision calculation. The other integrals I_{2} and I_{3} have values between about zero and one and each keep their initial value with a variation of less than 10^{−8}, 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 quasiintegrals, quasiactions, 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:
where
with m_{0} = 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 GD1. 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 RungeKutta 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} = 0 (i.e. a scale free potential) this choice would cover all the different types of orbit. Here with m_{0} = 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 , and (positions of the foci of the ellipsoidal coordinate system) to obtain the best conservation of the quasiintegral I_{2} and I_{3} 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 Δ_{1} = 0.49 and Δ_{3} = 0.33, which give the quasiintegrals with the smallest fluctuations. The integrals of motion I_{2} and I_{3} 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_{2} and I_{3}) 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_{2} and I_{3} implies that the measurements of I_{2} and I_{3} 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 I_{2}. The red dashed line indicates dispersion for the integral I_{3}. 
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_{2} or I_{3}. This shows that the tube orbits, which are the least eccentric orbits, have the most accurate quasiintegrals 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 I_{2} determinations for 5000 orbits with E = 1 versus their apocentre/pericentre ratio. Right: σ_{I3} uncertainties on the I_{3} 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 which are frequently used as quasiintegrals (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 quasiintegrals. 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_{2}) and of three (I_{3}) for the accuracy gain. The same gain is obtained on the position of the peak of the I_{2} and I_{3} 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_{2}, I_{3}) 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 threedimensional classification of orbits in terms of E, I_{2}, and I_{3} for bound orbits of a triaxial Stäckel potential and a crosssection for a given energy E. Figures 5b and c show the distribution of orbits with their accuracy colourcoded.
Fig. 4. Comparison of accuracies. Left: σ_{L⊥} versus σ_{I3} for 5000 orbits. Right: σ_{Lz} versus σ_{I2} for 5000 orbits. 
Fig. 5. Twodimensional classification of orbits. Each point (i_{2}, i_{3}) 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 colourcoded. Right: the σ_{I3} dispersion is colourcoded. 
Fig. 6. From top to bottom: shortaxis tube orbit, longaxis 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 Δ_{1} and Δ_{3} of the Stäckel potentials have been adjusted on these orbits.
Fig. 7. Poincaré sections. Red: orbit computation. Blue: from the analytic quasiintegrals. Left: (x, v_{x}) section for orbits within the (x, y) plane, (z = ż = 0). Right: (y, v_{y}) section for orbits within the (y, z) plane, x = ẋ = 0. 
It must be noted that we calculated the variations of I_{2} and I_{3} 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_{2} or I_{3} along that extension would be significantly smaller. Figure 8 shows that abrupt variations of the quasiintegrals 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 quasiintegrals are significantly more stationary.
Fig. 8. Variation of I_{2} (blue) and I_{3} (red) for three orbits. Over the short timelength of the order of an average rotation period (Δt = 3), the variations of the quasiintegrals are significantly smaller compared to the variations over long periods. The largest variations of the quasiintegrals are observed for I_{3}(∼ − 0.2) (bottom curve) of the box orbit (I_{2} ∼ 0). The sharp variations of I_{3} are related to the passage at the pericentre. The shortaxis tube (I_{2} ∼ I_{3} ∼ 0.1) and the longaxis tube (I_{2} ∼ −0.1 and I_{3} ∼ 0.5) orbits have smoother variations of the quasiintegrals. 
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_{2} and I_{3} should lead to better agreement between the analytic and numerical phase space sections (Fig. 7). A difficulty to obtaining better analytic quasiintegrals 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 quasiintegrals 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 nonaxisymmetric models.
A second possible application of the two integrals I_{2} and I_{3} is the generalisation of Shu’s distribution functions for stellar discs with slight ellipticity. In this case, in Shu’s threedimensional axisymmetric version (Eqs. (8)–(9) in Bienaymé et al. 2015), it is sufficient to replace the L_{z} angular momentum with its generalisation I_{2}, 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_{2} 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^{2})
In order to test the relevance of the quasiintegrals 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:
Osipkov (1979) and Merritt (1985) showed that for a spherical potential, any distribution function of the form f(E + L^{2}/2α^{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
The ratio of radialtotangential 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
the radius of the circular orbit r_{circ} of energy E is given by
and we deduced the following quantity related to it
We then set for the distribution function:
with only two free parameters, n and a. We can see that with the potential being logarithmic, the term on the righthand side separates into two terms. A first term contains only the variable r, which gives for the associated density law: ν(r)∝r^{−n}. 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 .
Table 1 gives the velocity dispersions for some values of the a parameter when v_{0} = 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 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(v_{r}, v_{t}) from topleft to bottomright, 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 twobytwo. To do this, we introduced a dependency on the angular momentum L_{z} and we wrote:
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^{−n}, 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 isocontours (logarithmically equally spaced) for distribution function (Eq. (28)) for b and c values of Table 2 (raw and column corresponding). 
σ_{r} σ_{ ϕ} σ_{θ} velocity dispersions at z = 0.
σ_{r} σ_{ϕ} σ_{θ} velocity dispersions at x = y = 0 towards the poles.
4.1.3. f(E, I_{2}, I_{3})
Finally, we considered the triaxial logarithmic potential (see Eq. (21)) with m_{0} = 0.3, q_{y} = 0.95, and q_{z} = 0.85 and we defined the distribution function (with n = 3):
For a triaxially symmetrical potential, the quasiintegrals I_{2} and I_{3} are workable generalisations of and (L^{2} − Lz^{2})/2. The figures presented here are made using for I_{2} and I_{3} 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 nonzero, the potential is no longer a Stäckel potential. The integrals I_{2} and I_{3} are approximated and the distribution function is therefore not strictly stationary. Moreover, the presence of a nonzero core radius, m_{0} = 0.3, means that the potential is no longer strictly scalefree. 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^{−3}.
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 isocontours of the latter model that is dominated with radial orbits.
Fig. 11. Velocity dispersions versus xdistance (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 longdashed lines indicate tangential dispersions σ_{ϕ} and σ_{θ}, respectively. 
Fig. 12. Density for radially biased model (b = c = +0.3). Shown above are equally spaced isocontours 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 quasiintegral orbits are the least wellpreserved 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:
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} = 0.3. This allows the accuracy of the distribution function to be tested and the usefulness of the quasiintegral I_{2} and I_{3} 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 (BlandHawthorn & 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} = 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 F_{x}, F_{y}, F_{z} 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 quasiintegral 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 F_{x}, F_{y}, F_{z} 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, velocitydominant 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, velocitydominant 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 quasiinvariant integral can be achieved by adjusting a Stäckel potential for each of the orbits of a galactic potential. Here, we show that quasiintegrals 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 quasiactions. 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 quasiintegrals 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 quasiintegrals, 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 quasiintegrals are analytic and easier to calculate. As for the actions, they can be evaluated here by a simple quadrature using the expressions of quasiintegrals 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 quasiintegral 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 quasiintegrals 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 isocontours 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; LyndenBell 1962; Tsiganov 2012). These potentials also depend on a generic function and therefore have a certain generality that merits further study.
We remark that the most studied Stäckel potentials have analytic expressions of their three integrals with symmetries in confocal ellipsoidal coordinate systems (Ollongren 1962; LyndenBell 1962; de Zeeuw 1985a). On the other hand, those with symmetries in a parabolic coordinate system have been little studied (Ollongren 1962; LyndenBell 1962; Tsiganov 2012) and left aside.
Acknowledgments
The author would like to thank the International Space Science Institute, Berne, Switzerland, for providing financial support and meeting facilities.
References
 Bienaymé, O., & Traven, G. 2013, A&A, 549, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bienaymé, O., Robin, A. C., & Famaey, B. 2015, A&A, 581, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J. 2012, MNRAS, 426, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, G. 1960, Z. Astrophys., 49, 273 [NASA ADS] [Google Scholar]
 de Zeeuw, T. 1985a, MNRAS, 216, 273 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, T. 1985b, MNRAS, 216, 599 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, P. T., & LyndenBell, D. 1985, MNRAS, 215, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Felberg, E. 1968, NASA Technical Report TR R287 [Google Scholar]
 Helmi, A., Veljanoski, J., Breddels, M. A., Tian, H., & Sales, L. V. 2017, A&A, 598, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hénon, M., & Heiles, C. 1964, AJ, 69, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Hietarinta, J. 1987, Phys. D Nonlinear Phenom., 28, 248 [NASA ADS] [Google Scholar]
 Kent, S. M., & de Zeeuw, T. 1991, AJ, 102, 1994 [NASA ADS] [CrossRef] [Google Scholar]
 Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 1962, MNRAS, 124, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995 [NASA ADS] [CrossRef] [Google Scholar]
 Malhan, K., Ibata, R. A., & Martin, N. F. 2018, MNRAS, 481, 3442 [NASA ADS] [CrossRef] [Google Scholar]
 McGill, C., & Binney, J. 1990, MNRAS, 244, 634 [NASA ADS] [Google Scholar]
 McMillan, P. J., & Binney, J. J. 2008, MNRAS, 390, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D. 1985, AJ, 90, 1913 [Google Scholar]
 Morse, P. M., & Feshbach, H. 1953, Methods of Theoritical Physics, Ch. 5 (McGraw Hill: New York) [Google Scholar]
 Ollongren, A. 1962, Bull. Astron. Inst. Neth., 16, 241 [NASA ADS] [Google Scholar]
 Osipkov, L. P. 1979, Sov. Astron. Lett., 5, 42 [NASA ADS] [Google Scholar]
 Posti, L., & Helmi, A. 2019, A&A, 621, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Posti, L., Binney, J., Nipoti, C., & Ciotti, L. 2015, MNRAS, 447, 3060 [Google Scholar]
 Sanders, J. 2012, MNRAS, 426, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. L., & Binney, J. 2014, MNRAS, 441, 3284 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. L., & Binney, J. 2015, MNRAS, 447, 2479 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. L., & Binney, J. 2016, MNRAS, 457, 2107 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H. 1969, ApJ, 158, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Tsiganov, A. V. 2012, SIGMA, 8, 031 [Google Scholar]
 Vasiliev, E. 2019, MNRAS, 482, 1525 [NASA ADS] [CrossRef] [Google Scholar]
 Whittaker, E. T., & Watson, G. N. 1902, A course in Modern Analysis (Cambridge University Press) [Google Scholar]
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
with α, β, and γ, which are three fixed numbers with −γ ≤ ν ≤ −β ≤ μ ≤ −α ≤ λ.
The expressions for x, y, and z Cartesian coordinates are
The expressions for λ, μ, and ν coordinates are obtained from the relations
and they are the three solutions of the thirddegree polynomial equation
or
with X = Y − b/3. The equation has three solutions given that 4p^{3} + 27q^{2} < 0.
The solutions are , , , with and
For spherical potentials, we have α = β = γ and r^{2} = x^{2} + y^{2} + z^{2} = λ + α.
For prolate potentials, (λ, μ) define oblate spheroidal coordinates and β = γ. The x axis is the axis of axisymmetry and we have , , and .
For oblate potentials, (λ, ν) define prolate spheroidal coordinates and α = β. The z axis is the axis of axisymmetry and we have R^{2} = x^{2} + y^{2}, λ + ν = −α − γ + R^{2} + z^{2}, and λν = αγ − γR^{2} − αz^{2}.
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 I_{2}. The red dashed line indicates dispersion for the integral I_{3}. 

In the text 
Fig. 2. Left: σ_{I2} uncertainties on the I_{2} determinations for 5000 orbits with E = 1 versus their apocentre/pericentre ratio. Right: σ_{I3} uncertainties on the I_{3} 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. Twodimensional classification of orbits. Each point (i_{2}, i_{3}) 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 colourcoded. Right: the σ_{I3} dispersion is colourcoded. 

In the text 
Fig. 6. From top to bottom: shortaxis tube orbit, longaxis tube orbit, and box orbit. 

In the text 
Fig. 7. Poincaré sections. Red: orbit computation. Blue: from the analytic quasiintegrals. Left: (x, v_{x}) section for orbits within the (x, y) plane, (z = ż = 0). Right: (y, v_{y}) section for orbits within the (y, z) plane, x = ẋ = 0. 

In the text 
Fig. 8. Variation of I_{2} (blue) and I_{3} (red) for three orbits. Over the short timelength of the order of an average rotation period (Δt = 3), the variations of the quasiintegrals are significantly smaller compared to the variations over long periods. The largest variations of the quasiintegrals are observed for I_{3}(∼ − 0.2) (bottom curve) of the box orbit (I_{2} ∼ 0). The sharp variations of I_{3} are related to the passage at the pericentre. The shortaxis tube (I_{2} ∼ I_{3} ∼ 0.1) and the longaxis tube (I_{2} ∼ −0.1 and I_{3} ∼ 0.5) orbits have smoother variations of the quasiintegrals. 

In the text 
Fig. 9. Distribution function f(v_{r}, v_{t}) from topleft to bottomright, with a = −2, 1, 0, +1, +2 from tangentially to radially anisotropic velocity distributions. 

In the text 
Fig. 10. Density isocontours (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 xdistance (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 longdashed 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 isocontours 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 F_{x}, F_{y}, F_{z} 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 F_{x}, F_{y}, F_{z} 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.