© ESO, 2015

1. Introduction

Tracer objects orbiting the Milky Way can be used to infer the gravitational field at large Galacto-centric radii. Jeans modelling links the field with the available kinematical data, under the important assumption that these objects can be described by a collision-less system of test bodies in steady-state equilibrium (Jeans 1915). Integral to this approach is a phase-space distribution function (PDF), while the physical observables (e.g. the number density, the velocity dispersion ellipsoid, etc.) are secondary quantities that are PDF-dependent functionals on the phase space.

Galaxy mass can be estimated based on the radial velocity dispersion (RVD) data. In the literature, one can find mass values of $\mathrm{4.}{\mathrm{2}}_{-0.4}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$ (Deason et al. 2012a), $\mathrm{4.}{\mathrm{9}}_{-1.1}^{\mathrm{+}\mathrm{1.1}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$ (Kochanek 1996), $\mathrm{5.}{\mathrm{4}}_{-3.6}^{\mathrm{+}\mathrm{0.2}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$ (Wilkinson & Evans 1999), and $\mathrm{5.}{\mathrm{4}}_{-0.4}^{\mathrm{+}\mathrm{0.1}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$ (Sakamoto et al. 2003), all within 50 kpc; $\mathrm{4.}{\mathrm{0}}_{-0.7}^{\mathrm{+}\mathrm{0.7}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$, within 60 kpc (Xue et al. 2008); or (5.8−6.0) × 10¹¹ M_⊙ enclosed within 100 kpc (Klypin et al. 2002). On the other hand, depending on the model assumptions, the virial mass is (8−10−12) × 10¹¹ M_⊙ (Battaglia et al. 2005; Xue et al. 2008; Kafle et al. 2014) or even (18−25) × 10¹¹ (Sakamoto et al. 2003). The kinematics of an extended orphan stream indicates the mass of ~2.7 × 10¹¹ M_⊙ within 60 kpc (with disk+bulge mass of 1.3 × 10¹¹ M_⊙; Newberg et al. 2010; Sesar et al. 2013), significantly less than suggested by the above estimates within 50 kpc. This points to some model-dependent effects.

In this context, it is natural to ask what the lower bound is for the Galaxy mass indicated by the kinematics of the outermost tracers, paying more attention to the phase-space model rather than the particular mass model.

With the simplest working hypothesis of the absence of the extended halo, the Galactic gravitational field at large distances, to a fair degree of approximation, would be that of a point mass. In this approximation, a single total mass parameter is decisive for the field asymptotics of any compact mass distribution. Bahcall & Tremaine (1981) proposed in this approximation for the neighbouring galaxies, a mass estimator $\frac{\mathit{C}}{\mathit{G}}{}^{\mathrm{⟨}}\mathit{v}_{\mathit{z}}^{\mathrm{2}}\hspace{0.17em}{\mathit{R}}^{\mathrm{⟩}}$, with the averaging performed over distant tracers at various projected radii R, where C is a constant. More recently, Watkins et al. (2010) considered a spherical symmetric counterpart of this estimator, $\frac{\mathit{C}}{\mathit{G}}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\hspace{0.17em}{\mathit{r}}^{\mathit{\gamma }}\mathrm{⟩}$, with an arbitrary power of the radial distance r.

Mass estimators based on Jeans theory, irrespective of the adopted mass model, are related to a PDF restricted, in particular, by some indirect constraints appearing because of the assumptions made about the secondary quantities. The restrictions usually concern the flattening of the velocity dispersion ellipsoid β. This quantity is poorly known for peripheral tracers. Introducing a variable β leads to difficulties in solving the Jeans problem. To overcome this, β is often assumed to be a position independent parameter. On the other hand, this is too much constraining an assumption, since any limitation on β indirectly imposes restrictions on the function space admissible for PDF’s. We conjectured (Bratek et al. 2014) that the lower bound for the total mass may increase in response to these constraints, while there is no definite upper bound. Consequently, the mass is likely to be overestimated.

Recently, there is a growing interest in methods of determining the general assumption-free PDF’s from the kinematical data. Magorrian (2014) proposed a framework in which the gravitational potential is inferred from a discrete realization of the unknown distribution function using snapshots of stellar kinematics. Our previous article (Bratek et al. 2014) is placed within this field of interest. Therein, we proposed a method of determining PDF’s from a given spherically symmetric RVD profile, without imposing any constraints on the secondary quantities. Even in the simplest case of a point mass approximation, our method allowed us to faithfully reconstruct the shape of the RVD profile, including its low-size and variable features. By considering various PDF’s giving rise to RVD’s overlapping with that observed at larger radii, we showed that there is no upper bound for the total mass, while there is a sharp lower bound, which is slightly below 2.0 × 10¹¹ M_⊙. For lower masses, no PDF could be found to account for the measured RVD within the acceptable limits.

The lower bound may also suggest a Galaxy mass lower than given in the literature. A natural question arises as to whether a low mass is physically reasonable and could appear in other models. This cannot be excluded. A mass (2.4−2.6) × 10¹¹ would be consistent with past results in a three-component mass model (Merrifield 1992), with the best estimate in the point mass field (Little & Tremaine 1987), and most remarkably, with the recent value inferred from the kinematics of the Orphan stream (Newberg et al. 2010; Sesar et al. 2013). The aim of the present paper is to verify the reliability of the low-mass solution, when the point mass potential is replaced by a more realistic potential.

When considering a three-component model in place of a point mass potential, one has to take the Galactic rotation curve into account as a constraint on the mass distribution profile. Based on this, McMillan (2011) found the density distribution for various Galactic components. In particular, the author considered a model consisting of a stellar disk of mass 6.43 × 10¹⁰ M_⊙ and a Navarro-Frenk-White (NFW) dark matter halo (Navarro et al. 1997). The mass function of the corresponding density distribution is such that the mass enclosed within a sphere of radius r = 20 kpc is M₍₂₀₎ = 2.58 × 10¹¹ M_⊙, while that at r = 16.8 kpc is M_(16.8) = 2.2 × 10¹¹ M_⊙. A mass function of this order of magnitude at r = 20 kpc is not only the property of the NFW profile, and should be expected for any spherically symmetric dark halo that dominates the dynamics of the disk. This is because the Keplerian mass function M(r) = rv²(r) /G evaluated at r = 20 kpc gives a value of 2.2 × 10¹¹ M_⊙ for the velocity v = 220 km s^-1. The NFW profile (as well as many other dark halo profiles considered in the literature) is non-integrable, and its mass function is divergent in the limit r → ∞. Then, the requirement of a fixed finite mass determines the size of the halo. It follows for a model consisting of a low-mass stellar disk, and a dark halo with a total mass of about the lower bound that we have found, that the matter distribution should be enclosed entirely within the sphere of radius ≈20 kpc. In this case, the potential would quickly become Keplerian at r> 20 kpc and the conclusions drawn within the point mass potential would remain valid.

A model in which the dark matter halo does not dominate the dynamics of the inner Galaxy is even more interesting because, in this situation, a slightly more massive disk contribution significantly breaks the spherical symmetry of the potential at low radii. Then the question arises as to what extent the higher multipoles of the potential could affect the motion of the halo tracer objects. We deal with this issue in the subsequent sections.

The lower bound referred to above coincides with the sum of the dynamical mass ≈1.5 × 10¹¹ M_⊙ inferred from the rotation curve in disk model (Jałocha et al. 2014; Sikora et al. 2012) and the mass (1.2−6.1) × 10¹⁰ M_⊙ of the hot gaseous halo surrounding Milky Way (Gupta et al. 2012). The gravitational potential of these components can be interpreted as a perturbation of a point mass potential, which breaks the spherical symmetry at low radii. For complicated potentials, however, the distribution integral on the phase space cannot be explicitly constructed because the first integrals characterizing admissible orbits are not known in an explicit form. To overcome this difficulty, a numerical simulation can be performed.

In this paper, we present an example of this kind of simulation in which the test bodies represent the tracer objects orbiting the Galaxy. The initial conditions for the simulation in the perturbed field should be chosen close to a known stationary solution of a Jeans problem in the non-perturbed field. For this purpose we make use of a PDF to be found similar to that found in Bratek et al. (2014), but it is not obvious if this initial PDF and the resulting RVD profile would be stable against such a perturbation. Running a simulation with an initial PDF consistent with the observed high RVD could lead to an RVD profile with a value that is quickly and steadily decreasing. If this happened, this would mean that the initial approximation was far from a stationary solution in the new potential. The main goal of the present paper is to investigate this stability issue and thus also the reliability of the point mass approximation; that is, we test whether predictions for the RVD and the mass in the new potential are comparable to those made in the point mass approximation, that is, we test the structural stability of these predictions.

The structure of this paper is the following. In Sect. 2 we recall the main ideas behind the Keplerian ensemble method of obtaining a PDF. In Sect. 3, we use a representative of possible PDF’s to set the initial conditions for our n-body simulation in the modified potential. Next, we test the stability of the resulting RVD profile, and then our conclusions follow.

2. First approximation of PDF from the Keplerian ensemble method

Mathematical preliminaries¹. An elliptic orbit of a test body bound in the field of a point mass M is fully characterized by 5 integrals of motion: the Euler angles (Φ,Θ,Ψ), determining the orbit orientation; the eccentricity e, describing the orbit flattening; and the dimensionless energy parameter $\mathit{ϵ}\mathrm{=}\mathrm{-}\frac{\mathit{RE}}{\mathit{GM}}$, describing the size of the large semi-axis (R is an arbitrary unit of length while E is the energy per unit mass). We call a spherically symmetric collection of confocal ellipses a Keplerian ensemble. From the Jeans (1915) theorem, for this ensemble in a steady state equilibrium it suffices to consider PDF’s as being functions of e and ϵ only. Accordingly, instead of r, θ, φ, v_r, v_θ, v_φ, we use phase coordinates u, θ, φ, e, ϵ, ψ defined by $\begin{array}{ccc}& & \mathit{r}\mathrm{=}\mathit{R}\hspace{0.17em}\mathit{u,} {\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{=}\frac{\mathit{GM}}{\mathit{R}}\mathrm{(}\frac{\mathrm{2}}{\mathit{u}}\mathrm{-}\frac{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathit{ϵ}\hspace{0.17em}{\mathit{u}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}\hspace{0.17em}{\mathit{ϵ}}^{\mathrm{)}}\mathit{,}\\ & & \mathrm{\left(}{\mathit{v}}_{\mathit{\theta }}\mathit{,}{\mathit{v}}_{\mathit{\phi }}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{u}}\hspace{0.17em}\sqrt{\frac{\mathit{GM}}{\mathit{R}}}{\mathrm{(}{\frac{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathit{ϵ}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{\left(}\sin \mathit{\psi ,}\cos \mathit{\psi }\mathrm{\right)}\mathit{.}\end{array}$(1)For physical reasons, we assume all orbits to be confined entirely within a region u ∈ (u_a,u_b) bounded by two spheres of radii u_a and u_b. With this approach, all orbits, with pericentra that are too low (i.e. violating the point mass approximation) or with apocentra that are too high (e.g. beyond Local Group members), are excluded. Consequently, the space of parameters (e,ϵ) gets restricted to a domain S: $\frac{\mathrm{1}\mathrm{+}\mathit{e}}{\mathrm{2}{\mathit{u}}_{\mathit{b}}}\mathit{<}\mathit{ϵ}\mathit{<}\frac{\mathrm{1}\mathrm{-}\mathit{e}}{\mathrm{2}{\mathit{u}}_{\mathit{a}}}$ and $\mathrm{0}\mathrm{\le }\mathit{e}\mathit{<}\frac{{\mathit{u}}_{\mathit{b}}\mathrm{-}{\mathit{u}}_{\mathit{a}}}{{\mathit{u}}_{\mathit{b}}\mathrm{+}{\mathit{u}}_{\mathit{a}}}\mathit{<}\mathrm{1}$. On integrating out the angles θ,φ,ψ, the principal integral ^∫f(r,v) d³r d³v reduces (to within an unimportant constant factor) to $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\underset{{\mathit{u}}_{\mathit{a}}}{\overset{{\mathit{u}}_{\mathit{b}}}{\mathrm{\int }}}{\mathit{\nu }}_{\mathit{u}}\mathrm{\left[}\mathit{f}\mathrm{\right]}\hspace{0.17em}\mathrm{d}\mathit{u,} {\mathit{\nu }}_{\mathit{u}}\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{\equiv }\underset{}{{\mathrm{\int }}^{\mathit{S}\mathrm{\left(}\mathit{u}\mathrm{\right)}}}\frac{\mathit{e}\hspace{0.17em}\mathrm{d}\mathit{e}\hspace{0.17em}\mathrm{d}\mathit{ϵ}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{e,ϵ}\mathrm{\right)}}{\sqrt{\mathit{ϵ}\hspace{0.17em}\mathrm{(}\mathit{ϵ}\mathrm{-}{\frac{\mathrm{1}\mathrm{-}\mathit{e}}{\mathrm{2}\mathit{u}}}^{\mathrm{)}}\mathrm{(}\frac{\mathrm{1}\mathrm{+}\mathit{e}}{\mathrm{2}\mathit{u}}\mathrm{-}{\mathit{ϵ}}^{\mathrm{)}}}}\mathrm{·}\end{array}$(2)The integration domain S(u) ⊂ S is a u-dependent quadrilateral region $\mathit{ϵ}\mathrm{\in }{}^{\mathrm{\right(}}\mathrm{Max}{}^{\mathrm{\right(}}\frac{\mathrm{1}\mathrm{-}\mathit{e}}{\mathrm{2}\mathit{u}}\mathit{,}\frac{\mathrm{1}\mathrm{+}\mathit{e}}{\mathrm{2}{\mathit{u}}_{\mathit{b}}}^{\mathrm{\left)}}\mathit{,}\mathrm{Min}{}^{\mathrm{\right(}}\frac{\mathrm{1}\mathrm{+}\mathit{e}}{\mathrm{2}\mathit{u}}\mathit{,}\frac{\mathrm{1}\mathrm{-}\mathit{e}}{\mathrm{2}{\mathit{u}}_{\mathit{a}}}^{\mathrm{\left)}}^{\mathrm{\left)}}$,$\mathit{e}\mathrm{\in }{}^{\mathrm{\right(}}\mathrm{0}\mathit{,}\frac{{\mathit{u}}_{\mathit{b}}\mathrm{-}{\mathit{u}}_{\mathit{a}}}{{\mathit{u}}_{\mathit{a}}\mathrm{+}{\mathit{u}}_{\mathit{b}}}^{\mathrm{\left)}}$, each point of which corresponds to a spherically symmetric pencil of confocal elliptical orbits intersecting a sphere of a certain radius u. The functional ν_u [ f ] has the interpretation of the probability density for the variable r/R to fall within the spherical shell u<r/R<u + du.

Given a PDF f(e,ϵ), the expectation value ⟨ g ⟩ _r for an observable g = g(e,ϵ,u) inside that shell equals $\mathrm{⟨}\mathit{g}{\mathrm{⟩}}_{\mathit{r}}\mathrm{=}\frac{{\mathit{\nu }}_{\mathit{u}}\mathrm{\left[}\mathit{fg}\mathrm{\right]}}{{\mathit{\nu }}_{\mathit{u}}\mathrm{\left[}\mathit{f}\mathrm{\right]}}\mathrm{·}$(3)In particular, given M and f(e,ϵ,u), the model RVD profile $\mathrm{⟨}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathit{/}\mathit{G}$ is obtained with $\mathit{g}\mathrm{\left(}\mathit{e,ϵ,u}\mathrm{\right)}\mathrm{\equiv }\frac{\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}{\mathit{G}\hspace{0.17em}\mathit{M}}\mathrm{=}\mathrm{2}\mathrm{-}\frac{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}\mathit{ϵ}\hspace{0.17em}\mathit{u}}\mathrm{-}\mathrm{2}\hspace{0.17em}\mathit{ϵ}\hspace{0.17em}\mathit{u}\hspace{0.17em}$ substituted for g in Eq. (3)².

But we are concerned with the inverse problem: given an RVD profile matching the observations, we want to derive a distribution function f the RVD profile would follow from. This problem can be solved as follows. First, we consider an auxiliary function h(e,ϵ) such that f = h² (then f ≥ 0, as required for a probability density) and make a series expansion in polynomials orthogonal on S, i.e. $\mathit{h}\mathrm{\left(}\mathit{e,ϵ}\mathrm{\right)}\mathrm{\approx }{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{D}}_{\mathit{d}}}{\mathit{h}}_{\mathit{k}}\hspace{0.17em}{\mathrm{𝒬}}_{\mathit{k}}\mathrm{\left(}\mathit{e,ϵ}\mathrm{\right)}\mathit{.}$(4)The ’s are constructed with the help of a Gramm-Schmidt orthogonalization method on S. Next, given a mass parameter M, we find an optimum sequence of expansion coefficients h_k by minimizing a discrepancy measure between a) the RVD profile from measurements, ${\overline{\mathit{p}}}_{\mathit{r}}\mathrm{\equiv }\overline{\mathit{r}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}\mathit{/}\mathit{G}$, where the averaging is taken over all halo tracer objects within a spherical shell of some width and a given radius r; and b) the model ⟨ g ⟩ _r profile calculated from Eq. (3) with the help of the function h corresponding to the optimum h_k’s. With these h_k’s, the discrepancy measure can be reduced further by replacing M with a better fit value, e.g. $\mathit{M}\mathrm{\to }{\mathit{M}}_{\mathit{bf}}\mathrm{=}\frac{{\sum }_{\mathit{r}}\mathrm{⟨}\mathit{g}{\mathrm{⟩}}_{\mathit{r}}{\overline{\mathit{p}}}_{\mathit{r}}}{{\sum }_{\mathit{r}}\mathrm{⟨}\mathit{g}{\mathrm{⟩}}_{\mathit{r}}^{\mathrm{2}}}$ if the ${\sum }_{\mathit{r}}{\left(\mathit{M}\mathrm{⟨}\mathit{g}{\mathrm{⟩}}_{\mathit{r}}\mathrm{-}{\overline{\mathit{p}}}_{\mathit{r}}\right)}^{\mathrm{2}}$ norm is used.

This way, a PDF f(e,ϵ) consistent with the RVD measurements can be reconstructed, provided M is large enough. For M above a limiting value M_cut, there is always a PDF for which the RVD profile is perfectly accounted for, while below this limit no satisfactory fit can be found. For M>M_cut, increasing the number D_d of the basis polynomials , efficiently decreases the fit residuals, but for D_d high enough the residuals appear to tend to some small non-zero limit. For M<M_cut the fit residuals remain very large, regardless of D_d, and rapidly increase with decreasing M. This shows that M_cut is the lower bound for the mass in the point mass approximation.

2.1. Profile of RVD from measurements

Without transverse velocity components, the radial motions of kinematic tracers cannot be unambiguously transformed from the local standard of rest frame (LSR) to the Galacto-centric frame. However, for a spherically symmetric distribution of tracers one can try to assume a β(r) profile or find a self-consistent one by iterations. The particular model of β(r) affects the RVD significantly only inside a spherical region of several r_⊙ in diameter. We must bear in mind, however, a twofold influence of the particular model of β on the total mass determination: both β(r) itself and the so obtained β-dependent RVD profile enter the spherical Jeans equations.

2.1.1. Formula relating the LSR radial motion measurements to the Galacto-centric RVD

Here, we consider a spherically symmetric ensemble of test bodies described by some PDF and the resulting β(r), and then also $\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ (averaging over spherical shells). For a test body with a velocity vector v in the Galacto-centric coordinate frame, the radial and tangential components of v are v_r = v°e_r, v_θ = v°e_θ, v_φ = v°e_φ, with ${{e}}_{\mathit{r}}\mathrm{=}\frac{{r}}{\mathrm{\left|}{r}\mathrm{\right|}}\mathit{,}{{e}}_{\mathit{\theta }}\mathit{,}{{e}}_{\mathit{\phi }}$ forming an ortho-normal basis tangent to the lines of constant spherical coordinates r,θ,φ. Although v can be determined for closer objects, only its projection $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}{v}\mathrm{◦}{{e}}_{\mathit{\varrho }}$ onto the line of sight determined by the unit vector ${{e}}_{\mathit{\varrho }}\mathrm{=}\frac{{r}\mathrm{-}{{r}}_{\mathrm{\odot }}}{\mathrm{|}{r}\mathrm{-}{{r}}_{\mathrm{\odot }}\mathrm{|}}$ can be measured for all objects. This is the only kinematical information available at large distances, suitable for constraining the total Galactic mass. It is connected with direct measurements of the LSR relative velocity v_ϱ along the direction e_ϱ through the relation $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}{\mathit{v}}_{\mathit{\varrho }}\mathrm{+}{{v}}_{\mathrm{\odot }}\mathrm{◦}{{e}}_{\mathit{\varrho }}$.

Assuming a β(r), we can relate $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ to $\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ through the following identity true both for r<r_⊙ and r>r_⊙: $\begin{array}{ccc}& & \mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\left(\mathrm{1}\mathrm{-}\frac{\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{4}}\mathit{H}\left(\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{\odot }}\right)\right)\\ & & \mathit{H}\mathrm{(}{\mathit{x}}^{\mathrm{)}}\mathrm{=}\mathrm{1}\mathrm{+}{\mathit{x}}^{-2}\mathrm{-}\frac{{\mathrm{(}{{\mathit{x}}^{\mathrm{2}}\mathrm{-}\mathrm{1}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathrm{2}\hspace{0.17em}{\mathit{x}}^{\mathrm{3}}}\ln \left|\frac{\mathit{x}\mathrm{+}\mathrm{1}}{\mathit{x}\mathrm{-}\mathrm{1}}\right|\mathrm{·}\end{array}$(5)For the isotropic dispersion, β(r) = 0, we can formally make the (incorrect) identification $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{v}}_{\mathit{r}}\end{array}\mathrm{=}{\mathit{v}}_{\mathit{r}}$ globally, without making any error in equating the resulting dispersions $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{⟩}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$. In general, we can only infer $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\mathrm{~}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ at large radii if β is asymptotically bound (which may not hold for nearly circular orbits). This is because $\mathit{H}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{~}\frac{\mathrm{8}}{\mathrm{3}{\mathit{x}}^{\mathrm{2}}}$ asymptotically, hence $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\mathrm{~}\left(\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathrm{2}{\mathit{r}}_{\mathrm{\odot }}^{\mathrm{2}}}{\mathrm{3}{\mathit{r}}^{\mathrm{2}}}\right)\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ for r large enough. The relation Eq. (5) was given by Dehnen et al. (2006, Eq. (6) therein) and by Battaglia et al. (2005, the correct version can be found in the erratum Battaglia et al. 2006). For a mathematical completness we present our independent derivation of this formula in Appendix A.

While determining a PDF f = f(r,v) from the $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ observable, a self-consistent β(r) may be looked for by iterations. The first recursion step makes the assumption $\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ as if β(r) = 0, and a first approximation to f is obtained, from which a β(r) prediction for the next iteration step is calculated. Substituted in Eq. (5), the β(r) gives rise to a new $\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$. The process is repeated until a stable β(r) is reached. However, the distinction between $\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ and $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ is practically unimportant unless the lower radii region is considered. In preparing the RVD profile below, we may neglect this distinction.

Fig. 1 Profile of RVD ${\mathit{G}}^{-1}\mathrm{⟨}\mathit{r}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$ for tracers with $\frac{\mathrm{1}}{\mathrm{2}\mathit{G}}\mathit{r}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\lesssim \mathrm{3.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$. The horizontal bars represent the effective radial bin size of the moving average. The vertical bars indicate the spread in the profile due to the inclusion/exclusion of random subsets of tracers. A detailed description of how this profile was obtained is given in Bratek et al. (2014). For comparison, the grey solid circles show the ${\mathit{G}}^{-1}\mathit{r}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$ values calculated based on the velocity dispersion points determined in Xue et al. (2008) and Deason et al. (2012b).

2.1.2. Measurements data

In our previous work (Bratek et al. 2014), we determined a RVD profile (Fig. 1), which we now assume as the basis for generating the initial conditions for the simulation in Sect. 3.1, and use it as a reference profile for comparison with the simulation results in Sect. 3.4. In preparing this profile, we assume the Sun’s distance to the Galactic centre R_° = 8.5 ± 0.4 kpc, the local disk rotation speed V_° = 240 ± 16 km s^-1 (Bovy et al. 2009; Koposov et al. 2010) and the components of the Sun’s velocity vector with respect to the LSR (U,V,W) = (11.1 ± 1.7,12.24 ± 2.5,7.25 ± 0.9) km s^-1 (Schönrich et al. 2010). We used the following position-velocity data: the halo giant stars (Dohm-Palmer et al. 2001; Starkenburg et al. 2009) from the Spaghetti Project Survey (Morrison et al. 2000); the blue horizontal branch stars (Clewley et al. 2004) from the United Kingdom Schmidt Telescope observations and SDSS; the field horizontal branch and A-type stars (Wilhelm et al. 1999) from the Beers et al. (1992) survey; the globular clusters (Harris 1996) and the dwarf galaxies (Mateo 1998). The data were recalculated to epoch J2000 when necessary. In addition, we included the ultra-faint dwarf galaxies, such as Ursa Major I and II, Coma Berenices, Canes Venatici I and II, Hercules (Simon & Geha 2007), Bootes I, Willman 1 (Martin et al. 2007), Bootes II (Koch et al. 2009), Leo V (Belokurov et al. 2008), Segue I (Geha et al. 2009), and Segue II (Belokurov et al. 2009). To eliminate a possible decrease in the RVD at lower radii due to circular orbits in the disk, we excluded tracers in the neighbourhood (R/ 20)² + (Z/ 4)²< 1 (in units of kpc) of the mid-plane. We also did not take into account: a) a distant Leo T, located at r> 400 kpc; b) Leo I, rejected for reasons largely discussed in Bratek et al. (2014); c) a single star for which $\mathit{r}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{G}\mathrm{\right)}\mathit{>}\mathrm{5.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$; and d) four additional objects for which $\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{G}\mathrm{\right)}\gtrsim \mathrm{3.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$ (these are: 88-TARG37, Hercules, J234809.03-010737.6 and J124721.34+384157.9). As shown with the help of a simple asymptotic estimator (Bratek et al. 2014), had we not excluded d) the total expected mass would have been increased by only a factor of ≈1.16.

3. Simulation of RVD in a background field

We model the Galactic potential as consisting of: a disk-like component (accounting for the Galactic rotation curve) and a hot gaseous halo, Ψ = Ψ_disk + Ψ_gas, of total mass M_ref = 1.8 × 10¹¹ M_⊙. As a starting point for further analyses, we construct an initial PDF by applying the method of Sect. 2 to the RVD profile in Fig. 1, assuming M = M_ref, which is close to the lower bound for this RVD. In Fig. 2 we show the resulting PDF on the (e,ϵ) plane.

Fig. 2 A distribution function f(e,ϵ) consistent with the RVD profile shown in Fig. 1. The function was obtained with the help of the Keplerian ensemble method, assuming Mref = 1.8 × 1011 M⊙, Rua = 18 kpc, and Rub = 240 kpc. The contour plot shows (f(e,ϵ) /fS)1 / 10, with fS being the maximum value of f(e,ϵ) on the triangular domain.

3.1. Setting the initial conditions

The first stage towards determining the initial conditions corresponding to the PDF f(e,ϵ) shown in Fig. 2 involves generating a random set of initial radii ${\mathrm{ℐ}}_{\mathrm{0}}\mathrm{=}\mathrm{\{}{\mathit{u}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$ in the range u_a<u_i<u_b (N is the number of all test bodies), and with the number density ν_u [ f ] from Eq. (2). This task can be done with the help of the inverse cumulative probability function $\mathit{\chi }\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{\equiv }{{}^{\mathrm{\int }}\mathit{u}_{\mathit{a}}}_{\mathit{u}}^{}\mathit{\nu }\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}$ of the probability density ν(u).

In the next stage, we need to ensure the spherical symmetry of the initial state. We assign to ℐ₀ a set $\mathrm{\left\{}\mathrm{\right(}{\mathit{\theta }}_{\mathit{i}}\mathit{,}{\mathit{\phi }}_{\mathit{i}}\mathrm{)}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$ of spherical coordinates of directions uniformly distributed on the unit sphere. This gives us the initial positions ${\mathrm{ℐ}}_{\mathrm{1}}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}{\mathit{u}}_{\mathit{i}}\mathit{,}{\mathit{\theta }}_{\mathit{i}}\mathit{,}{\mathit{\phi }}_{\mathit{i}}\mathrm{)}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$.

Next, to obtain the initial velocities, we choose random parameters (e,ϵ,ψ) consistently with the initial PDF, assigning to each (u_k,θ_k,φ_k) ∈ ℐ₁ an elliptic orbit; a particular test body would follow in the point mass field. To this end, we consider triples of random numbers (e,ϵ,X) uniformly distributed in their respective range: e ∈ (0,1), ϵ ∈ (0,1 / (2u_a)) and X ∈ (0,f_S), with X being an auxiliary variable and f_S = max { f(e,ϵ):(e,ϵ) ∈ S }. For each u_k ∈ ℐ₀, we carry on generating random triples (e,ϵ,X) until we encounter one (labelled with a subscript k) for which both X<f(e_k,ϵ_k) and (e_k,ϵ_k) ∈ S(u_k). This procedure yields a set of random pairs ${\mathrm{ℐ}}_{\mathrm{2}}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}{\mathit{e}}_{\mathit{i}}\mathit{,}{\mathit{ϵ}}_{\mathit{i}}\mathrm{)}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$ with a non-uniform number density distribution f(e,ϵ) and each confined to a u_i-dependent region S(u_i). To each u_k, we also assign its respective random angle ψ_k uniformly distributed in the range (0,2π ] and fixing the plane of the corresponding ellipse. In effect, we obtain a set ${\mathrm{ℐ}}_{\mathrm{2}}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}{\mathit{e}}_{\mathit{i}}\mathit{,}{\mathit{ϵ}}_{\mathit{i}}\mathit{,}{\mathit{\psi }}_{\mathit{i}}\mathrm{)}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$ and form the set $\mathrm{ℐ}\mathrm{=}{\mathrm{ℐ}}_{\mathrm{1}}\mathrm{\oplus }{\mathrm{ℐ}}_{\mathrm{2}}\mathrm{=}\mathrm{\left\{}\mathrm{\right(}{\mathit{u}}_{\mathit{i}}\mathit{,}{\mathit{\theta }}_{\mathit{i}}\mathit{,}{\mathit{\phi }}_{\mathit{i}}\mathit{,}{\mathit{e}}_{\mathit{i}}\mathit{,}{\mathit{ϵ}}_{\mathit{i}}\mathit{,}{\mathit{\psi }}_{\mathit{i}}\mathrm{)}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}$. Finally, by applying the transformation Eq. (1) to each element of ℐ, we obtain the required set of initial positions and velocities in spherical coordinates, leading to an initial randomly generated RVD overlapping well with that in Fig. 1 in the region of interest.

Fig. 3 A sequence A,B,...S,T of simulated RVD models shown at distinct simulation instants with red lines. The models were obtained from the blue line RVD (evolved from the initial PDF in the Ψdisk + Ψgas potential of total mass Mref) by rescaling the horizontal and vertical directions so as to overlap with the black line reference RVD (from measurements) as good as possible. The blue line PDF does not differ much from the violet line PDF, which assumes β = 0. The light gray region is the RVD uncertainty defined by the vertical bars in Fig. 1. The bottom figure in each panel shows a β(r) profile corresponding to its respective RVD model. The decrease in β(r) towards negative values in the lower radii region r< 40 kpc is a model artifact discussed in the text.

3.2. Gravitational potential

The Ψ_disk part of Ψ is described by the thin disk model ${\mathrm{\Psi }}_{\mathrm{disk}}\mathrm{\left(}\mathit{\rho ,\zeta }\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{4}\mathit{G}\hspace{0.17em}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\hspace{0.17em}\frac{\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\hspace{0.17em}\mathit{K}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathit{\sigma }\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}\mathrm{\right)}}{\sqrt{\mathrm{(}\mathit{\rho }\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathit{\rho }\end{array}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{\zeta }}^{\mathrm{2}}}}\mathit{,}$(6)with σ(ρ) being the column mass density of a finite-width disk found by recursions from the Galactic rotation curve in Jałocha et al. (2014). Here, $\mathit{k}\mathrm{=}\sqrt{\frac{\mathrm{4}\hspace{0.17em}\mathit{\rho }\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}}{\mathrm{(}\mathit{\rho }\mathrm{+}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{\zeta }}^{\mathrm{2}}}}$ and K is the elliptic integral of the first kind defined in Gradshtein et al. (2007). Most of the mass is enclosed within the inner disk ρ< 20 kpc: M₂₀ = 1.49 × 10¹¹ M_⊙, while M₃₀ = 1.51 × 10¹¹ M_⊙. The outer ρ> 30 kpc disk’s contribution to Eq. (6) is thus negligible and we can limit the integration to $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{\in }\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{30}\mathrm{\right)}\mathrm{kpc}$. To reduce the computation time, we tabulated the integral Eq. (6) at mesh points { ρ_j,ζ_k }, obtaining a smooth Ψ_disk by means of the interpolating series $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{\Psi }}_{\mathrm{disk}}\end{array}\mathrm{\left(}\mathit{\rho ,\zeta }\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{p,q,r}}{\mathit{\omega }}_{\mathit{pqr}}\hspace{0.17em}{\mathit{\rho }}^{\mathit{p}}{\mathit{\zeta }}^{\mathit{q}}\mathrm{(}{\mathit{\rho }}^{\mathrm{2}}\mathrm{+}{\mathit{\zeta }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{-}\mathit{r}\mathit{/}\mathrm{2}}$ with the coefficients ω_pqr found by the least-squares method, minimizing the discrepancy between Eq. (6) and the series evaluated at the mesh points. Within the desired accuracy, we found this approximation procedure to be numerically more efficient than the usual two-dimensional interpolation.

Based on the O_VIIK_α absorption-line strengths in the spectra of galactic nuclei and galactic sources, Gupta et al. (2012) found large amounts of baryonic mass in the form of hot gas surrounding the Galaxy. Assuming a homogeneous sphere model, they found the electron density n_e of 2.0 × 10^-4/ cm³ and the path length L of 72 kpc. Among other parameters, the total mass of the gas depends on the gas metalicity and the oxygen-to-helium abundance. For a reasonable set of parameters, they found the total mass to be 1.2−6.1 × 10¹⁰ M_⊙. We may assume M_gas = 3.0 × 10¹⁰ M_⊙ consistent with these values. More recently, applying the same observational method, Miller & Bregman (2013) found the mass function M(r) of the circum-galactic hot gas using a modified density profile $\begin{array}{ccc}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{0}}{\mathrm{(}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathit{c}}{\mathrm{)}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{-}\mathrm{3}\mathit{\lambda }\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$We use it as the source of the spherical component Ψ_gas(r), with the parameters n₀ = 0.46 cm^-3, r_c = 0.35 kpc and λ = 0.58 allowable by the best fit to the measurements. Then the integrated mass is M_gas = 3 × 10¹⁰ M_⊙ at r = 100 kpc. For λ ≤ 1, the mass function is divergent and the integration must be cut off at some radius, which is to some extent arbitrary. The cutoff at 100 kpc falls within the limits 18 kpc and 200 kpc on the minimum and maximum mass of the halo considered in Miller & Bregman (2013).

3.3. Numerical solution of the equations of motion

We consider a test particle of mass m in cylindrical coordinates (ρ,ϕ,ζ), moving in an axi-symmetric gravitational field described by the potential Ψ. In this symmetry, the angular momentum component J_ϕ = mρ²ϕ′(t) is conserved. On account of ϕ being a monotone function of the time t for orbits with J_ϕ ≠ 0, we may regard ϕ as the independent parameter. In this parametrization, the equations of motion reduce to $\begin{array}{ccc}{\mathit{\rho }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathrm{𝒥}}{\mathit{v}}_{\mathit{\rho }}\mathit{,} {\mathit{v}}_{\mathit{\rho }}^{\mathrm{\prime }}\mathit{\varphi }\mathrm{=}\frac{\mathrm{𝒥}}{\mathit{\rho }}\mathrm{-}\frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathrm{𝒥}}\hspace{0.17em}{\mathit{\partial }}_{\mathit{\rho }}\mathrm{\Psi }\mathrm{\left(}\mathit{\rho ,\zeta }\mathrm{\right)}\mathit{,}\\ {\mathit{\zeta }}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\varphi }\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathrm{𝒥}}{\mathit{v}}_{\mathit{\zeta }}\mathit{,} {\mathit{v}}_{\mathit{\zeta }}^{\mathrm{\prime }}\mathit{\varphi }\mathrm{=}\mathrm{-}\frac{{\mathit{\rho }}^{\mathrm{2}}}{\mathrm{𝒥}}\hspace{0.17em}{\mathit{\partial }}_{\mathit{\zeta }}\mathrm{\Psi }\mathrm{\left(}\mathit{\rho ,\zeta }\mathrm{\right)}\mathit{,}\end{array}$(7)with being the angular momentum per unit mass, and v_ρ and v_ζ the velocity variables in the ortho-normal basis of the coordinate lines ρ,ζ. We solve Eq. (7) numerically by using a 4-order Runge-Kutta method with adaptive step size, controlled so as to keep below some small threshold value the relative change |Δℰ / ℰ| in the energy per unit mass $\mathrm{ℰ}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{}^{\mathrm{\right(}}{\mathit{v}}_{\mathit{\rho }}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\zeta }}^{\mathrm{2}}\mathrm{+}{\mathrm{𝒥}}^{\mathrm{2}}\mathit{/}{\mathit{\rho }}^{\mathrm{2}}^{\mathrm{\left)}}\mathrm{+}\mathrm{\Psi }\mathrm{\left(}\mathit{\rho ,\zeta }\mathrm{\right)}$. The relative change of energy along each trajectory during our simulation is then always smaller than 10^-6: |ℰ(t) / ℰ(0) − 1| < 10^-6 for all t, and this precision suffices for the purposes of this work.

3.4. Results and discussion

Using the numerical procedure of Sect. 3.3, we obtained 3665 trajectories of test bodies starting from the initial conditions of Sect. 3.1 and bound in the potential Ψ_disk + Ψ_gas defined in Sect. 3.2. The initial state agrees with the initial PDF (Fig. 2) of a stationary solution of Jeans’ problem in the point mass potential and is consistent with the reference RVD (Fig. 1). Using these trajectories, we determined the RVD evolution from the initial one and track it through a sequence of snapshots taken at various instants, as shown in Fig. 3 (with a step size of ≈1 Gyr). Each snapshot can be regarded as an independent RVD model used to estimate the Galaxy mass by comparing the evolved RVD with the reference RVD. In this approach, the mass we assign to Ψ_disk + Ψ_gas becomes a function of the simulation time, while the extent of that time has no physical meaning.

At this point, it is appropriate to bring attention to some features of the initial PDF that persist as model artifacts during the simulation. Namely, in the lower radii region, the evolved RVD values are reduced relative to the reference RVD. The first reason is that for the initial PDF from the point mass approximation, the velocities are too high in a fraction of objects in the modified potential, owing to a more extended mass distribution, and either quickly populate more distant regions or are not bound (in preparing the evolved RVD profiles we considered only bound trajectories). Consequently, the higher velocity values do not contribute in this region and the RVD values are reduced. The other model artifact is due to a cutoff in the PDF domain introduced in the Keplerian ensemble method to automatically prevent test bodies from penetrating the interior of a central spherical region where the point mass approximation is violated. As so, there is no limit to the number of almost nearly circular orbits that the external neighbourhood of this region can accommodate. Too many elongated orbits cannot occur there for geometrical reasons, while the admissible elongated orbits only enter this region with their pericentric sides (where radial motions are nearly vanishing). In consequence, the overall mean RVD in this neighbourhood is reduced below the observed values. Because the initial PDF has been identified with the PDF of the Keplerian ensemble, a qualitatively similar reduction mechanism in the evolved RVD comes about in the modified potential, reflecting in the β(r) reduced towards more negative values in the lower radii region. However, more circular motions in this region could be interpreted as consistent with a contribution from a cold disk. With a better initial PDF this model effect could be eliminated, but it seems of no importance for the accuracy of the total mass determination for which the region of greater radii is more important. A similar cutoff mechanism may increase the number density in the neighbourhood of the upper boundary (which we assumed to be of 240 kpc). Namely, the high RVD values observed for moderate radii, and modelled in the point mass approximation by more elongated elliptic orbits, are reduced to zero again close to the upper boundary. This reduction in the RVD appears naturally; the elongated orbits enter this region with their apocentric sides where radial motions are almost vanishing and where test bodies spent a relatively longer time, and this effect can be amplified by increasing the number density of bodies on more circular orbits. Because the model RVD is compared with the measurements at moderate radii, this effect can again be neglected.

Fig. 4 Mass estimator μMref as a function of the simulation time, with factors μ obtained from best-fit model RVD at various threshold radii RT.

Now, we return to the main issue. As mentioned earlier, we want to verify the expectation that the evolved PDF should be in a sense close to the initial PDF, independent of the simulation instant if the point mass approximation well describes the real situation at higher distances. If the RVD evolved in the modified potential turned out to be collapsing to much smaller values or change its shape completely, then this would mean that the mass estimate based on the initial PDF was wrong and inconsistent with the new evolved PDF.

As seen in Fig. 3, although the evolved RVD changes with the simulation time, it decreases a little in some regions and grows again later. The RVD generally remains high in the larger radii region. Similarly, the characteristic maxima in the initial RVD are not destroyed, but they oscillate. Besides the evolved RVD profile (blue line) corresponding to the mass M_ref, in each snapshot we also shown a corrected RVD (red line). We consider this corrected RVD our model RVD profile, obtained by multiplying M_ref and the radial variable with suitable factors close to unity, respectively μ and α, so as to make the model RVD coincide with the reference RVD as well as possible in the sense of the least squares³. During the simulation run, the length factor α varied in the range (0.85;1.02) with the mean 0.92 ± 0.03, while the mass factor μ varied in the range (1.03;1.37) with the mean 1.18 ± 0.07 (see Fig. 4). This gives the total mass estimate of (2.12 ± 0.13) × 10¹¹ M_⊙, oscillating in the range (1.85;2.47) × 10¹¹ M_⊙.

4. Conclusions

The lower bound for the Galaxy mass of ≈2.1 × 10¹¹ M_⊙ obtained within the Keplerian ensemble framework is sufficient to preserve the values and qualitative features of the RVD profile from halo tracers within 150 kpc during a numerical simulation run in a modified potential. In this sense the numerically evolved RVD is stable. These results also substantiate structural stability of the point mass approximation (with more general unconstrained PDF’s), showing that the lower range for Galaxy mass estimates is reliable. A possible correction factor 1.16 to account for the four halo objects rejected in Sect. 2.1.2 would give a value of (2.5 ± 0.2) × 10¹¹, consistent with earlier estimates of 2.4 × 10¹¹ M_⊙ (Merrifield 1992; Little & Tremaine 1987) or with the value of (2.6−2.7) × 10¹¹ M_⊙ recently inferred from the kinematics of the Orphan stream (Newberg et al. 2010; Sesar et al. 2013) within ~60 kpc.

The crucial role in our analysis is played by the general unconstrained phase space. We stress that the phase-space model is not less important than the mass model, and focussing more attention on generic phase spaces might help to reduce the missing mass problem.

The dark matter halo is the most hypothetical and less constrained Galactic component. The size of the halo is usually given in terms of the virial radius defined as a radius of the sphere, which has an average density larger, by a factor of Δ = 200, than the critical density of the Universe (see e.g. McMillan 2011, because the convention for the parameter Δ is not unique). This leads to the virial mass of the halo of the order of 10¹² M_⊙. On the other hand, in Sect. 1 we argued that the dark matter halo with a low mass of about the lower bound M ≈ 2 × 10¹¹ M_⊙ and the size of about r ≈ 20 kpc is also allowed by the RVD profile observations, whereas, as we mentioned in Sect. 1, in the literature the mass estimates based on the RVD profile Galaxy differ from each other by a factor larger than two. This shows that the total mass of the halo is poorly constrained by the motions of the distant tracers.

We present here an extreme example of a model without NDM halo, and we have shown the model to be stable in a sense that it accounts for the measured RVD at each simulation instant. This shows that our model can be thought of as a collisionless system close to a steady state. The possibility of accounting for the RVD observations without NDM halo shows that either the halo is unnecessary for the understanding of the motions of the kinematical tracers, or that other observational features (e.g. the measurements of the β(r) function) are needed to define constraints on the phase space, which would allow to disambiguate between various halo mass profiles.

References

Appendix A: Derivation of a relation between the observable $\mathrm{⟨}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{v}\end{array}\begin{array}{c}\mathrm{2}\\ \mathit{r}\end{array}\mathrm{⟩}$ and the radial dispersion $\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$

Let v be the Galacto-centric velocity vector. Expressed in terms of its radial V_r and transversal components V_θ,V_φ it reads ${v}\mathrm{=}{}^{\mathrm{[}}\mathit{V}_{\mathit{r}}\sin \mathit{\theta }\cos \mathit{\phi }\mathrm{+}{\mathit{V}}_{\mathit{\theta }}\cos \mathit{\theta }\cos \mathit{\phi }\mathrm{-}{\mathit{V}}_{\mathit{\phi }}\sin \mathit{\phi ,}\mathrm{(}{\mathit{V}}_{\mathit{r}}\sin \mathit{\theta }{\mathrm{+}}^{}$V_θcosθ)sinφ + V_φcosφ,V_rcosθ − V_θsinθ^]. As e_r = $\frac{{r}}{\left|{r}\right|}\mathrm{=}\left[\sin \mathit{\theta }\cos \mathit{\phi ,}\sin \mathit{\theta }\sin \mathit{\phi ,}\cos \mathit{\theta }\right]$ and r_⊙ = [r_⊙,0,0] the l.o.s versor is ${{e}}_{\mathit{\rho }}\mathrm{\equiv }\frac{{r}\mathrm{-}{{r}}_{\mathrm{\odot }}}{\left|{r}\mathrm{-}{{r}}_{\mathrm{\odot }}\right|}\mathrm{=}\frac{\left[\mathit{r}\sin \mathit{\theta }\cos \mathit{\phi }\mathrm{-}{\mathit{r}}_{\mathrm{\odot }}\mathit{,r}\sin \mathit{\theta }\sin \mathit{\phi ,r}\cos \mathit{\theta }\right]}{\sqrt{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathrm{\odot }}^{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{r}{\mathit{r}}_{\mathrm{\odot }}\sin \mathit{\theta }\cos \mathit{\phi }}}$. We take the mean values $\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{\equiv }\mathrm{⟨}{{}^{\mathrm{\right(}}{{e}}_{\mathit{r}}\mathrm{◦}{v}^{\mathrm{\left)}}}^{\mathrm{2}}\mathrm{⟩}$ and $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{⟩}\mathrm{\equiv }\mathrm{⟨}{{}^{\mathrm{\right(}}{{e}}_{\mathit{\rho }}\mathrm{◦}{v}^{\mathrm{\left)}}}^{\mathrm{2}}\mathrm{⟩}$ over thin spherical shells and consider them as functions of the Galacto-centric distance r. For a spherically symmetric system, we define $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}$ as $\begin{array}{ccc}\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }\hspace{0.17em}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\hspace{0.17em}\mathrm{d}\mathit{\theta }\sin \mathit{\theta }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\hspace{0.17em}\mathrm{d}\mathit{\phi }\hspace{0.17em}{⟨{\mathrm{(}{{{e}}_{\mathit{\rho }}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathrm{◦}{v}}^{\mathrm{)}}}^{\mathrm{2}}⟩}_{\mathrm{int}}\mathit{.}& & \end{array}$Here, ⟨ ·⟩ _int is the averaging over the velocities weighted by a spherically symmetric PDF f(r,v(r)), normalized so as ν(r) ⟨ (·) ⟩ _int ≡ ^∫(·)f(r,v(r)) d³v, with ν denoting the number density. The scalar product squared ${{}^{\mathrm{\right(}}{{e}}_{\mathit{\rho }}\mathrm{\left(}\mathit{r,\theta ,\phi }\mathrm{\right)}\mathrm{◦}{v}^{\mathrm{\left)}}}^{\mathrm{2}}$ is a homogenuous form of second degree in the velocities V_r,V_θ,V_φ with coefficients being functions of r,θ,φ. With a direct inspection, one can notice that the integration over θ,φ of the coefficients standing at V_rV_θ, V_θV_φ, V_φV_r, gives zero (the velocity products are independent of θ,φ). Thus, upon integration over the velocities, we can focus only on the terms involving dispersions $\mathrm{⟨}{\mathit{V}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, $\mathrm{⟨}{\mathit{V}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, $\mathrm{⟨}{\mathit{V}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}$. Furthermore, it also follows from the spherical symmetry that $\mathrm{⟨}{\mathit{V}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{⟨}{\mathit{V}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and, trivially, that the ratios $\mathrm{⟨}{\mathit{V}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{⟩}\mathit{/}\mathrm{⟨}{\mathit{V}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$, $\mathrm{⟨}{\mathit{V}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{⟩}\mathit{/}\mathrm{⟨}{\mathit{V}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$ define the same function of r. In accordance with the common convention in the theory of spherical Jeans equations, we express this function in terms of the flattening of the dispersion ellipsoid, β(r). Then, by making the substitution $\mathrm{⟨}{\mathit{V}}_{\mathit{\theta }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{(}\mathit{r}\mathrm{\left)}\mathrm{\right)}\mathrm{⟨}{\mathit{V}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and $\mathrm{⟨}{\mathit{V}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{(}\mathit{r}\mathrm{\left)}\mathrm{\right)}\mathrm{⟨}{\mathit{V}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, we obtain that $\begin{array}{ccc}\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{\left(}\mathit{r}\mathrm{\right)}\left(\mathrm{1}\mathrm{-}\frac{\mathit{\beta }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{r}}_{\mathrm{\odot }}^{\mathrm{2}}}\mathrm{·}\mathrm{ℐ}\mathrm{\left(}\mathit{\alpha }\mathrm{\right(}\mathit{r}\mathrm{\left)}\mathrm{\right)}\right)\mathit{,}& & \end{array}$where $\mathit{\alpha }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\equiv }\frac{\mathrm{2}\mathit{r}{\mathit{r}}_{\mathrm{\odot }}}{{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{r}}_{\mathrm{\odot }}^{\mathrm{2}}}\mathit{<}\mathrm{1}$ for r ≠ r_⊙ and $\begin{array}{ccc}\mathrm{ℐ}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\hspace{0.17em}\mathrm{d}\mathit{\theta }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\hspace{0.17em}\mathrm{d}\mathit{\phi }\hspace{0.17em}\hspace{0.17em}\frac{\sin \mathit{\theta }\mathrm{(}{{\cos }^{\mathrm{2}}\mathit{\theta }{\cos }^{\mathrm{2}}\mathit{\phi }\mathrm{+}{\sin }^{\mathrm{2}}\mathit{\phi }}^{\mathrm{)}}}{\mathrm{1}\mathrm{-}\mathit{\alpha }\sin \mathit{\theta }\cos \mathit{\phi }}\mathrm{·}& & \end{array}$We recall that all integrals that are zero by symmetries have already been omitted in the expression for ℐ(α). Because the requirements for the integration of a functional series term by term and its limit are met for 0 ≤ α< 1, the integral ℐ(α) can be calculated by a Taylor series expansion in α (note that owing to the vanishing of the integrals ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\cos }^{\mathit{m}}\mathit{\phi }\hspace{0.17em}\mathrm{d}\mathit{\phi }$ with odd m, only even powers of α are present in the series). On reducing the summands with the help of the Pythagorean trigonometric identity, the remaining non-zero coefficients in the power series in α arrange to products of elementary definite integrals $\mathrm{ℐ}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\sum }_{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{+}\mathrm{\infty }}{\mathit{\alpha }}^{\mathrm{2}\mathit{n}}{}^{\mathrm{\right(}}{\mathit{S}}_{\mathit{n}}{\mathit{C}}_{\mathit{n}}\mathrm{-}{\mathit{S}}_{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{C}}_{\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{\left)}}$, where ${\mathit{S}}_{\mathit{n}}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\pi }}{\sin }^{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\theta }\mathrm{=}\mathrm{2}\frac{\mathrm{\left(}\mathrm{2}\mathit{n}\mathrm{\right)}\mathrm{!}}{\mathrm{(}\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{!}}$; ${\mathit{C}}_{\mathit{n}}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\cos }^{\mathrm{2}\mathit{n}}\mathit{\phi }\hspace{0.17em}\mathrm{d}\mathit{\phi }\mathrm{=}\mathrm{2}\mathit{\pi }\frac{\mathrm{(}\mathrm{2}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}{\mathrm{\left(}\mathrm{2}\mathit{n}\mathrm{\right)}\mathrm{!}}$. Now, ${\mathit{C}}_{\mathit{n}}{\mathit{S}}_{\mathit{n}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}$ and ${\mathit{C}}_{\mathit{n}\mathrm{+}\mathrm{1}}{\mathit{S}}_{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{3}}$. Hence, $\mathrm{ℐ}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\alpha }}{\sum }_{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{2}\hspace{0.17em}{\mathit{\alpha }}^{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{3}}}{}^{\mathrm{\right(}}\mathrm{-}\mathrm{2}\mathit{\alpha }\mathrm{+}{\sum }_{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}}{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{1}}^{\mathrm{\left)}}$, where we have subtracted the excess term 2α in the second series after renaming n → n + 1. Both of the infinite series are Taylor series expansions of $\ln \frac{\mathrm{1}\mathrm{+}\mathit{\alpha }}{\mathrm{1}\mathrm{-}\mathit{\alpha }}$, therefore, $\begin{array}{ccc}\mathrm{ℐ}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\alpha }}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{1}\mathrm{-}{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{2}{\mathit{\alpha }}^{\mathrm{3}}}\ln \left(\frac{\mathrm{1}\mathrm{+}\mathit{\alpha }}{\mathrm{1}\mathrm{-}\mathit{\alpha }}\right)\mathit{,} \mathrm{0}\mathrm{\le }\mathit{\alpha }\mathit{<}\mathrm{1.}& & \end{array}$Next, using the earlier expression for $\mathrm{⟨}\stackrel{}{{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathrm{⟩}$ and substituting the definition of α(r) in place of α, we finally obtain Eq. (5).

All Figures

Fig. 1 Profile of RVD ${\mathit{G}}^{-1}\mathrm{⟨}\mathit{r}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$ for tracers with $\frac{\mathrm{1}}{\mathrm{2}\mathit{G}}\mathit{r}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\lesssim \mathrm{3.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$. The horizontal bars represent the effective radial bin size of the moving average. The vertical bars indicate the spread in the profile due to the inclusion/exclusion of random subsets of tracers. A detailed description of how this profile was obtained is given in Bratek et al. (2014). For comparison, the grey solid circles show the ${\mathit{G}}^{-1}\mathit{r}\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}$ values calculated based on the velocity dispersion points determined in Xue et al. (2008) and Deason et al. (2012b).

In the text

	Fig. 2 A distribution function f(e,ϵ) consistent with the RVD profile shown in Fig. 1. The function was obtained with the help of the Keplerian ensemble method, assuming Mref = 1.8 × 1011 M⊙, Rua = 18 kpc, and Rub = 240 kpc. The contour plot shows (f(e,ϵ) /fS)1 / 10, with fS being the maximum value of f(e,ϵ) on the triangular domain.
In the text

Fig. 3 A sequence A,B,...S,T of simulated RVD models shown at distinct simulation instants with red lines. The models were obtained from the blue line RVD (evolved from the initial PDF in the Ψdisk + Ψgas potential of total mass Mref) by rescaling the horizontal and vertical directions so as to overlap with the black line reference RVD (from measurements) as good as possible. The blue line PDF does not differ much from the violet line PDF, which assumes β = 0. The light gray region is the RVD uncertainty defined by the vertical bars in Fig. 1. The bottom figure in each panel shows a β(r) profile corresponding to its respective RVD model. The decrease in β(r) towards negative values in the lower radii region r< 40 kpc is a model artifact discussed in the text.

In the text

	Fig. 4 Mass estimator μMref as a function of the simulation time, with factors μ obtained from best-fit model RVD at various threshold radii RT.
In the text