Integrals of motion for non-axisymmetric potentials

Olivier Bienayme

arXiv:1906.02078·astro-ph.GA·July 17, 2019

Integrals of motion for non-axisymmetric potentials

Olivier Bienayme

PDF

TL;DR

This paper develops explicit integrals of motion for non-axisymmetric galactic potentials, enabling better modeling of stellar populations and gravitational fields with a focus on triaxial and Stäckel potentials.

Contribution

It introduces an analytic formulation of integrals of motion with explicit potential dependence applicable to non-axisymmetric potentials, including approximate solutions for general cases.

Findings

01

Accurately models stationary stellar populations using derived integrals.

02

Recovers force fields with satisfactory accuracy.

03

Less precise in reconstructing mass density distributions.

Abstract

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.

Figures40

Click any figure to enlarge with its caption.

Tables4

Table 1. Table 1: Velocity moments for Eq. 27 distribution function with v 0 = 1 subscript 𝑣 0 1 v_{0}=1 and n = + 3 𝑛 3 n=+3 .

$a$	$σ_{r}$	$\frac{σ_{t}}{\sqrt{2}}$	$\frac{σ_{r}}{σ_{t} \sqrt{2}}$	$< \| v_{t} \| >$	$σ (\| v_{t} \|)$
-10	0.13	0.70	0.19	0.99	0.09
-2	0.30	0.67	0.45	0.93	0.21
-1	0.41	0.64	0.64	0.87	0.29
0	0.577	0577	1	0.724	0.378
+1	0.72	0.49	1.49	0.55	0.41
+2	0.80	0.42	1.91	0.44	0.40
+10	0.92	0.28	3.35	0.22	0.32

Table 2. Table 2: b 𝑏 b and c 𝑐 c co-efficients for models with density plotted in Figure 10 . Each row and column of the table correspond to those of Figure 10 , and of Tables 3 and 4 .

$+ 4$ , $+ 0$	+2, $- 2$	$+ 0$ , $- 4$
$+ 2$ , $+ 2$	+0, +0	$- 2$ , $- 2$
$+ 0$ , $+ 4$	$- 2$ , $+ 2$	$- 4$ , $+ 0$

Table 3. Table 3: σ r subscript 𝜎 𝑟 \sigma_{r} σ ϕ subscript 𝜎 italic-ϕ \sigma_{\phi} σ θ subscript 𝜎 𝜃 \sigma_{\theta} velocity dispersions at z = 0 𝑧 0 z=0 . Each row and column of the table correspond to those of Figure 10 , and of Tables 2 and 4

$0.68 0.28 0.28$	$0.33 0.23 0.92$	$0.22 0.22 0.95$
$0.80 0.42 0.42$	$0.58 0.58 0.58$	$0.30 0.67 0.67$
$0.68 0.68 0.28$	$0.33 0.92 0.23$	$0.22 0.95 0.22$

Table 4. Table 4: σ r subscript 𝜎 𝑟 \sigma_{r} σ ϕ subscript 𝜎 italic-ϕ \sigma_{\phi} σ θ subscript 𝜎 𝜃 \sigma_{\theta} velocity dispersions at x = y = 0 𝑥 𝑦 0 x=y=0 towards the poles. Each row and column of the table correspond to those of Figure 10 , and of Tables 2 and 3

$0.87 0.35 0.35$	$0.80 0.42 0.42$	$0.58 0.58 0.58$
$0.80 0.42 0.42$	$0.58 0.58 0.58$	$0.30 0.67 0.67$
$0.58 0.58 0.58$	$0.30 0.67 0.67$	$0.21 0.69 0.69$

Equations76

Φ (λ, μ, ν) = - 4 \frac{ζ ( λ ) ( λ + α ) ( λ + β ) ( λ + γ )}{( λ - μ ) ( λ - ν )} - 4 \frac{η ( μ ) ( μ + α ) ( μ + β ) ( μ + γ )}{( μ - ν ) ( μ - λ )} - 4 \frac{θ ( ν ) ( ν + α ) ( ν + β ) ( ν + γ )}{( ν - λ ) ( ν - μ )} .

Φ (λ, μ, ν) = - 4 \frac{ζ ( λ ) ( λ + α ) ( λ + β ) ( λ + γ )}{( λ - μ ) ( λ - ν )} - 4 \frac{η ( μ ) ( μ + α ) ( μ + β ) ( μ + γ )}{( μ - ν ) ( μ - λ )} - 4 \frac{θ ( ν ) ( ν + α ) ( ν + β ) ( ν + γ )}{( ν - λ ) ( ν - μ )} .

ζ (λ) = - 1/4 \frac{Φ ( λ , - β , - γ )}{λ + α} .

ζ (λ) = - 1/4 \frac{Φ ( λ , - β , - γ )}{λ + α} .

η (μ) = - 1/4 \frac{Φ ( - α , μ , - γ )}{μ + β},

η (μ) = - 1/4 \frac{Φ ( - α , μ , - γ )}{μ + β},

θ (ν) = - 1/4 \frac{Φ ( - α , - β , ν )}{ν + γ} .

θ (ν) = - 1/4 \frac{Φ ( - α , - β , ν )}{ν + γ} .

Φ (λ, μ, ν) = \frac{( λ + β ) ( λ + γ )}{( λ - μ ) ( λ - ν )} Φ (λ, - β, - γ) + \frac{( μ + α ) ( μ + γ )}{( μ - ν ) ( μ - λ )} Φ (- α, μ, - γ) + \frac{( ν + α ) ( ν + β )}{( ν - λ ) ( ν - μ )} Φ (- α, - β, ν),

Φ (λ, μ, ν) = \frac{( λ + β ) ( λ + γ )}{( λ - μ ) ( λ - ν )} Φ (λ, - β, - γ) + \frac{( μ + α ) ( μ + γ )}{( μ - ν ) ( μ - λ )} Φ (- α, μ, - γ) + \frac{( ν + α ) ( ν + β )}{( ν - λ ) ( ν - μ )} Φ (- α, - β, ν),

H = Φ (λ, μ, ν) + 1/2 v_{x}^{2} + 1/2 v_{y}^{2} + 1/2 v_{z}^{2}

H = Φ (λ, μ, ν) + 1/2 v_{x}^{2} + 1/2 v_{y}^{2} + 1/2 v_{z}^{2}

J = Ψ (λ, μ, ν) + f_{J} (L^{2}, v_{x}, v_{y}, v_{z})

J = Ψ (λ, μ, ν) + f_{J} (L^{2}, v_{x}, v_{y}, v_{z})

K = Ξ (λ, μ, ν) + f_{K} (L_{x}^{2}, L_{y}^{2}, L_{z}^{2}, v_{x}, v_{y}, v_{z})

K = Ξ (λ, μ, ν) + f_{K} (L_{x}^{2}, L_{y}^{2}, L_{z}^{2}, v_{x}, v_{y}, v_{z})

I_{2} = \frac{α ^{2} H + α J + L}{γ - α} = ψ (λ, μ, ν) + 1/2 (α - β) v_{x}^{2} + 1/2 \frac{( α - β )}{( α - γ )} L_{y}^{2} + 1/2 L_{z}^{2}

I_{2} = \frac{α ^{2} H + α J + L}{γ - α} = ψ (λ, μ, ν) + 1/2 (α - β) v_{x}^{2} + 1/2 \frac{( α - β )}{( α - γ )} L_{y}^{2} + 1/2 L_{z}^{2}

I_{3} = \frac{γ ^{2} H + γ J + L}{α - γ} = ξ (λ, μ, ν) + 1/2 (γ - β) v_{z}^{2} + 1/2 L_{x}^{2} + 1/2 \frac{( γ - β )}{( γ - α )} L_{y}^{2}

I_{3} = \frac{γ ^{2} H + γ J + L}{α - γ} = ξ (λ, μ, ν) + 1/2 (γ - β) v_{z}^{2} + 1/2 L_{x}^{2} + 1/2 \frac{( γ - β )}{( γ - α )} L_{y}^{2}

ψ = \frac{( μ + α ) ( ν + α )}{α - γ} [\frac{( λ + β ) ( λ + γ )}{( λ - μ ) ( λ - ν )} Φ (λ, - β, - γ) + \frac{( μ + γ ) ( λ + α )}{( μ - ν ) ( μ - λ )} Φ (- α, μ, - γ) + \frac{( ν + β ) ( λ + α )}{( ν - λ ) ( ν - μ )} Φ (- α, - β, ν)],

ψ = \frac{( μ + α ) ( ν + α )}{α - γ} [\frac{( λ + β ) ( λ + γ )}{( λ - μ ) ( λ - ν )} Φ (λ, - β, - γ) + \frac{( μ + γ ) ( λ + α )}{( μ - ν ) ( μ - λ )} Φ (- α, μ, - γ) + \frac{( ν + β ) ( λ + α )}{( ν - λ ) ( ν - μ )} Φ (- α, - β, ν)],

ξ = \frac{( μ + γ ) ( λ + γ )}{γ - α} [\frac{( λ + β ) ( ν + γ )}{( λ - μ ) ( λ - ν )} Φ (λ, - β, - γ) + \frac{( μ + α ) ( ν + γ )}{( μ - ν ) ( μ - λ )} Φ (- α, μ, - γ) + \frac{( ν + α ) ( ν + β )}{( ν - λ ) ( ν - μ )} Φ (- α, - β, ν)] .

ξ = \frac{( μ + γ ) ( λ + γ )}{γ - α} [\frac{( λ + β ) ( ν + γ )}{( λ - μ ) ( λ - ν )} Φ (λ, - β, - γ) + \frac{( μ + α ) ( ν + γ )}{( μ - ν ) ( μ - λ )} Φ (- α, μ, - γ) + \frac{( ν + α ) ( ν + β )}{( ν - λ ) ( ν - μ )} Φ (- α, - β, ν)] .

Φ (λ, μ, ν) = \frac{( μ + γ )}{( μ - ν )} Φ (λ, μ, - γ) + \frac{( ν + β )}{( ν - μ )} Φ (λ, - β, ν) + \frac{( β - γ )}{( μ - ν )} Φ (λ, - β, - γ),

Φ (λ, μ, ν) = \frac{( μ + γ )}{( μ - ν )} Φ (λ, μ, - γ) + \frac{( ν + β )}{( ν - μ )} Φ (λ, - β, ν) + \frac{( β - γ )}{( μ - ν )} Φ (λ, - β, - γ),

ψ (λ, μ, ν) = \frac{( λ + α ) ( μ + γ ) ( ν + α )}{( μ - ν ) ( α - γ )} Φ (λ, μ, - γ) - \frac{( λ + α ) ( μ + α ) ( ν + β )}{( μ - ν ) ( α - γ )} Φ (λ, - β, ν) - [\frac{( λ + α ) ( λ + β ) ( μ + γ ) ( ν + α )}{( λ - μ ) ( μ - ν ) ( α - γ )} + \frac{( λ + α ) ( λ + γ ) ( μ + α ) ( ν + β )}{( ν - λ ) ( μ - ν ) ( α - γ )} + \frac{( λ + β ) ( λ + γ ) ( μ + α ) ( ν + α )}{( λ - μ ) ( ν - λ ) ( α - γ )}] Φ (λ, - β, - γ),

ψ (λ, μ, ν) = \frac{( λ + α ) ( μ + γ ) ( ν + α )}{( μ - ν ) ( α - γ )} Φ (λ, μ, - γ) - \frac{( λ + α ) ( μ + α ) ( ν + β )}{( μ - ν ) ( α - γ )} Φ (λ, - β, ν) - [\frac{( λ + α ) ( λ + β ) ( μ + γ ) ( ν + α )}{( λ - μ ) ( μ - ν ) ( α - γ )} + \frac{( λ + α ) ( λ + γ ) ( μ + α ) ( ν + β )}{( ν - λ ) ( μ - ν ) ( α - γ )} + \frac{( λ + β ) ( λ + γ ) ( μ + α ) ( ν + α )}{( λ - μ ) ( ν - λ ) ( α - γ )}] Φ (λ, - β, - γ),

ξ (λ, μ, ν) = \frac{( λ + γ ) ( μ + γ )}{( μ - ν ) ( γ - α )} [(ν + γ) Φ (λ, μ, - γ) - (ν + β) Φ (λ, - β, ν) - (γ - β) Φ (λ, - β, - γ)] .

ξ (λ, μ, ν) = \frac{( λ + γ ) ( μ + γ )}{( μ - ν ) ( γ - α )} [(ν + γ) Φ (λ, μ, - γ) - (ν + β) Φ (λ, - β, ν) - (γ - β) Φ (λ, - β, - γ)] .

ψ (λ, μ, ν) = \frac{( λ + α ) ( μ + α )}{α - γ} Φ (λ, μ, ν) - \frac{( λ + α ) ( μ + γ )}{α - γ} Φ (λ, μ, - γ) - (λ + β) Φ (λ, - β, - γ),

ψ (λ, μ, ν) = \frac{( λ + α ) ( μ + α )}{α - γ} Φ (λ, μ, ν) - \frac{( λ + α ) ( μ + γ )}{α - γ} Φ (λ, μ, - γ) - (λ + β) Φ (λ, - β, - γ),

ξ (λ, μ, ν) = \frac{( λ + γ ) ( μ + γ )}{γ - α} [Φ (λ, μ, ν) - Φ (λ, μ, - γ)] .

ξ (λ, μ, ν) = \frac{( λ + γ ) ( μ + γ )}{γ - α} [Φ (λ, μ, ν) - Φ (λ, μ, - γ)] .

ξ (λ, ν) = (λ + γ) [Φ (λ, ν) - Φ (λ, - γ)] .

ξ (λ, ν) = (λ + γ) [Φ (λ, ν) - Φ (λ, - γ)] .

ξ (λ, ν) = \frac{( λ + γ )}{( λ - ν )} [(ν + γ) Φ (λ, - γ) - (ν + α) Φ (- α, ν)]

ξ (λ, ν) = \frac{( λ + γ )}{( λ - ν )} [(ν + γ) Φ (λ, - γ) - (ν + α) Φ (- α, ν)]

Φ = (β + γ) x^{2} + (γ + α) y^{2} + (α + β) z^{2} + (- α - β - γ + x^{2} + y^{2} + z^{2})^{2}

Φ = (β + γ) x^{2} + (γ + α) y^{2} + (α + β) z^{2} + (- α - β - γ + x^{2} + y^{2} + z^{2})^{2}

Φ (x, y, z) = \frac{1}{2} lo g (m_{0}^{2} + m^{2})

Φ (x, y, z) = \frac{1}{2} lo g (m_{0}^{2} + m^{2})

m^{2} = x^{2} + \frac{y ^{2}}{q _{y}^{2}} + \frac{z ^{2}}{q _{z}^{2}}

m^{2} = x^{2} + \frac{y ^{2}}{q _{y}^{2}} + \frac{z ^{2}}{q _{z}^{2}}

f (E, L) = exp [- (E + L^{2} /2 α^{2})] .

f (E, L) = exp [- (E + L^{2} /2 α^{2})] .

\frac{σ _{r}^{2}}{σ _{t}^{2}} = \frac{σ _{r}^{2}}{\frac{1}{2} < v _{t}^{2} >} = 1 + \frac{r ^{2}}{α ^{2}} .

\frac{σ _{r}^{2}}{σ _{t}^{2}} = \frac{σ _{r}^{2}}{\frac{1}{2} < v _{t}^{2} >} = 1 + \frac{r ^{2}}{α ^{2}} .

Φ (r) = v_{0}^{2} lo g r,

Φ (r) = v_{0}^{2} lo g r,

E = v_{0}^{2} lo g (r_{c i r c}) + 1/2 v_{0}^{2},

E = v_{0}^{2} lo g (r_{c i r c}) + 1/2 v_{0}^{2},

r_{c} (E) = exp [E / v_{0}^{2} - 1/2] = r exp [\frac{1}{2} \frac{v _{r}^{2} + v _{t}^{2}}{v _{0}^{2}} - \frac{1}{2}] .

r_{c} (E) = exp [E / v_{0}^{2} - 1/2] = r exp [\frac{1}{2} \frac{v _{r}^{2} + v _{t}^{2}}{v _{0}^{2}} - \frac{1}{2}] .

f (E, L) = exp [- \frac{n}{v _{0}^{2}} (E + a \frac{L ^{2}}{r _{c} ( E ) ^{2}})],

f (E, L) = exp [- \frac{n}{v _{0}^{2}} (E + a \frac{L ^{2}}{r _{c} ( E ) ^{2}})],

f = exp {- n [E + \frac{1}{r _{E}^{2}} (b L_{z}^{2} + c (L^{2} - L_{z}^{2}))]} .

f = exp {- n [E + \frac{1}{r _{E}^{2}} (b L_{z}^{2} + c (L^{2} - L_{z}^{2}))]} .

f = exp [- n (E + \frac{1}{r _{E}^{2}} (b I_{2} + c I_{3}))] .

f = exp [- n (E + \frac{1}{r _{E}^{2}} (b I_{2} + c I_{3}))] .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

11institutetext: Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, 11 rue de l’Université, F-67000 Strasbourg, France

Integrals of motion for non-axisymmetric potentials

O. Bienaymé 11

Abstract

Aims.

Methods.

Results.

Abstract

*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 box-type orbits.

Key Words.:

Methods: numerical – Galaxies: kinematics and dynamics

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 symmetry111 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; Lynden-Bell, 1962; de Zeeuw, 1985a). On the other hand, those with symmetries in a parabolic coordinate system have been little studied (Ollongren, 1962; Lynden-Bell, 1962; Tsiganov, 2012) and left aside.. 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 Section 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 Section 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 Section 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 $\zeta(\lambda)$ , $\eta(\mu),$ and $\theta(\nu)$ :

[TABLE]

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 $\Phi$ consists in equating both potentials along these three different axes. Other possibilities are also presented.

If $\Phi$ is a continuous function at the origin, the formulation in Eq.1 implies that $\Phi(-\alpha,-\beta,-\gamma)=0$ . Evaluating $\Phi$ at positions $(\lambda,\mu=-\beta,\nu=-\gamma)$ , we deduce $\zeta(\lambda)$ :

[TABLE]

Similarly, $\eta(\mu)$ and $\theta(\nu)$ are obtained evaluating $\Phi\left(-\alpha,-\beta,\nu\right)$ and $\Phi\left(-\alpha,\mu,-\gamma\right)$ :

[TABLE]

Then, substituting these expressions of $\zeta,\eta,$ and $\theta$ within the expression of $\Phi$ (Eq.1), we obtain

[TABLE]

a valid expression if $\Phi(-\alpha,-\beta,-\gamma)=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 equations 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 & Lynden-Bell (1985, equations 7 and 14), and for convenience write the exact integrals in compact form as

[TABLE]

and

[TABLE]

According to de Zeeuw & Lynden-Bell (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 ( $\alpha=\beta$ ) with the $z$ axis for the rotational symmetry, $I_{2}=1/2\,L_{z}^{2}$ . 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 $\zeta,\eta,$ and $\theta$ , 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 $\Phi(-\alpha,-\beta,-\gamma)=0$ ):

[TABLE]

Other choices of the three positions are possible but require more algebra. However, a simple case can be deduced, again evaluating $\Phi$ at three positions; determining $\zeta$ , $\eta$ , and $\theta$ ; and substituting in $\Phi$ , $\psi,$ and $\xi$ 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):

[TABLE]

A simpler combination is obtained after substitution of $\Phi(\lambda,-\beta,\nu)$ deduced from Eq. 13:

[TABLE]

We note that within the plane $z=0$ , $\xi\left(\lambda,\mu,\nu=-\gamma\right)$ is null. Moreover, along the $x$ axis, $\psi\left(\lambda,-\beta,-\gamma\right)=\left(-\beta+\alpha\right)\Phi\left(\lambda,-\beta,-\gamma\right),$ is null in the case of an axisymmetric potential with $\alpha=\beta$ .

Other expressions of the three integrals depending on $\Phi$ 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 $(\lambda,\mu,\nu)$ position.

In the case of an axisymmetric potential with $\alpha=\beta$ , 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, equation A7).

[TABLE]

For an axisymmetric potential with setting $\alpha=\beta$ , 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

[TABLE]

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

[TABLE]

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 $\alpha=-2$ , $\beta=-1$ and $\gamma=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) and Sanders & Binney (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:

[TABLE]

where

[TABLE]

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_{\Phi}$ correspond to density flattening values of $q_{\rho}\sim 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_{\rho}=1.22\pm$ 0.23 from the analysis of globular cluster orbits while Malhan & Ibata (2018) found a flattening of the density $q_{\rho}=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 (Fehlberg, 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 $\Delta_{1}=\sqrt{\beta-\alpha}$ , $\Delta_{2}=\sqrt{\gamma-\beta}$ and $\Delta_{3}=\sqrt{\gamma-\alpha}$ (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 $\Delta_{1}=0.49$ and $\Delta_{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, $\sigma_{I2}$ or $\sigma_{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\times 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 $\sigma_{I2}$ and $\sigma_{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 $\sigma_{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_{\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_{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 $\sigma_{I_{3}}$ and $\sigma_{L_{\perp}}$ . For $\sigma_{I2}$ and $\sigma_{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 of a triaxial Stäckel potential and a cross-section for a given energy $E$ . Figures 5b and 5c show the distribution of orbits with their accuracy colour-coded.

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 $\Delta_{1}$ and $\Delta_{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 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.

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.

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, 2018). 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.

4.1 Prescribed halo distribution functions

4.1.1 $f(E,L^{2}$ )

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:

[TABLE]

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

[TABLE]

The ratio of radial-to-tangential velocity dispersions changes continuously with distance to the centre, and the parameter $\alpha$ defines the transition radius. In the specific case of a logarithmic potential, this scale factor $\alpha$ 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

[TABLE]

the radius of the circular orbit $r_{circ}$ of energy E is given by

[TABLE]

and we deduced the following quantity related to it

[TABLE]

We then set for the distribution function:

[TABLE]

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: $\nu(r)\propto 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 $\sigma_{r}=\sigma_{t}/\sqrt{2}=v_{0}/\sqrt{n}$ .

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=-\infty$ or $a=+\infty$ correspond to distributions with purely tangential or purely radial orbits, respectively. The tangential velocity has two components, according to $\theta$ and $\phi$ . $|v_{t}|$ refers to the tangential velocity modulus, while $\sigma_{t}/\sqrt{2}$ refers to the dispersion of a single tangential component. Figure 9 represents the corresponding $f(v_{r},v_{t})$ distributions.

4.1.2 $f(E,L_{z},L_{\perp}$ )

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 $\sigma_{r}$ , $\sigma_{\phi},$ and $\sigma_{\theta},$ different two-by-two. To do this, we introduced a dependency on the angular momentum $L_{z}$ and we wrote:

[TABLE]

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 $\theta$ (the angle between the vertical axis $z$ and the vector radius $\@vec{r}$ , and $\phi$ is the azimuth angle). The velocity distribution does not depend on the distance to the centre $r$ but only on the angle $\theta$ .

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_{\phi},v_{\theta})$ . For $b=c,$ the distribution function has a spherical density (see Eq. 28). In these cases, the velocity dispersions $\sigma_{\phi}$ and $\sigma_{\theta}$ are equal and do not depend on $r$ , $\phi,$ or $\theta$ (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 Figure 10 (also the first and third columns). It should be noted that if $b<c$ , in the $z=0$ plane the vertical velocity dispersion $\sigma_{\theta}$ 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 $\sigma_{\phi}/\sigma_{r}$ and $\sigma_{\theta}/\sigma_{r}$ ratios arbitrarily chosen within the $z=0$ plane. Outside the $z=0$ plane towards the poles $\theta=0$ or $\pi/2$ , the velocity distributions along the three main axes vary, and, for symmetry reasons, the dispersions $\sigma_{\theta}$ and $\sigma_{\phi}$ become equal to the poles (Table 4).

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$ ):

[TABLE]

For a triaxially symmetrical potential, the quasi-integrals $I_{2}$ and $I_{3}$ are workable generalisations of $L_{z}^{2}/2$ and $(L^{2}-Lz^{2})/2$ . The figures presented here are made using for $I_{2}$ and $I_{3}$ the expressions given by the 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 $\sim 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 $\sigma_{\phi}\simeq\sigma_{\theta}$ . 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).

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:

[TABLE]

These equations allow us to calculate the force field associated with the gravitational potential $\Phi$ . 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 $\Phi$ (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 $\nabla\Phi,$ 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 $\sim 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, $\rho_{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 $\rho_{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 $\rho_{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.

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

Acknowledgements.

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

Appendix A $(x,y,z)$ $(\lambda,\mu,\nu)$ transformation

Stäckel potentials are defined using ellipsoidal coordinates. The $(\lambda,\mu,\nu)$ coordinates define confocal ellipsoids and hyperboloids and are roots of the equation

[TABLE]

with $\alpha$ , $\beta$ , and $\gamma,$ which are three fixed numbers with $-\gamma\leq\nu\leq-\beta\leq\mu\leq-\alpha\leq\lambda$ .

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

[TABLE]

The expressions for $\lambda$ , $\mu,$ and $\nu$ coordinates are obtained from the relations

[TABLE]

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

[TABLE]

or

[TABLE]

with $X=Y-b/3.$ The equation has three solutions given that $4p^{3}+27q^{2}<0$ .

The solutions are $Y_{1}=r\sin\frac{\varphi}{3}$ , $Y_{2}=r\sin\left(\frac{\varphi}{3}+\frac{2\pi}{3}\right)$ , $Y_{3}=r\sin\left(\frac{\varphi}{3}+\frac{4\pi}{3}\right)$ , with $r=2\sqrt{-\frac{p}{3}}$ and $\sin\varphi=-\frac{3q}{pr}$ .

For spherical potentials, we have

$\alpha=\beta=\gamma$ and $r^{2}=x^{2}+y^{2}+z^{2}=\lambda+\alpha$ .

For prolate potentials, $(\lambda,\mu)$ define oblate spheroidal coordinates and

$\beta=\gamma$ . The $x$ axis is the axis of axisymmetry and we have $\tilde{z}^{2}=y^{2}+z^{2}$ ,

$\lambda+\mu=-\alpha-\gamma+x^{2}+\tilde{z}^{2}$ , and $\lambda\mu=\alpha\gamma-\gamma x^{2}-\alpha\tilde{z}^{2}$ .

For oblate potentials, $(\lambda,\nu)$ define prolate spheroidal coordinates and

$\alpha=\beta$ . The $z$ axis is the axis of axisymmetry and we have $R^{2}=x^{2}+y^{2}$

$\lambda+\nu=-\alpha-\gamma+R^{2}+z^{2}$ , and $\lambda\nu=\alpha\gamma-\gamma R^{2}-\alpha z^{2}$ .

Bibliography34

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Bienaymé et al. (2015) Bienaymé, O., Robin, A. C., & Famaey, B. 2015, A&A, 581, A 123
2Bienaymé & Traven (2013) Bienaymé, O., & Traven, G. 2013, A&A, 549, A 89
3Binney (2012) Binney, J. 2012, MNRAS, 426, 1324
4Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529
5Bovy (2015) Bovy, J. 2015, Ap JS, 216, 29
6Contopoulos (1960) Contopoulos, G. 1960, Z Ap, 49, 273
7de Zeeuw (1985 a) de Zeeuw, T. 1985 a, MNRAS, 216, 273
8de Zeeuw (1985 b) de Zeeuw, T. 1985 b, MNRAS, 216, 599

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Integrals of motion for non-axisymmetric potentials

Abstract

Abstract

Key Words.:

1 Introduction

2 Triaxial Stäckel potentials and their integrals of motion

3 Tests

4 Distribution functions and Jeans equations

4.1 Prescribed halo distribution functions

4.1.1 f(E,L2f(E,L^{2}f(E,L2)

4.1.2 f(E,Lz,L⊥f(E,L_{z},L_{\perp}f(E,Lz​,L⊥​)

4.1.3 f(E,I2,I3f(E,I_{2},I_{3}f(E,I2​,I3​)

4.2 Test with Jeans equations

5 Conclusion

Acknowledgements.

Appendix A (x,y,z)(x,y,z)(x,y,z) (λ,μ,ν)(\lambda,\mu,\nu)(λ,μ,ν) transformation

4.1.1 $f(E,L^{2}$ )

4.1.2 $f(E,L_{z},L_{\perp}$ )

4.1.3 $f(E,I_{2},I_{3}$ )

Appendix A $(x,y,z)$ $(\lambda,\mu,\nu)$ transformation