Integrals of motion for non-axisymmetric potentials

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\"ackel 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 box-type orbits.

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.
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.

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 double-precision 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 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: 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 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 = 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 ∆ 1 = √ β − α, ∆ 2 = √ γ − β and ∆ 3 = √ γ − α (positions of the foci of the ellipsoidal coordinate system) to obtain the best conservation of the quasi-integral 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 quasi-integrals 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. 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 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).
It is interesting to compare these integrals with the angular momenta |L z | and L ⊥ = L 2 x + L 2 y 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 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. Figure 4 shows the correlation between the accuracies σ I 3 and σ L ⊥ . For σ I2 and σ L z 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 three-dimensional classification of orbits in terms of E, I 2 , and I 3 for bound orbits   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.
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 quasi-integrals are significantly more stationary.
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 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.

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,   (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. 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 2 and I 3 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 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.

Prescribed halo distribution functions
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: 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 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 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 right-hand 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 σ t / √ 2 refers to the dispersion of a single tangential component. Figure 9 represents the corresponding f (v r , v t ) distributions.

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 , , , , , Fig. 9. Distribution function f (v r , v t ) from top-left to bottom-right, with a = −2, 1, 0, +1, +2 from tangentially to radially anisotropic velocity distributions. Table 2. b and c co-efficients for models with density plotted in Fig. 10.  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).

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 quasi-integrals I 2 and I 3 are workable generalisations of L 2 z /2 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 non-zero, 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 non-zero core radius, m 0 = 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 −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 iso-contours of the latter model that is dominated with radial orbits.
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).

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 quasi-integral 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 (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 = 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. 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.
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.

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 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.