Envelopes and vertical amplitudes of disk-crossing orbits

TL;DR

This paper develops a unified formula for the envelopes of disk-crossing orbits in galactic disks, valid across all vertical amplitudes, linking previous models and relying only on observable surface density.

Contribution

It extends existing models to describe regular orbits with arbitrary vertical amplitudes using a continuous, observable-based formula.

Findings

Unified formula valid for all vertical amplitudes

Connects previous thin-disk and high-amplitude models

Depends only on observable surface density

Abstract

We recently found that regular orbits in axially symmetric galactic disks have their envelopes $Z(R)$ accurately described by the relation $Z(R)\propto[\Sigma_I(R)]^{-1/3}$, if their amplitudes are comparable to the disk thickness, where $\Sigma_I$ is the surface density of the disk (integrated over its whole vertical range). Moreover, the usual adiabatic approximation gives a good description of the orbits' envelopes for low vertical amplitudes. However, these two approaches are apparently disconnected, since their expressions differ qualitatively. Our purpose in this paper is to fill this gap by extending these previous formulae to regular orbits with arbitrary vertical amplitudes inside the disk. We compare existing $Z(R)$ estimates: the razor-thin disk case, the adiabatic approximation (low-amplitude orbits in three-dimensional disks), and the integrated surface-density estimate (high-amplitude orbits in three-dimensional disks) in order to establish a connection between them. The formula presented here links the aforementioned results in an elegant and continuous way, being valid for vertical amplitudes throughout the whole vertical extension of the disk and with an expression which has the same form for all regimes. The advantage of the present formalism is the dependence of $Z(R)$ only on observable quantities, namely the disk's vertically integrated surface density, without the need to obtain the gravitational potential for the system.

Equations

$$I_{3}^{\ (AA)}=Z(R)\big{[}\Phi_{zz}(R,0)\big{]}^{1/4}$$

$$Z(R)\propto\big{[}\Phi_{zz}(R,0)\big{]}^{-1/4}$$

$$Z(R)\propto\big{[}\Sigma_{I}(R)\big{]}^{-1/3}$$

$${\Sigma }_I(R)={\int }_{-\zeta }^{\zeta }\rho (R,z)dz,$$

$${\nabla }^2\Phi =\frac{{\partial }^2\Phi }{\partial z^2}+\frac{1}{R}\frac{d{v}_c^2}{dR}=4\pi G\rho$$

$$Z(R)\propto\big{[}\Phi_{zz}(R,0)\big{]}^{-1/4}\propto\big{[}\Sigma_{0}(R)\big{]}^{-1/3}$$

$${\Sigma }_0(R)={\int }_{-Z(R)}^{Z(R)}\rho (R,z)dz.$$

$${\Sigma }_0(R)\approx 2{\rho }_0(R)Z(R).$$

$${\Sigma }_0(R)\propto Z(R){\Phi }_{zz}(R,0),$$

$${\Phi }_{zz}(R,0)\propto \frac{{\Sigma }_0(R)}{Z(R)}.$$

$$Z^4(R)\propto \frac{1}{{\Phi }_{zz}(R,0)}\propto \frac{Z(R)}{{\Sigma }_0(R)},$$

$$Z(R)\propto\big{[}\Sigma_{0}(R)\big{]}^{-1/3}.$$

$$\Sigma (R)={\int }_{-Z(R)}^{Z(R)}\rho (R,z)dz,$$

$$\frac{Z(R)}{Z(R^{\prime })}={\left[\frac{\Sigma (R^{\prime })}{\Sigma (R)}\right]}^{1/3}$$

$$Z(R)=A\big{[}\Sigma(R)\big{]}^{-1/3},$$

$$Z_o\to Z_1\to Z_2\to ...\to Z_N.$$

$$a_{n}=\operatorname*{max}_{R\in I}\,\frac{\big{|}Z_{n}(R)-Z_{n-1}(R)\big{|}}{\big{|}Z_{n}(R)\big{|}}$$

$$\Phi =-\frac{GM}{\sqrt{R^2+(a+(z^2+b^2{)}^{1/2}{)}^2}},$$

$${\Sigma }_{\Psi }(R)=f(R){\int }_{-Z}^Z\rho (R,z)dz,$$

Full text

††pubyear: 2019††pagerange: Envelopes and vertical amplitudes of disk-crossing orbits – Envelopes and vertical amplitudes of disk-crossing orbits

1 Introduction

It is widely accepted nowadays that general Hamiltonian dynamical systems have their phase space composed of both regular and chaotic regions. This fact was confirmed by a variety of numerical experiments (e.g. Henon & Heiles, 1964; Contopoulos, 2002; Hunter, 2003, 2005; Ramos-Caro et al., 2008). Seminal studies showed this pattern decades ago in galactic dynamics (Contopoulos, 1960; Ollongren, 1962; Contopoulos, 1963; Henon & Heiles, 1964, see Contopoulos 2002 for a thorough review), starting the “quest” for the third integral of motion $I_{3}$ in axially symmetric disk systems. Integrability by means of polynomial invariants was studied in Hietarinta (1987), whereas Stäckel potentials (Kuzmin, 1952, 1956; de Zeeuw, 1985, 1988; Dejonghe & de Zeeuw, 1988; de Zeeuw & Hunter, 1990) became the standard example of integrable galactic models.

The current theoretical understanding of the third-integral formalism in axially symmetric disk galaxies is based on the existence of adiabatic invariants around equatorial circular orbits (Binney & Tremaine, 2008; Binney, 2010; Binney & McMillan, 2011), as well as on the order-by-order expansion of the invariant (Contopoulos, 1960, 1963, 2002). There was also a recent attention devoted to the numerical obtention of action-angle coordinates for nearly integrable regions in phase space (McMillan & Binney, 2008; Binney, 2012a; Sanders, 2012; Sanders & Binney, 2014), as well as to efficient methods of obtaining polynomial approximate integrals (Bienaymé & Traven, 2013) and closed-formula expressions for $I_{3}$ (Bienaymé et al., 2015; Vieira & Ramos-Caro, 2016). These many complementary approaches improved the development of dynamical models for the Galaxy (de Zeeuw, 1987; Binney, 2010; Binney & McMillan, 2011; Binney, 2012b), which are a key ingredient to analyze and interpret the new kinematic data obtained from Galaxy surveys (Binney & Sanders, 2014), such as Gaia DR2 (Gaia Collaboration et al., 2018).

Let us assume an axially symmetric system. An approximate third integral gives us, in particular, the form of the orbit’s amplitude $Z$ as a function of radius $R$ on the $R$–$z$ plane (its “envelope" $Z(R)$). In spite of the significant advances in recent years, closed-form analytical approximations for the non-classical integral $I_{3}$ (and for the corresponding envelope) are still very scarce, and the existent results have limited range of validity. The usual adiabatic approximation based on decoupled radial and vertical motions (Binney & Tremaine, 2008), for instance, is only valid very close to the equatorial plane (Binney & McMillan, 2011; Vieira & Ramos-Caro, 2014). It was recently found, based on the study of orbits in razor-thin disks (Vieira & Ramos-Caro, 2015, 2016), that a simple relation can describe the envelopes $Z(R)$ of disk-crossing orbits whose amplitudes are near the disk’s vertical edge (Vieira & Ramos-Caro, 2014). The orbit’s amplitude scales with a fixed power of the integrated surface density. However, the formula fails for intermediate and low vertical excursions. In summary, the two above limits, with near-equatorial and near-edge vertical amplitudes, can be understood in terms of straightforward physical arguments. But, in order to have a quantitative physical description of the behaviour of off-equatorial orbits with varying vertical amplitudes, a connection between these two regimes must be done.

The purpose of this paper is to fill this gap. We present an extension of the previous formulae, which connects the two regimes described above in an unified and elegant way. We obtain an expression which describes the envelopes $Z(R)$ of disk-crossing (tube) orbits with vertical amplitudes spanning the whole vertical extension of the disk, and which gives the aforementioned known approximations in their respective ranges of validity. Moreover, this expression for $Z(R)$ depends only on the disk’s dynamical surface density, integrated up to the orbit’s vertical amplitude (see Eqs. (13)–(14)). In this way, the density profile completely determines the envelopes of tube orbits for any vertical range of amplitudes inside the disk.

The paper is structured as follows: Section 2 summarizes the two known cases which give good predictions for the envelopes $Z(R)$: the adiabatic approximation (2) (orbits with very small vertical amplitudes) and the surface-density estimate (3) for the envelopes of orbits near the disk edge. It also presents the unified framework which connects the adiabatic approximation and the near-edge estimate via an elegant procedure, based solely on the disk’s surface density (integrated up to the orbit’s amplitude). Section 3 briefly comments on the method of successive approximations used to numerically calculate the envelopes given by the integral equation resulting from our new framework, and Section 4 presents numerical results which confirm our predictions. We present our conclusions in Section 5.

2 $Z(R)$ envelopes on the $R$–$z$ plane

We work in cylindrical coordinates ($R,z,\varphi$) and consider axially symmetric disks with a density profile $\rho(R,z)$ concentrated around the equatorial plane $z=0$. As mentioned in the Introduction, there are different approaches to estimate the orbits’ envelopes on the meridional plane $R$–$z$. We know that the epicyclic approximation gives us an adiabatic invariant of the form (when evaluated at the envelope $Z(R)$)

[TABLE]

as approximate third integral for the system (Binney & Tremaine, 2008), where $\Phi$ is the gravitational potential, $\Phi_{zz}=\partial^{2}\Phi/\partial z^{2}$, and $Z$ is the orbit’s vertical amplitude. This expression results in (Binney & Tremaine, 2008)

[TABLE]

for the vertical amplitudes of nearly equatorial orbits. However, it is well known that this approximation has a very limited range of applicability, being valid only for orbits whose vertical excursions are very small (Binney & McMillan, 2011). In particular, their vertical amplitudes must be much smaller than the disk thickness (Vieira & Ramos-Caro, 2014). Variants of this formula were proposed in order to overcome this problem, as for instance the corrections to the adiabatic approximation of Binney & McMillan (2011). In spite of the fact that this correction has a wider range of validity and can describe accurately the corresponding off-equatorial orbits, it introduces additional phenomenological parameters which must be adjusted a posteriori (Binney & McMillan, 2011, see also Vieira & Ramos-Caro 2014).

On the other hand, Vieira & Ramos-Caro (2014) obtained a simple expression for the envelopes of disk-crossing orbits whose vertical amplitudes are comparable to the disk thickness $\zeta$:

[TABLE]

where

[TABLE]

and $\rho$ is the system’s density distribution (which may contain different disk components or even a sparse dark-matter halo). This formula was based on our results for razor-thin disks (Vieira & Ramos-Caro, 2015, 2016) and depends only on the integrated surface density $\Sigma_{I}$, which we call for brevity “disk’s surface density”. Therefore, Eq. (3) can be evaluated from observational data alone (without the knowledge of the gravitational potential). However, although we know that expressions (2) and (3) are both valid in their respective ranges of applicability, there is a variety of disk-crossing orbits which “live” in the intermediate regions, that is, whose vertical amplitudes are spread throughout the range $0<Z<\zeta$. These orbits cannot be described by neither of the expressions (2) or (3), since the assumptions which led to each of these equations cease to be valid for intermediate-amplitude orbits.

2.1 Connecting the formulae

In order to keep the framework of writing the orbits’ envelopes for arbitrary vertical amplitudes only in terms of the disk density, we must find a link between expressions (2) and (3). Doing this in an unified way will allow us to obtain an expression which is valid for the whole range of vertical amplitudes inside the disk, and which has the correct limits for almost-zero and near-edge amplitudes. The connection is done once we note that Poisson equation can be written in the form (Binney & Tremaine, 2008)

[TABLE]

on the equatorial plane, where $v_{c}$ is the speed of a circular orbit of radius $R$. Since most disk galaxies have very small velocity gradients far from the central region (Sofue & Rubin, 2001), the term involving $v_{c}$ may be neglected if the density distribution is large enough on the equatorial plane. We would like to have

[TABLE]

in the domain of validity of the adiabatic approximation (2) (i.e. for very small vertical amplitudes), where $\Sigma_{0}$ is an appropriately defined integrated density. The most natural choice would be

[TABLE]

Since $Z(R)$ is small, we may consider the integrand in the above expression as a constant function of $z$ (evaluated at $z=0$), $\rho(R,z)\approx\rho_{0}(R)$, which gives us

[TABLE]

With the aid of Eq. (5), and neglecting the term involving the velocity gradient, we obtain

[TABLE]

and thus

[TABLE]

Substituting in Eq. (2), we obtain

[TABLE]

that is,

[TABLE]

Therefore, Eqs. (2) and (3) will have the same form if the “surface density” appearing in the formula (3) for $Z(R)$ is, instead of (4),

[TABLE]

where the limit of integration $Z(R)$ is the orbit’s vertical amplitude (which depends on the galactocentric radius $R$). Note that, when the orbit’s amplitude is comparable to the disk thickness, the two formulae (3) and (14) give practically the same result. We call the $\Sigma(R)$ given by Eq. (14) the “orbit’s surface density”. The orbit’s amplitude as a function of $R$ will be given by the implicit equation

[TABLE]

for $R$, $R^{\prime}$ in the orbit’s radial range, with $\Sigma$ given as a function of $Z$ by expression (13). We claim that this expression is valid through the whole vertical range of the disk, that is, for regular orbits with amplitudes in the whole range $0<Z<\zeta$.

We might approximate the vertical amplitude $Z(R)$ appearing in (13) by an average value along the orbit, if the amplitudes do not vary much. In this case, Eq. (14) is given explicitly as a function of $R$, once we have an estimate for the vertical amplitude, and the equations would give an explicit formula for $Z(R)$. However, this procedure can produce additional errors in the prediction for the envelope $Z(R)$, since we neglect the $R$-dependence of the integration interval: the ratio $\Sigma/\Sigma_{I}$ is then more likely to keep a constant value along the orbit. The prediction for $Z(R)$ will therefore be very similar to Eq. (3), and the correspondence would be exact if the density profile had a separable form, $\rho(R,z)=\Sigma_{I}(R)\xi(z)$, with $\int\xi dz=1$. For flattened disks, the aforementioned plane-parallel case is only approximately valid in general, so an important piece of information is being underestimated if we neglect the $R$-dependence of the integration limits. Moreover, Eq. (3) is only valid for amplitudes which are close to the disk’s vertical edge (Vieira & Ramos-Caro, 2014), so in the discussion which follows we will always refer to the implicit expression (14), where the integration limits in (13) are the actual values of the orbit’s envelope. We test below this expression for a variety of orbits in disk-like potentials, with amplitudes varying in the range $0<Z\lesssim\zeta$.

3 An iterative procedure for $Z(R)$

Equation (14) does not give us a direct method to obtain $Z(R)$, since the integration limits of $\Sigma$, Eq. (13), are precisely the function sought. We would have then to solve an integral equation to obtain the orbit’s envelope. One way to circumvent this issue is to consider an iterative procedure by means of a Neumann series (the well-known “method of successive approximations”, see Arfken & Weber, 1972), in which a function $Z_{o}(R)$ is chosen as an initial guess, inserted on the right-hand side of the equation for $Z(R)$ (see Eqs. (13)–(14), with $\pm Z_{o}(R)$ being the integration limits for $\Sigma$):

[TABLE]

where $A$ is a constant and $\Sigma(R)$ is calculated from Eq. (13). This first step will generate a new function $Z_{1}(R)$, which will be used in the second step as the new integration limit for $\Sigma$ in Eq. (15), giving rise to a function $Z_{2}(R)$. After $N$ steps, we will have a sequence of functions

[TABLE]

Ideally, the sequence of functions $(Z_{n})$ will be uniformly convergent in the range of radii relevant for the orbit, and its limit will be the solution $Z(R)$ to Eqs. (13)–(14). In that case, the sequence ($a_{n}$) given by

[TABLE]

will converge to zero, where $I$ is a previously defined interval corresponding to the orbit’s radial range. We therefore look for an approximate solution to these equations by obtaining a value of $N$ for which $a_{N}<\epsilon$, where $\epsilon$ is a tolerance parameter. The iterated function $Z_{N}$ will therefore be an approximation to the actual solution $Z(R)$ within the error $\epsilon$. This procedure is, therefore, dependent only on the initial guess $Z_{o}(R)$ for $Z$; it does not rely upon an a priori knowledge of the envelope. In this way, it can be used to make predictions about the orbit’s shape knowing (for instance) only one point on the zero-velocity curve. We will follow this procedure to numerically calculate the orbits’ envelopes for specific example disk potentials.

4 Tests with simple disk potentials

We perform in this Section numerical experiments in specific disk-like potentials in order to test the validity of prediction (13)–(14), for a variety of orbits with amplitudes spread in the range $0<Z(R)<\zeta$. We find in all cases that expression (14), with $\Sigma(R)$ given by the orbit’s surface density (13), is a good approximation to the real orbits’ envelopes, with errors below 1%. The prediction (14) is better than both the adiabatic approximation (2) and the estimate based on the disk’s surface density (3), even in their ranges of validity (very low vertical amplitudes and amplitudes near the disk edge, respectively).

We choose as reference the point on the orbit with the largest vertical amplitude $Z$. The predictions for the envelopes given by the different approximations are made to agree with the vertical amplitude of the numerically integrated orbit at this point, and we compare the three formulae (2), (3), and (14) along the orbit’s radial range. Equations (2) and (3) give explicitly $Z$ as a function of $R$ and can be readily compared with the numerically calculated envelope. On the other hand, Eq. (14), with $\Sigma(R)$ given by the orbit’s surface density (13), does not give $Z(R)$ explicitly. In order to obtain the predicted envelope we apply the iterative method described in Section 3, starting with $Z_{o}(R)$ given by the adiabatic approximation (2), Z_{o}(R)\propto\big{[}\Phi_{zz}(R,0)\big{]}^{-1/4}. We perform 10 iterations ($N=10$) for each orbit in order to obtain an accurate approximation to Eqs. (13)–(14). Very good approximations already appear for smaller $N$.

4.1 Miyamoto-Nagai disk

We consider as the main example the Miyamoto-Nagai potential (Miyamoto & Nagai, 1975; Binney & Tremaine, 2008)

[TABLE]

which is commonly used to model the disk component of Galaxy mass models (Allen & Santillán, 1991; Irrgang et al., 2013; Barros et al., 2016) and also to analyze the Saggitarius dwarf galaxy debris (Johnston et al., 1995; Helmi, 2004; Law et al., 2009; Ibata et al., 2013; Deg & Widrow, 2013). The disk’s vertical edge is estimated as $\zeta/a\approx 3b/a$ if $b\ll a$ (Vieira & Ramos-Caro, 2014). We analyze various orbits for different values of energy $E$, angular momentum $L_{z}$, and initial conditions. We consider only the cases with a sufficiently flattened density profile ($b/a\lesssim 1/10$).

We see in Figs. (1)–(3) representative examples of orbits with large radial extent and different vertical amplitudes (when compared with the disk thickness). The orbits’ envelopes are accurately described by Eqs. (13)–(14). If the vertical amplitudes are very small, being in the region where the usual adiabatic approximation (2) is valid (as in Fig. 1), both expressions (2) and (14) give good approximations for the envelopes, as expected from the above considerations. Moreover, as it can be seen on the right panel of Fig. 1, our proposed envelope (14) gives an improvement to the adiabatic approximation, since it takes into account the whole integration interval for each galactocentric radius, and not only an equatorial-plane quantity (showing differences even when the range of integration is very small). In other words, the dependence of $\rho$ on $z$ in the integrand of (13) is important in the description of nearly equatorial orbits, since it gives better predictions than the equatorial-plane approximation (8).

This effect is also seen in orbits with intermediate amplitudes lying inside the disk ($Z(R)<\zeta$, see Fig. 2). The prediction from the orbit’s surface density, Eqs. (13)–(14), gives results which are significantly better than either the adiabatic approximation (2) or the disk’s surface density estimate (3), as illustrated on the right panel of Fig. 2. For amplitudes near the disk edge $\zeta=3b$, Eqs. (13)–(14) give predictions close to the disk’s surface density estimate (3), but with a small improvement, as exemplified in Fig. 3. As $Z$ approaches the disk edge, the difference between predictions (3) and (14) becomes negligible. The prediction from (14) always lies between the adiabatic approximation (2) and the disk’s surface-density approximation (3) for the Miyamoto-Nagai disk and it is, in general, a better description than both of these formulae.

The errors in the orbits’s amplitudes, calculated from Eq. (14) with the procedure described in Section 3, are less than 1% and are systematically smaller than the corresponding errors of formulae (2) and (3) (see Figs. 1–3). These errors are of the same order of magnitude as the errors appearing in the comparison of formula (3) with flattened Stäckel models (Vieira & Ramos-Caro, 2014). Indeed, numerical experiments with Stäckel potentials (the Kuzmin-Kutuzov potential of Dejonghe & de Zeeuw, 1988) show that the predictions of Eqs. (13)–(14) are also the best approximation for the orbits’ envelopes.

Tests including a sparse spherical halo do not influence the predictions, as expected from the numerical experiments of Vieira & Ramos-Caro (2014) concerning near-edge amplitudes (although in this case we must consider in Eq. (13) the total dynamical density of the system “disk+halo”). The present expression is therefore an improvement for the envelopes’ estimates of regular orbits inside the disk whenever the disk component is dominant, i.e. when the orbit is far from the central (spheroidal) bulge component.

5 Conclusions and final remarks

We presented in this paper an extension of the modeling of the envelopes of disk-crossing orbits on the system’s meridional plane. It is based on the disk’s surface density distribution, which was previously presented for the cases of off-equatorial orbits around razor-thin disks (Vieira & Ramos-Caro, 2015, 2016) and of orbits in three-dimensional disks with vertical amplitudes comparable to the disk thickness (Vieira & Ramos-Caro, 2014). The generalization discussed here encompasses the two known limits for the envelope estimates: The limit of very low vertical amplitudes (usual adiabatic approximation, Eq. (2), see for example Binney & Tremaine, 2008) and the case of vertical amplitudes near the disk thickness (the integrated surface density approach, Eqs. (3) and (4), see also Vieira & Ramos-Caro, 2014).

The proposed formula, Eq. (14) with $\Sigma$ given by Eq. (13), is valid for regular tube orbits with arbitrary vertical amplitudes $0<Z<\zeta$ (where $\zeta$ is the disk thickness), i.e. for arbitrary vertical amplitudes inside the disk. All numerical experiments gave us relative errors smaller than 1% for the predictions for $Z(R)$, which were systematically smaller than the predictions from the adiabatic approximation (2) and the predictions via the disk’s surface density (3). In view of the present formalism and of previous results obtained for near–$\zeta$ vertical amplitudes (Vieira & Ramos-Caro, 2014), this formula is expected to be valid for all regions in which the disk component is dominant.

It is worth noting that our formalism deals only with the envelope $Z(R)$ of the orbit, that is, with the “$z$-limits” of the orbit, since our approach is based on an extension of the vertical action variable close to the equatorial plane. The inner and outer “$R$-limits” of the orbit (the “$R(z)$-envelopes”) could be obtained, in principle, by a procedure which generalizes the results obtained via the radial action variable on the equatorial plane. Such procedure is beyond the scope of this work.

A direct extension of the above formalism allows us to analyze the case of very flattened, three-dimensional axially symmetric disk systems in modified theories of gravity which admit a Hamiltonian flow for test particles, $H=p^{2}/2+\Psi$, where $\Psi$ is a modified gravitational potential. This case was presented in Vieira & Ramos-Caro (2014) for orbits whose vertical amplitude is comparable to the disk thickness. In an analogous manner, we obtain for arbitrary vertical amplitudes that formula (14) is valid with $\Sigma$ given by

[TABLE]

where $Z=Z(R)$ is the orbit’s vertical amplitude at radius $R$ and $f(R)$ is the radial dependence of the modification to the potential, given by \partial\Psi/\partial|z|\big{|}_{z=0}=f\,\Sigma when the disk component of the system is modelled as razor-thin. There are two example cases which fit into this category of theories: MOND (Bekenstein & Milgrom, 1984) gives us $\Sigma_{\Psi}=\Sigma/\mu$, with $\Sigma$ given by Eq. (13) and $\mu$ is MOND’s interpolating function, and RGGR (Rodrigues et al., 2010) gives us $\Sigma_{\Psi}=[1-V_{\infty}^{2}/\Phi_{N}(R,0)]\,\Sigma$, where $\Sigma$ is given by Eq. (13), $V_{\infty}$ is the asymptotic circular velocity of test particles on the equatorial plane, and $\Phi_{N}$ is the corresponding Newtonian potential. More details about this formalism and about the case of near-edge vertical amplitudes can be found at Vieira & Ramos-Caro (2014).

The problem of finding the approximate third integral of motion (in terms of phase-space coordinates) which gives the envelopes (14) is a much more difficult task. Although the corresponding formula was recently found for orbits crossing razor-thin disks (Vieira & Ramos-Caro, 2016), this integral for three-dimensional disks would have to involve the orbit’s surface density (13), which is a function of $R$ and $z$ (and not only on $R$, as it is in the razor-thin case). Moreover, it is not clear how the kinetic and the surface-density dependent terms would behave in the intermediate regions between small and high $z$. This question is still open and deserves a further study.

Acknowledgments

We thank Alberto Saa for fruitful discussions. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. RSSV acknowledges the support of CAPES and of São Paulo Research Foundation (FAPESP), grant 2010/00487-9. The authors acknowledge the support of FAPESP grant 2013/09357-9.