Issue 
A&A
Volume 581, September 2015



Article Number  A123  
Number of page(s)  9  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201526516  
Published online  21 September 2015 
Quasi integral of motion for axisymmetric potentials
^{1}
Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS,
UMR 7550, 11 rue de l’Université,
67000
Strasbourg,
France
email:
olivier.bienayme@unistra.fr
^{2}
Institut Utinam, CNRS UMR 6213, Université de FrancheComté, OSU
THETA FrancheComtéBourgogne, Observatoire de Besançon, BP 1615, 25010
Besançon Cedex,
France
Received: 12 May 2015
Accepted: 20 July 2015
We present an estimate of the third integral of motion for axisymmetric threedimensional potentials. This estimate is based on a Stäckel approximation and is explicitly written as a function of the potential. We tested this scheme for the Besançon Galactic model and two other dischalo models and find that orbits of disc stars have an accurately conserved third quasi integral. The accuracy ranges from of 0.1% to 1% for heights varying from z = 0 kpc to z = 6 kpc and Galactocentric radii R from 5 to 15 kpc. We also tested the usefulness of this quasi integral in analytic distribution functions of disc stellar populations: we show that the distribution function remains approximately stationary and that it allows to recover the potential and forces by applying Jeans equations to its moments.
Key words: methods: numerical / Galaxy: kinematics and dynamics
© ESO, 2015
1. Introduction
Many approaches exist for the galactic modelling of stellar dynamics. With the availability of a large amount of new accurate data for the velocity and position of stars in the extended solar neighbourhood, the development of the most precise tools for analysing the kinematics of stars is becoming inescapable.
Numerical modelling techniques are the most direct ones, such as for the Scharwschild method, the madetomeasure method (Syer & Tremaine 1996), or of course ab initio Nbody and hydrodynamical simulations (Renaud et al. 2013). Spectral analysis (Papaphilippou & Laskar 1998; Valluri & Merritt 1998) also offers useful numerical tools for recognising and identifying resonances. Numerical resolution of the collisionless Boltzmann equation is another direct way to model, but it is still limited by available numerical resources (Yoshikawa et al. 2013; Colombi et al. 2015).
Analytical techniques are more explicit but imply analytical approximations. In the case of axisymmetric galactic potentials, it is known that orbits are generally constrained by an effective third integral of motion in addition to the energy and angular momentum. From the Jeans theorem, therefore, equilibrium models of axisymmetric galaxies should be represented by distribution functions that only depend on these three integrals. The modelling of stellar distribution functions can thus be achieved if the effective third integral can be approximated analytically or numerically.
A long list of works in this direction already exists, and very useful bibliographies may be found in de Zeeuw (1985), de Zeeuw & LyndenBell (1985), De Bruyne et al. (2000), Binney (2012), and Sanders & Binney (2014), among others. We concentrate here on the use of Stäckel potentials that have been shown in many cases to be efficient for modelling different families of orbits. Published works on this subject present a wide variety of approaches, and we can distinguish between the local and global approaches. The local modelling consists in locally fitting the true potential with a Stäckel potential. This allows nearly exact modelling of orbits and integrals of motion to be obtained in the immediate neighbourhood of the position considered.
Global modelling, either of the potential or of the orbits over a predefined phasespace volume, has opened the way to many different technicalities over the past 50 yr. In each case the goal was to use a Stäckel third integral as an approximate of the third quasi integral of the fitted potential, as in Wayman (1959), van de Hulst (1962), Ollongren (1962), Hori (1962), and more recently, Manabe (1979), Batsleer & Dejonghe (1994), or Famaey & Dejonghe (2003).
A recent novel application of the Stäckel potential approximations has been to model the local stellar kinematics and the gravitational potential in the solar neighbourhood at large z distances up to 1−2 kpc from the galactic plane (Bienaymé et al. 2014; Piffl et al. 2014). These recent studies put strict constraints on the vertical variation in the gravitational potential and in the local density of the dark matter halo. Recently, Sanders & Binney (2014) have proposed a global Stäckel fitting of 3D potentials, and they give algebraic expressions that allow recovery of the integrals of motion, expressed as action variables, and they present numerical applications.
In the present paper, we proceed to a Stäckel potential fitting by using a simple expression for the integral of motion that explicitly depends on the potential. Our study has similarities with the works recently published by Binney (2012) or by Sanders & Binney (2014), proposing different formulations of a third integral. Although there are also advantages in working with action integrals, the present approach is, however, much simpler and more straightforward to apply.
The paper proceeds as follows. In Sect. 2, we give a new expression for the third integral, and its derivation is described in the Appendix. In Sect. 3 we examine, for three potentials, the constancy of this third integral along orbits. Section 4 shows, in the case of the Besançon Galactic model, how moments of distribution functions based on this third integral are in accordance with Jeans equations. Section 5 finally considers the application to the collisionless Boltzmann equation.
2. Quasi integral of motion
We let V(R,z) be an axisymmetric gravitational potential. Besides the energy E and the vertical component of the angular momentum L_{z}, which are integrals of motion, we define here, as an approximate third integral of motion, ${\mathit{I}}_{\mathrm{s}}\mathrm{=}\mathrm{\Psi}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{L}}^{\mathrm{2}}\mathrm{}{\mathit{L}}_{\mathit{z}}^{\mathrm{2}}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}$(1)where L and L_{z} are the total and vertical angular momenta, z_{0} a fixed parameter at fixed E and L_{z}, v_{z} the vertical velocity, and with $\mathrm{\Psi}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}\mathrm{}\left[\mathit{V}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{}\mathit{V}\mathrm{\left(}\sqrt{\mathit{\lambda}}\mathit{,}\mathrm{0}\mathrm{\right)}\right]\hspace{0.17em}\frac{\mathrm{(}\mathit{\lambda}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}$(2)and $\begin{array}{ccc}\mathit{\lambda}\mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathrm{}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\hspace{0.17em}\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathrm{}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\hspace{0.17em}{\mathit{R}}^{\mathrm{2}}\hspace{0.17em}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{.}& \end{array}$(3)This approximate integral I_{s} is an exact integral in the case of Stäckel potentials expressed in an ellipsoidal coordinate system of focus z_{0}, which is thus a known quantity. Its derivation is explained in the Appendix.
The only free quantity for determining I_{s}, when V(R,z) is known, is the parameter z_{0}. Here, it will be numerically adjusted by minimizing σ_{Is} the dispersion of the quasi integral I_{s} along all the orbits with the same energy and vertical angular momentum. Thus z_{0} will itself be a function of the two integrals E and L_{z}.
In a preliminary analysis, which was not developed further, we used the generic differential equation of Stäckel potentials that can be rewritten to determine z_{0} locally as a function of the second derivatives of the potential at any positions (see Eq. (10) in Bienaymé 2009; de Zeeuw 1985). However, we find out that such a procedure gives us significantly fewer accurate results than fitting a single z_{0} for a given orbit or family of orbits. Here, we look for the z_{0} value that minimizes the dispersion of I_{s} along each orbit of a family of orbits with the same E and L_{z}. We note, however, that for any (E, L_{z})family of orbits and with z_{shell} being the maximum zextent of the shell orbit, the value of z_{0} that gives the best fit is close to the value of z_{0} obtained by applying Eq. (10) of Bienaymé (2009) at position (R_{c},z_{shell}/ 2) (with R_{c}(L_{z}) the radius of the circular orbit with angular momentum L_{z}).
The degree of the approximation of the quasi integral I_{s} can be controlled in several ways: (1) inspection of the surfaces of section; (2) conservation of orbital weights and of the spatial density; and (3) conservation of I_{s} along the orbits to validate the labelling of orbits. The conservation of orbital weights and the conservation of I_{s} are the criteria that can be best expressed numerically.
In the following sections, we test the degree of conservation of this quasi integral along orbits in the cases of three nonStäckel potentials. We also determine its efficiency to build stationary stellar disc distribution functions and to measure the potential from the Jeans equations.
3. Testing the conservation of the integral along orbits
We examined the conservation of the quasiintegral of motion along orbits of stars corresponding to disc stellar populations. For these disc components, stellar motions have restricted oscillations in the radial and vertical directions. For the smallest oscillations, the phasespace domain explored by stars is the closest to a Stäckel potential, and the quasi integral tends to be constant.
We considered stars with identical angular momentum L_{z} (with E_{c} the energy of the corresponding circular orbit) and identical energy E (we define ΔE = E − E_{c}). For each star, we computed the mean value of I_{s} and its dispersion σ_{Is} from ~10 000 steps along the orbit, and we determined the value of z_{0} that minimizes σ_{Is}. To ensure that σ_{s} is determined sufficiently well, the total time integration for each orbit is about 300 dynamical times (i.e. 300 galactic rotations). The numerical integration is a RungeKuttaFehlberg of order 7(8) (Fehlberg 1968). We find that stellar orbits with the same values of L_{z} and ΔE all have σ_{Is} minimum for similar values of z_{0}. For different energies and angular momenta, the adjusted z_{0} will be different.
To allow an easier presentation of our results, we normalized I_{s} as ${\mathit{I}}_{\mathrm{3}}\mathrm{=}\mathrm{}\frac{{\mathit{I}}_{\mathrm{s}}}{\mathrm{\Delta}\mathit{E}}{\left(\mathrm{1}\mathrm{+}\frac{{\mathit{R}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\left(}{\mathit{L}}_{\mathit{z}}\mathrm{\right)}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\right)}^{1}\mathit{,}$(4)still an integral of motion. I_{3} = 0 corresponds to orbits confined within the galactic plane, while I_{3} ~ 1 corresponds to shell orbits with minimum radial variations when they cross the plane z = 0.
In the following sections, we consider three potentials: a logarithmic potential with a disc component, a flattened logarithmic potential, and a potential similar to the Besançon Galactic model potential. In all the three cases, the constancy of the third integral is satisfied at 0.2 to 0.5 per cent for orbits corresponding to thin or thick stars and to a 1 to 2 per cent for orbits with larger vertical oscillations.
3.1. Logarithmic potential of MyamotoNagai type
The logarithmic potential of MyamotoNagai type $\mathrm{\Phi}\mathrm{=}{\mathit{v}}_{\mathrm{\infty}}^{\mathrm{2}}\mathrm{log}{\left({\mathit{R}}^{\mathrm{2}}\mathrm{+}{\left(\mathit{a}\mathrm{+}\sqrt{{\mathit{b}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}$(5)has a positive density (Zotos 2011). It has a spherical halo and a disc component (with a + b the halo core radius and b the disc thickness). This potential presents analogy with the family of densitypotential pairs of logarithmic MyamotoNagai type given by Bienaymé (2009).
We consider the orbits with the angular momentum L_{z} = 8.5 (v_{∞} = 1, a = 0.5, b = 0.2, R_{c} = 8.5) and with E = E_{c} + ΔE. The initial conditions of the computed orbits are v_{R} = 0, with R_{initial} equally spaced. These initial conditions do not include resonant orbits that do not cross the galactic plane perpendicularly.
Table 1 summarizes the results. For each energy ΔE, the adjusted z_{0} that minimizes σ_{I3} is given. The maximum vertical extension z_{max} corresponds to an orbit with I_{3} close to 1. The dispersions of I_{3} always remain small, the mean dispersion being smaller than one per cent. However, for some resonant orbits, the dispersion is larger, of the order of 2 per cent. We also note that at large energies corresponding to halo stars the variation of I_{3} still remains small. This must be linked to the spherical shape of the potential at large z, easily modelled by a Stäckel potential.
Dischalo logarithmic potential.
In the next section we consider a different requirement using a flattened potential.
Fig. 1 Top panel: meridional projection for three orbits within the logarithmic potential (Sect. 3.2) with L_{z} = 8.5 and ΔE = 0.2. Zero velocity curve and some coordinate curves of the elliptic coordinates (with z_{0} = 5.9) are drawn. Middle panel: analytical determination of the envelopes of the orbits. Bottom panel: surfaces of section for the same three orbits. Black crosses: numerically computed orbits. Red lines: the corresponding sections obtained from Eqs. (1)−(3). Blue line: surface of section of the orbit confined in the midplane (i.e. v_{z} = 0). 
Fig. 2 Left: quasi integral I_{3} versus z_{max} for stars with ΔE = 6400 and R_{c}(L_{z}) = 8500 pc for the BGM potential. Errors bars are σ_{I3}. Right: idem with ΔE = 12 800. 
Fig. 3 Left: dispersion σ_{I3} of the integral I_{3} along orbits versus z_{max} for stars with ΔE = 6400 and R_{c}(L_{z}) = 8500 for the BGM potential. Right: idem with ΔE = 12 800. 
3.2. Logarithmic potential
To examine the impact of the flattening of the dark halo on the dispersions σ_{I3}, we consider the traditional logarithmic potential Richstone (1980): $\mathrm{\Phi}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{log}\mathrm{(}{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathit{/}{\mathit{q}}^{\mathrm{2}}\mathrm{)}\mathit{.}$(6)We set v_{c} = 1 and q = 0.8 to correspond to a flattening of the isodensity q_{ρ} = 0.65.
We computed orbits with L_{z} = 8.5 (R_{c} = 8.5, E_{c} = 1). Table 2 summarizes the results. Surfaces of section for L_{z} = 1, ΔE = 0.05, and 0.5 are shown in Fig. 9 of Bienaymé & Traven (2013), and we note the lack of resonant orbits. This may explain the excellent conservation of I_{3} at a level of 0.2 per cent for stellar disc orbits. Orbits with larger extension up to z_{max} = 9 (E = 0.4) have I_{3} conserved at 2 per cent.
The inspection (see Fig. 1) of the surface of section and of the meridional projection of three orbits (ΔE = 0.2 and L_{z} = 8.5) allows visualization of the relative agreement between the nearly exact numerical orbits and the contours obtained from the approximated third integral. For the meridional projection, the red dashed lines are ellipsoids and hyperboloids from the ellipsoidal coordinates with z_{0} = 5.9 (lines are drawn from the corner of the orbit envelopes). The blue continuous lines in Fig. 1 (middle panel) are the envelopes that are determined semianalytically from E, L_{z} and the quasi integral I_{3} that are explicitly known: we determine surfaces of section at various fixed R (resp. fixed z) and find the z (resp. R) extension limits of the orbit from the condition ∂I_{3}/∂v_{r} = 0 (resp. ∂I_{3}/∂v_{z} = 0). We note small differences between the envelopes of numerical orbits (top panel) and semianalytical envelopes (middle panel). We also note that the analytical hypersurfaces defined from E, I_{3} do not have the topology of a torus everywhere. The bottom panel of Fig. 1 shows a section of the phase space (z = 0) for the same three orbits from the numerically computed orbits and from the analytical section given by Eq. (1). If v_{c} = 200 km s^{1}, the maximum v_{R} of these three orbits is ~120, 50, and 20 km s^{1}.
3.3. The Besançon Galactic mass model
We apply our test to the gravitational potential of the Besançon Galactic model (BGM; Bienaymé et al. 1987; Robin et al. 2003; Czekaj et al. 2014), which is a more realistic one than the two previous ones, although we limited our analysis to its axisymmetric version. This model has a dynamical consistency at the solar Galactic radius position in the sense that the thickness of the stellar components are constrained by the potential and by the vertical velocity dispersions. However, outside of the solar neighbourhood, the vertical kinematics and the thickness of stellar components are constrained by observations and not by a dynamical consistency of the model. The motivation of this work is to improve the BGM dynamical consistency by using stationary distribution functions to model the stellar kinematics. Besides the integrals E and L_{z}, we therefore require an approximate third integral to describe the stellar kinematics.
We determined the gravitational potential of the Besançon Galactic model according to the characteristics of its components described in Robin et al. (2003) and Czekaj et al. (2014). The dark matter component is represented by a spherical component. The orbits with large vertical extensions have shapes that are mainly determined by the spherical halo. Since the exact shape of the dark halo remains uncertain in practice and is a controversial question, we consider a less favourable situation using a flattened dark matter component.
The dark matter halo density used in the BGM is ${\mathit{\rho}}_{\mathrm{dm}}\mathrm{~}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{core}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{\xb7}$(7)We flatten the halo, replacing the quantity R with $\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}\mathit{/}{\mathit{q}}^{\mathrm{2}}}$ in the expression of the potential. With R_{core} = 2700 pc and q = 0.8, the flattening of the density is q_{ρ} = 0.62. (For this value of q, negative densities only occur in a region far from our domain of interest.)
We also modify the pointmass bulge (Bienaymé et al. 1987) using the potential law ${\mathrm{\Phi}}_{\mathit{B}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{{\mathit{R}}_{\mathrm{bulbe}}^{\mathrm{2}}\mathrm{+}{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$ with R_{bulbe} = 2000 pc.
Logarithmic potential (see legend of Table 1).
We summarize the results for orbits with L_{z} corresponding to R_{c} = 8500 pc in Table 3. Results are listed by families of different energies ΔE. Figure 2 shows I_{3} versus z_{max}, the maximum vertical extension of orbits in cases ΔE = 6400 or 12 800 and R_{c} = 8500 pc. Error bars are the dispersion of I_{3} along each orbit. The dispersion of I_{3} is zero for orbits confined in the midplane. Shell orbits have I_{3} maximum and σ_{I3} close to a few 10^{4} (Fig. 3). This is unexpected since there is no reason that the potential along the “thin” orbits should have exactly the Stäckel form. On the other hand, when z_{max} ~ 2000 pc, the dispersion of I_{3} is maximum and increases sharply for all energies (see Fig. 3), also corresponding in Fig. 2 to a change of shape of the distribution of points. This is related to the presence of (1, 1) resonant orbits at this vertical height that are poorly modelled by our Stäckel adjustment. For these orbits the dispersion of the quasi integral is less than 5 per cent. The median dispersion of I_{3} for orbits with z_{max} smaller than 2 kpc remains very small ~0.002.
We performed the Stäckel adjustment for other galactic radii R_{c} from 1.5 to 15 kpc. For each pair of values (E = Ec + ΔE, L_{z}), we determine the z_{0} value that minimizes σ_{I3}. Figure 4 plots the fitted z_{0} (crosses) for ΔE = 100, 200, 400...12 800. Lines of constant ΔE are plotted. We obtain a tabulation of z_{0} that gives us a function of two integrals of motion z_{0}(E,L_{z}) that can be used to build a distribution function for stellar populations.
The median dispersions of I_{3} remain low except at galactic radii smaller than 3 kc where many resonant orbits dominate the phase space. Figure 5 plots the histograms of σ_{I3} for R_{c} from 3.5 to 15.5 kpc.
Fig. 4 BGM potential: z_{0} best fit versus R_{c}(L_{z}) and ΔE. Black lines ΔE = 100 and 200; red lines: 400 and 800; green lines: 1600 and 3200; blue lines: 64 000 and 128 000. 
Fig. 5 Histograms of I_{3} dispersion along orbits for z_{max} intervals delimited by 0, 1000 pc, 2000 pc, and beyond (respectively black, dotted red, dashed blue lines). For each of the three z_{max} intervals, the medians of σ_{I3} are respectively 5 × 10^{4}, 1.5 × 10^{3}, and 5 × 10^{3}. 
4. Distribution function and Jeans equations
Fig. 6 Left: dark lines: radial forces K_{R} of the BGM at several z above the galactic plane (0.5, 1, 2, and 3 kpc, bottom to top). Thin green lines: K_{R} recovered from a Jeans equation in the case of a thick disc. Thin red lines: the case of a thin disc. Right: relative errors on the recovered radial force K_{R} versus R at various z (0.5, 1, 2, and 3 kpc; resp. green, black, red, and blue lines) for thin disc DF (dotted lines) and thick disc DF (continuous lines). 
We tested the efficiency of using the quasi integral I_{3} to model the distribution function of disc stars using a Shu distribution function (DF) generalized to 3D axisymmetric potentials. It writes as (here, we correct a typographic error in Bienaymé 1999): $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{E,}{\mathit{L}}_{\mathit{z}}\mathit{,}{\mathit{I}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}& & \frac{\mathrm{2}\mathrm{\Omega}\mathrm{\left(}{\mathit{R}}_{\mathrm{c}}\mathrm{\right)}}{\mathrm{2}\mathit{\pi \kappa}\mathrm{\left(}{\mathit{R}}_{\mathrm{c}}\mathrm{\right)}}\frac{\mathrm{\Sigma}\mathrm{\left(}{\mathit{L}}_{\mathit{z}}\mathrm{\right)}}{{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}}\mathrm{exp}\left[\mathrm{}\frac{\mathrm{(}\mathit{E}\mathrm{}{\mathit{E}}_{\mathrm{circ}}\mathrm{)}}{{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}}\right]\\ & & \mathrm{\times}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi}}}\frac{\mathrm{1}}{{\mathit{\sigma}}_{\mathit{z}}}\mathrm{exp}\left[\mathrm{}\mathrm{(}\mathit{E}\mathrm{}{\mathit{E}}_{\mathrm{circ}}\mathrm{)}\left(\frac{\mathrm{1}}{{\mathit{\sigma}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{}\frac{\mathrm{1}}{{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}}\right){\mathit{I}}_{\mathrm{3}}\right]\mathit{.}\end{array}$(8)with R_{c}(L_{z}) the radius of the circular orbit with the angular momentum L_{z}, Ω the angular velocity, κ the epicyclic frequency, and E_{circ} the energy of a circular orbiting star at radius R_{c}. For sufficiently small velocity dispersions, the number density distribution, Σ(L_{z}) = Σ_{0}exp(−R_{c}/R_{ν}), is close to Σ(R) = Σ_{0}exp(−R/R_{ν}). We set constant the parameters σ_{R,z}(L_{z}) and R_{ν} = 2.5 kpc, which is close to the scale length of the number density distribution. The DF allows us to reproduce the triaxiality and tilt of the velocity ellipsoid and to model nearly exponential density disc distribution. It could be easily modified to reproduce any reasonable radial density law. We set a constant velocity dispersion (R_{σR,z} = ∞) for a thin disc (σ_{R},σ_{z}) = (40 km s^{1}, 20 km s^{1}) and for a thick disc (σ_{R},σ_{z}) = (60 km s^{1}, 40 km s^{1}).
The distribution function can also be written in a more readable form as $\mathit{f}\mathrm{\left(}\mathit{E,}{\mathit{L}}_{\mathit{z}}\mathit{,}{\mathit{I}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\mathit{g}\mathrm{\left(}{\mathit{L}}_{\mathit{z}}\mathrm{\right)}\mathrm{exp}\left[\mathrm{}\frac{{\mathrm{\mathcal{E}}}_{\mathit{R}}}{{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}}\right]\mathrm{exp}\left[\mathrm{}\frac{{\mathrm{\mathcal{E}}}_{\mathit{z}}}{{\mathit{\sigma}}_{\mathit{z}}^{\mathrm{2}}}\right]\mathit{,}$(9)where ℰ_{R} and ℰ_{z} are integrals of motion that can be easily deduced from Eq. (8): $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathrm{\mathcal{E}}}_{\mathit{R}}\mathrm{=}\mathrm{(}\mathit{E}\mathrm{}{\mathit{E}}_{\mathrm{circ}}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{}{\mathit{I}}_{\mathrm{3}}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathrm{\mathcal{E}}}_{\mathit{z}}\mathrm{=}\mathrm{(}\mathit{E}\mathrm{}{\mathit{E}}_{\mathrm{circ}}\mathrm{)}\hspace{0.17em}{\mathit{I}}_{\mathrm{3}}\mathit{.}\end{array}$(10)They are respectively related to the amount of radial or vertical motions that can be controlled versus the parameters σ_{R} and σ_{z}. Within a Stäckel potential and for an orbit with angular momentum L_{z}, they are respectively at position (R_{c}(L_{z}),z = 0) the radial and vertical kinetic energies. Shell orbits have ℰ_{R} = 0.
We determine the moments of the distribution function (Eq. (8)) and recover the radial and vertical forces, K_{R} and K_{z}, from the Jeans equations for a stationary axisymmetric potential. The moments are linked through the Jeans equations where the time derivative are not exactly zero since I_{3} is just an approximate integral: $\begin{array}{ccc}& & \frac{\mathit{\partial}\hspace{0.17em}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial R}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}^{\mathrm{2}}}\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial z}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}{\mathit{v}}_{\mathit{z}}}\mathrm{\right)}\mathrm{+}\mathit{\nu}\left(\frac{\overline{){\mathit{v}}_{\mathit{R}}^{\mathrm{2}}}\mathrm{}\overline{){\mathit{v}}_{\mathit{\phi}}^{\mathrm{2}}}}{\mathit{R}}\mathrm{+}\frac{\mathit{\partial}\mathrm{\Phi}}{\mathit{\partial R}}\right)\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial}\hspace{0.17em}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{z}}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial R}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}{\mathit{v}}_{\mathit{z}}}\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial z}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{\right)}\mathrm{+}\mathit{\nu}\left(\frac{\overline{){\mathit{v}}_{\mathit{R}}{\mathit{v}}_{\mathit{z}}}}{\mathit{R}}\mathrm{+}\frac{\mathit{\partial}\mathrm{\Phi}}{\mathit{\partial z}}\right)\mathrm{=}\mathrm{0.}\end{array}$By considering that I_{3} is close to an exact integral and neglecting the time derivative terms, we transform the Jeans equations as $\begin{array}{ccc}& & \frac{\mathit{\partial}}{\mathit{\partial R}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}^{\mathrm{2}}}\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial z}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}{\mathit{v}}_{\mathit{z}}}\mathrm{\right)}\mathrm{+}\mathit{\nu}\left(\frac{\overline{){\mathit{v}}_{\mathit{R}}^{\mathrm{2}}}\mathrm{}\overline{){\mathit{v}}_{\mathit{\phi}}^{\mathrm{2}}}}{\mathit{R}}\mathrm{+}{\left(\frac{\mathit{\partial}\mathrm{\Phi}}{\mathit{\partial R}}\right)}_{\mathrm{est}\mathit{.}}\right)\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial}}{\mathit{\partial R}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{R}}{\mathit{v}}_{\mathit{z}}}\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial z}}\mathrm{\left(}\mathit{\nu}\overline{){\mathit{v}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{\right)}\mathrm{+}\mathit{\nu}\left(\frac{\overline{){\mathit{v}}_{\mathit{R}}{\mathit{v}}_{\mathit{z}}}}{\mathit{R}}\mathrm{+}{\left(\frac{\mathit{\partial}\mathrm{\Phi}}{\mathit{\partial z}}\right)}_{\mathrm{est}\mathit{.}}\right)\mathrm{=}\mathrm{0.}\end{array}$Thus, under the assumption that the approximate integral I_{3} is exact, Eqs. (13), (14) are used to obtain an estimate of the radial and vertical forces, ${\mathit{K}}_{\mathit{R}}\mathrm{\approx}\mathrm{}{{}^{\mathrm{\right(}}\frac{\mathit{\partial \phi}}{\mathit{\partial R}}^{\mathrm{\left)}}}_{\mathrm{est}\mathit{.}}$${\mathit{K}}_{\mathit{z}}\mathrm{\approx}\mathrm{}{{}^{\mathrm{\right(}}\frac{\mathit{\partial \phi}}{\mathit{\partial z}}^{\mathrm{\left)}}}_{\mathrm{est}\mathit{.}}$. They can be compared to the exact forces to test the efficiency to recover the galactic potential using our Stäckel integral in a nonStäckel potential. Figure 6 shows the radial forces K_{R} of the BGM potential and its relative error at several z when recovered from Eqs. (13), (14).
Within the domain R = 3 to 16 kc and z smaller than 3 kpc, the relative error on the radial force K_{R} is smaller than one per cent for R> 6 kpc. The relative error is ten per cent at R = 3 kpc and increases below R = 3 kpc.
Figure 7 shows the exact and recovered K_{z} force versus z at several R (3.5, 4.5, 6.5, 8.5, 12.5). At large R> 6 kpc, the K_{z} forces are accurately recovered for the thick disc up to 6 kpc and up to 3 kpc for the thin disc. The vertical force K_{z} at the solar position R = 8 kpc is remarkably well recovered at 0.5% up to 6 kpc for the thick disc (7 scale heights; Fig. 8). For the thin disc, the accuracy is lost above 3 kpc (corresponding however to ten scale heights of the thin disc). At lower radius R from 3.5 to 5, the K_{z} force is only approximately recovered up to z = 2 kpc.
For comparison (Fig. 8), we also plot the K_{z} force recovered neglecting the cross term ⟨ v_{R}v_{z} ⟩ in Jeans equations, a classical assumption valid at low z. We note that this hypothesis remains valid up to 500 pc for the thin disc. At higher z, the recovered force diverges quickly from the exact one.
Fig. 7 Vertical force K_{z} versus z at R = (3.5, 4.5, 6.5, 8.5, 12.5 kpc) (top to bottom) from the BGM (black lines) and recovered K_{z} for a thick disc DF (blue dotted lines). 
Fig. 8 Vertical force K_{z} versus z at R = 8500 pc (black line) and recovered K_{z} from simplified Jeans equation assuming the separability of vertical and radial motions: red dotted line for a thin disc DF and blue dashed lines for a thick disc DF. 
5. Collisionless Boltzmann equation
The stationarity of the Jeans equation is a necessary condition for validating the degree of stationarity of a distribution function, but this is not a sufficient condition. We should also have to test the stationarity of all the other moments of the collisionless Boltzmann equation (CBE), since it is easy to find solutions for the Jeans equations that are not solutions of the CBE, so the Jeans equation test can be more optimistic than the indications solely obtained from the conservation of the integral of motion. For this reason we examine hereafter the stationarity of the CBE.
5.1. First estimate
We estimate the stationarity of the distribution function built with the approximated third integral by looking at the time variation of the DF (Eq. (8)) over a dynamical time (~ an orbit revolution). We have $\frac{\mathrm{d}\mathrm{ln}\mathit{f}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{\partial}\mathrm{ln}\mathit{f}}{\mathit{\partial t}}\mathrm{+}\left[\mathrm{ln}\mathit{f,H}\right]\mathrm{=}\mathrm{0}$(15)or $\frac{\mathit{\partial}\mathrm{ln}\mathit{f}}{\mathit{\partial t}}\mathrm{}\frac{\mathrm{\Delta}\mathit{E}}{{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}}\left(\frac{{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}}{{\mathit{\sigma}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{}\mathrm{1}\right)\frac{\mathrm{d}{\mathit{I}}_{\mathrm{3}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{0}\mathit{,}$(16)and the relative variation of f over a dynamical or longer time is determined by the variation σ_{I3} of the quasi integral I_{3} along orbits: $\left\frac{\mathrm{\Delta}\mathit{f}}{\mathit{f}}\right\mathrm{~}\frac{\mathrm{\Delta}\mathit{E}}{{\mathit{\sigma}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{\sigma}}_{\mathit{I}\mathrm{3}}\mathit{.}$(17)Orbits with small radial or vertical amplitudes (ΔE and σ_{I3} small) have the smallest variations in f. Thus the thinnest discs are the most accurately modelled.
The stationarity decreases quickly with z since we have $\mathrm{\Delta}\mathit{E}\mathit{/}{\mathit{\sigma}}_{\mathit{R}}^{\mathrm{2}}\mathrm{\propto}{\mathit{z}}^{\mathrm{~}\mathrm{0.5}}$. It also decreases because of the significant presence of resonant orbits between 1 to 2 kpc. However, at the solar Galactic radius R_{0}, orbits within the BGM with large vertical amplitude remain correctly modelled with a Stäckel potential thanks to small σ_{I3} of the order of 0.005. In this situation, the distribution function remains stationary at many scale heights above the galactic plane. For the thick disc, the stationarity (Eq. (17)) is about one per cent at six scale heights (z ~ 6 kpc). For the thin disc, it becomes insufficient at about ten scale heights (z ~ 3 kpc).
The time derivative in Eq. (15) is the time variation from an initial position of stars in phase space at t = 0 given by the initial condition of Eq. (9). How divergent are the orbits evolving from this supposedly nearequilibrium initial condition? We know that for potentials similar to the three potentials considered in this paper, the orbits are essentially regular and an effective third integral must exist, which could be eventually estimated with a higher accuracy using, for instance, highorder polynomials (Bienaymé & Traven 2013) or a torus fitting (Sanders & Binney 2014). This would, however, not be true in cases of significant presence of ergodicity, for instance close to the corotation of a barred rotating potential. It turns out that our approximate integrals oscillate along each orbit around a mean value, that a typical periodicity of these oscillations is the dynamical time, and that the amplitude of these oscillations is of the order of σ_{I3}. This is illustrated in Fig. 9 where the maxima of dI_{3}/ dt are of the order of σ_{I3}/t_{dyn} with t_{dyn} ~ 6.
Fig. 9 dI_{3}/ dt: time derivative of the quasi integral within the logarithmic potential for the three orbits shown in Fig. 1: (black I_{3} = 0.16, red I_{3} = 0.88, blue I_{3} = 1.06). 
5.2. Second estimate
It is also useful to consider the nonstationarity as a shift in positions or velocities of the modelled DF relative to a stationary DF. A simple 1D dimensional analysis allows it to be illustrated. With a quadratic potential Φ = az^{2} and a stationary DF as ${\mathit{f}}_{\mathrm{stat}\mathit{.}}\mathrm{=}\mathrm{exp}\mathrm{(}\mathrm{}\mathrm{(}\mathit{\phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{)}\mathit{/}{\mathit{\sigma}}^{\mathrm{2}}\mathrm{)}\mathit{,}$(18)the shifted DF in velocity by a factor v_{0} is $\mathit{f}\mathrm{=}\mathrm{exp}\left[\mathrm{}\frac{\mathrm{1}}{{\mathit{\sigma}}^{\mathrm{2}}}\left(\mathit{\phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{+}\frac{\mathrm{(}{\mathit{v}}_{\mathit{z}}\mathrm{}{\mathit{v}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}}\right)\right]\mathrm{\xb7}$(19)Such a DF could be considered satisfying if the shift is small even if we notice that the relative errors on f increase in the tails of the distribution.
We deduce from the CBE, $\frac{\mathit{\partial}\mathrm{ln}\mathit{f}}{\mathit{\partial t}}\mathrm{=}\mathrm{}\frac{{\mathit{K}}_{\mathit{z}}}{{\mathit{\sigma}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\mathit{,}$(20)which combined with Eq. (17), gives the variation in f over a dynamical time t_{dyn}: $\frac{\mathrm{\left}{\mathit{K}}_{\mathit{z}}\mathrm{\right}}{{\mathit{\sigma}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{v}}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{\Delta}\mathit{E}}{{\mathit{\sigma}}^{\mathrm{2}}}\hspace{0.17em}\frac{{\mathit{\sigma}}_{{\mathit{I}}_{\mathrm{3}}}}{{\mathit{t}}_{\mathrm{dyn}}}\mathrm{\xb7}$(21)With a quadratic potential and h, the scale height of the disc population, we estimate the velocity shift: ${\mathit{v}}_{\mathrm{0}}\mathrm{~}\frac{\mathrm{\Delta}\mathit{E}}{\mathrm{\left}{\mathit{K}}_{\mathit{z}}\mathrm{\right}}\hspace{0.17em}\frac{{\mathit{\sigma}}_{\mathit{I}\mathrm{3}}}{{\mathit{t}}_{\mathrm{dyn}}}\mathrm{~}\frac{\mathit{z}}{\mathit{h}}\mathit{\sigma}\hspace{0.17em}{\mathit{\sigma}}_{\mathit{I}\mathrm{3}}\mathit{.}$(22)For the thin disc DF at one scale height (h ~ 250 pc), the shift v_{0} is 0.2 km s^{1} and at 4 kpc (16 scale heights), it is 10 km s^{1} (to be compared to the 20 km s^{1} of the vertical velocity dispersions) where the application of the Jeans equation shows that the DF is not stationary. At one scale height of ~1 kpc for the thick disc, the shift is v_{0} ~ 1 km s^{1}. At 6 kpc, it is 7 km s^{1} (to be compared to the 40 km s^{1} of the vertical velocity dispersion), and the Jeans equation still validates the stationarity.
6. Conclusion
In this paper we have shown that Stäckel potentials can be used to fit a wide variety of disc stellar orbits within different axisymmetric potentials of disc galaxies. We also proposed a new and simple formulation for the third integral of motion of Galactic potentials, explicitly depending on the potential, which is an integral known to be exact in the case of Stäckel potentials.
The quality of the fit of orbits was quantitatively measured by looking at the conservation of the approximate third integral of motion besides the energy and angular momentum. By using the Besançon Galactic model, we showed that the third integral is conserved to a few thousandths in a wide volume around the solar neigbourhood. It is conserved to better than one per cent at 6 kpc above the Galactic plane at the solar position. However, the fit fails at low galactic radius (R smaller than 4 kpc) owing to the presence of a large number of resonant orbits.
In light of the future Gaia data, we also used distribution functions of disc stars, depending on the three integrals, and checked the ability of Jeans equation to recover the gravitational potential. We also considered the stationarity of the CBE.
In conclusion, Stäckel potentials can be good local approximations of general realistic potentials. They have been used many times to fit the global potential of galaxies (see for instance Dejonghe & de Zeeuw 1988) or the potential of our own Galaxy (Famaey & Dejonghe 2003). The fit of stellar orbits are also frequently used using Stäckel potentials with local or global fitting (Kent & de Zeeuw 1991). We showed here that they be can used very efficiently to build distribution function of disc stellar populations. A straightforward application will consist in using such modelling of DFs and in extending the dynamical consistency of the Besançon Galactic model to describe the kinematics of disc stellar populations.
References
 Batsleer, P., & Dejonghe, H. 1994, A&A, 287, 43 [NASA ADS] [Google Scholar]
 Bienaymé, O. 1999, A&A, 341, 86 [NASA ADS] [Google Scholar]
 Bienaymé, O. 2009, A&A, 500, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bienaymé, O., & Traven, G. 2013, A&A, 549, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bienaymé, O., Robin, A. C., & Crézé, M. 1987, A&A, 186, 359 [NASA ADS] [Google Scholar]
 Bienaymé, O., Famaey, B., Siebert, A., et al. 2014, A&A, 571, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J. 2012, MNRAS, 426, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Colombi, S., Sousbie, T., Peirani, S., Plum, G., & Suto, Y. 2015, MNRAS, in press [arXiv:1504.07337] [Google Scholar]
 Czekaj, M. A., Robin, A. C., Figueras, F., Luri, X., & Haywood, M. 2014, A&A, 564, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 De Bruyne, V., Leeuwin, F., & Dejonghe, H. 2000, MNRAS, 311, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Dejonghe, H., & de Zeeuw, T. 1988, ApJ, 333, 90 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, T. 1985, MNRAS, 216, 273 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, P. T., & LyndenBell, D. 1985, MNRAS, 215, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & Dejonghe, H. 2003, MNRAS, 340, 752 [NASA ADS] [CrossRef] [Google Scholar]
 Fehlberg, E. 1968, NASA technical report TR R287 [Google Scholar]
 Hori, G. 1962, PASJ, 14, 353 [NASA ADS] [Google Scholar]
 Kent, S. M., & de Zeeuw, T. 1991, AJ, 102, 1994 [NASA ADS] [CrossRef] [Google Scholar]
 Manabe, S. 1979, PASJ, 31, 369 [NASA ADS] [Google Scholar]
 LyndenBell, D. 1962, MNRAS, 124, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Ollongren, A. 1962, Bull. Astron. Inst. Netherlands, 16, 241 [Google Scholar]
 Papaphilippou, Y., & Laskar, J. 1998, A&A, 329, 451 [NASA ADS] [Google Scholar]
 Piffl, T., Binney, J., McMillan, P. J., et al. 2014, MNRAS, 445, 3133 [NASA ADS] [CrossRef] [Google Scholar]
 Renaud, F., Bournaud, F., Emsellem, E., et al. 2013, MNRAS, 436, 1836 [NASA ADS] [CrossRef] [Google Scholar]
 Richstone, D. O. 1980, ApJ, 238, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sanders, J. L., & Binney, J. 2014, MNRAS, 441, 3284 [NASA ADS] [CrossRef] [Google Scholar]
 Syer, D., & Tremaine, S. 1996, MNRAS, 282, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Valluri, M., & Merritt, D. 1998, ApJ, 506, 686 [NASA ADS] [CrossRef] [Google Scholar]
 van de Hulst, H. C. 1962, Bull. Astron. Inst. Netherlands, 16, 235 [NASA ADS] [Google Scholar]
 Wayman, P. A. 1959, MNRAS, 119, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshikawa, K., Yoshida, N., & Umemura, M. 2013, ApJ, 762, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Zotos, E. E. 2011, New Astron., 16, 391 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Quasi integral of motion
We refer the reader to de Zeeuw (1985) and de Zeeuw & LyndenBell (1985) for notations and for a detailed description of properties of Stäckel potentials.
An axisymmetric Stäckel potential is fully defined by two free functions h(λ) and h(ν). In the case of prolate spheroidal coordinates, ± z_{0} (${\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma}\mathrm{}\mathit{\alpha}$ ) are the foci of the confocal ellipsoids and hyperboloids used to define a system of coordinates (λ,ν) in which Stäckel potentials are more easily tractable.
Let $\mathit{V}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}\begin{array}{c}{}_{\mathrm{\u02dc}}\\ \mathit{V}\end{array}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}$ be the true potential that we approximate locally at the position (R_{1},z_{1}) = (λ_{1},ν_{1}) with a Stäckel potential Φ(λ,ν). We assume that z_{0} is already known.
By definition, we have for a Stäckel potential: $\mathrm{\Phi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{=}\mathrm{}\frac{\mathit{h}\mathrm{\left(}\mathit{\lambda}\mathrm{\right)}\mathrm{}\mathit{h}\mathrm{\left(}\mathit{\nu}\mathrm{\right)}}{\mathit{\lambda}\mathrm{}\mathit{\nu}}\mathit{,}$(A.1)and it is always possible to find a Säckel potential that coincides with any potential on the chosen coordinate surfaces λ = λ_{1} and ν = ν_{1}. It writes as $\mathrm{\Phi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{=}\frac{\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}\mathit{\lambda ,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{\right)}\mathrm{(}\mathit{\lambda}\mathrm{}{\mathit{\nu}}_{\mathrm{1}}\mathrm{)}\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{\right)}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathrm{}{\mathit{\nu}}_{\mathrm{1}}\mathrm{)}\mathrm{+}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,\nu}\mathrm{\right)}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathrm{}\mathit{\nu}\mathrm{)}}{\mathit{\lambda}\mathrm{}\mathit{\nu}}\mathrm{\xb7}$(A.2)Thus $\begin{array}{ccc}\mathit{h}\mathrm{\left(}\mathit{\lambda}\mathrm{\right)}& \mathrm{=}& \begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}\mathit{\lambda ,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{\right)}\mathrm{(}\mathit{\lambda}\mathrm{}{\mathit{\nu}}_{\mathrm{1}}\mathrm{)}\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{\right)}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathrm{}{\mathit{\nu}}_{\mathrm{1}}\mathrm{)}\mathrm{+}\mathit{C}\\ \mathit{h}\mathrm{\left(}\mathit{\nu}\mathrm{\right)}& \mathrm{=}& \mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,\nu}\mathrm{\right)}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathrm{}\mathit{\nu}\mathrm{)}\mathrm{+}\mathit{C}\end{array}$(A.3)and C an arbitrary constant.
We can express the third integral I_{s} associated to Φ as ${\mathit{I}}_{\mathrm{s}}\mathrm{=}\mathrm{\Psi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{z}}^{\mathrm{2}}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}{\mathit{v}}_{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\theta}}^{\mathrm{2}}\mathrm{)}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{R}}^{\mathrm{2}}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\right)\hspace{0.17em}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{+}\frac{\mathit{R}\hspace{0.17em}\mathit{z}\hspace{0.17em}{\mathit{v}}_{\mathit{R}}\hspace{0.17em}{\mathit{v}}_{\mathit{z}}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}$(A.4)or ${\mathit{I}}_{\mathrm{s}}\mathrm{=}\mathrm{\Psi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{L}}^{\mathrm{2}}\mathrm{}{\mathit{L}}_{\mathit{z}}^{\mathrm{2}}}{{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}$(A.5)with $\mathrm{\Psi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathit{\nu}\mathrm{+}\mathit{\gamma}\mathrm{)}\hspace{0.17em}\mathit{h}\mathrm{\left(}\mathit{\lambda}\mathrm{\right)}\mathrm{}\mathrm{(}\mathit{\lambda}\mathrm{+}\mathit{\gamma}\mathrm{)}\hspace{0.17em}\mathit{h}\mathrm{\left(}\mathit{\nu}\mathrm{\right)}}{\mathrm{(}\mathit{\gamma}\mathrm{}\mathit{\alpha}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathit{\lambda}\mathrm{}\mathit{\nu}\mathrm{)}}\mathrm{\xb7}$(A.6)We fix the remaining free constant C by setting h(ν = −γ) = 0, so Ψ is null at z = 0 in the plane of symmetry of the potential. Then, evaluated at (λ_{1},ν_{1}), this function simplifies as $\mathrm{\Psi}\mathrm{\left(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}\mathrm{}\frac{\mathrm{\left(}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\right(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{)}\mathrm{}\begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}\mathrm{}\mathit{\gamma}{\mathrm{\left)}}^{\mathrm{\right)}}\hspace{0.17em}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathrm{+}\mathit{\gamma}\mathrm{)}}{\mathit{\gamma}\mathrm{}\mathit{\alpha}}\mathrm{\xb7}$(A.7)If we set $\mathit{\gamma}\mathrm{=}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}$ and α = 0,
$\begin{array}{ccc}& & \begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{\left(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}{\mathit{\nu}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}\mathit{V}\mathrm{\left(}{\mathit{R}}_{\mathrm{1}}\mathit{,}{\mathit{z}}_{\mathrm{1}}\mathrm{\right)}\mathit{,}\\ & & \begin{array}{c}\mathrm{\u02dc}\\ \mathit{V}\end{array}\mathrm{(}{\mathit{\lambda}}_{\mathrm{1}}\mathit{,}\mathrm{}\mathit{\gamma}\mathrm{)}\mathrm{=}\mathit{V}\mathrm{(}\mathit{R}\mathrm{=}\sqrt{{\mathit{\lambda}}_{\mathrm{1}}}\mathit{,z}\mathrm{=}\mathrm{0}\mathrm{)}\mathit{,}\end{array}$and $\begin{array}{ccc}{\mathit{\lambda}}_{\mathrm{1}}\mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{R}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{(}{\mathit{R}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{R}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{.}& \end{array}$(A.10)Thus, from Eqs. (A.4) and (A.7) (similar to Eq. (2)), we obtain a simple expression for the third integral I_{s} at position (λ_{1},ν_{1}), which is exact in the case of a Stäckel potential and explicitly depends on the potential, the coordinates, and the velocities. We note that, in practice, the intermediate functions h(λ) and h(ν) do not need to be evaluated.
In summary, the proposed quasi integral in Eq. (1) results from the exact third integral of the orbits in the Stäckel potential of Eq. (A.2), which is equal to the true potential on the surfaces λ = λ_{1} and ν = ν_{1}. The integral is reevaluated locally for each point (λ,ν) = (λ_{1},ν_{1}).
The case of Stäckel potentials with oblate spheroidal coordinates leads to the same equations, but with z_{0} imaginary, ${\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}\mathit{<}\mathrm{0}$ and γ<α. In the z_{0} = ∞ limit, a Stäckel potential is separable in R and z, $\begin{array}{ccc}\mathit{V}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}{\mathit{V}}_{\mathrm{1}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{+}{\mathit{V}}_{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{,}& & \end{array}$and $\begin{array}{ccc}\mathrm{}{\mathit{I}}_{\mathrm{s}}\mathrm{=}\left[{\mathit{V}}_{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{}{\mathit{V}}_{\mathrm{2}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\right]\mathrm{+}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{2.}& & \end{array}$Thus the vertical energy is an integral of motion.
In the z_{0} = 0 limit, Stäckel potentials have the form $\begin{array}{ccc}\mathit{V}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}{\mathit{V}}_{\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\frac{{\mathit{V}}_{\mathrm{2}}\mathrm{\left(}\mathit{\theta}\mathrm{\right)}}{{\mathit{r}}^{\mathrm{2}}}\mathit{,}& & \end{array}$where r^{2} = R^{2} + z^{2} and θ is the polar angle (LyndenBell 1962, Table 1), and the integral of motion is $\begin{array}{ccc}\mathrm{}\hspace{0.17em}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{I}}_{\mathrm{s}}\mathrm{\to}\mathrm{\left[}{\mathit{V}}_{\mathrm{2}}\mathrm{\right(}\mathit{\theta}\mathrm{)}\mathrm{}{\mathit{V}}_{\mathrm{2}}\mathrm{(}\mathit{\pi}\mathit{/}\mathrm{2}\mathrm{\left)}\mathrm{\right]}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{L}}^{\mathrm{2}}\mathrm{}{{\mathit{L}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{.}& & \end{array}$
All Tables
All Figures
Fig. 1 Top panel: meridional projection for three orbits within the logarithmic potential (Sect. 3.2) with L_{z} = 8.5 and ΔE = 0.2. Zero velocity curve and some coordinate curves of the elliptic coordinates (with z_{0} = 5.9) are drawn. Middle panel: analytical determination of the envelopes of the orbits. Bottom panel: surfaces of section for the same three orbits. Black crosses: numerically computed orbits. Red lines: the corresponding sections obtained from Eqs. (1)−(3). Blue line: surface of section of the orbit confined in the midplane (i.e. v_{z} = 0). 

In the text 
Fig. 2 Left: quasi integral I_{3} versus z_{max} for stars with ΔE = 6400 and R_{c}(L_{z}) = 8500 pc for the BGM potential. Errors bars are σ_{I3}. Right: idem with ΔE = 12 800. 

In the text 
Fig. 3 Left: dispersion σ_{I3} of the integral I_{3} along orbits versus z_{max} for stars with ΔE = 6400 and R_{c}(L_{z}) = 8500 for the BGM potential. Right: idem with ΔE = 12 800. 

In the text 
Fig. 4 BGM potential: z_{0} best fit versus R_{c}(L_{z}) and ΔE. Black lines ΔE = 100 and 200; red lines: 400 and 800; green lines: 1600 and 3200; blue lines: 64 000 and 128 000. 

In the text 
Fig. 5 Histograms of I_{3} dispersion along orbits for z_{max} intervals delimited by 0, 1000 pc, 2000 pc, and beyond (respectively black, dotted red, dashed blue lines). For each of the three z_{max} intervals, the medians of σ_{I3} are respectively 5 × 10^{4}, 1.5 × 10^{3}, and 5 × 10^{3}. 

In the text 
Fig. 6 Left: dark lines: radial forces K_{R} of the BGM at several z above the galactic plane (0.5, 1, 2, and 3 kpc, bottom to top). Thin green lines: K_{R} recovered from a Jeans equation in the case of a thick disc. Thin red lines: the case of a thin disc. Right: relative errors on the recovered radial force K_{R} versus R at various z (0.5, 1, 2, and 3 kpc; resp. green, black, red, and blue lines) for thin disc DF (dotted lines) and thick disc DF (continuous lines). 

In the text 
Fig. 7 Vertical force K_{z} versus z at R = (3.5, 4.5, 6.5, 8.5, 12.5 kpc) (top to bottom) from the BGM (black lines) and recovered K_{z} for a thick disc DF (blue dotted lines). 

In the text 
Fig. 8 Vertical force K_{z} versus z at R = 8500 pc (black line) and recovered K_{z} from simplified Jeans equation assuming the separability of vertical and radial motions: red dotted line for a thin disc DF and blue dashed lines for a thick disc DF. 

In the text 
Fig. 9 dI_{3}/ dt: time derivative of the quasi integral within the logarithmic potential for the three orbits shown in Fig. 1: (black I_{3} = 0.16, red I_{3} = 0.88, blue I_{3} = 1.06). 

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.