Reliability of the measured velocity anisotropy of the Milky Way stellar halo

TL;DR

This paper assesses the reliability of measuring the velocity anisotropy of the Milky Way's stellar halo, highlighting limitations of current methods at large distances and emphasizing the importance of future Gaia data for accurate estimates.

Contribution

It demonstrates through mock analyses that local fitting methods are reliable only within 15 kpc and discusses the challenges of measuring velocity anisotropy at larger radii.

Findings

Local fitting methods are reliable only within 15 kpc.

Estimating tangential velocities becomes difficult at larger radii.

Future Gaia data will improve measurements beyond 20 kpc.

Abstract

Determining the velocity distribution of halo stars is essential for estimating the mass of the Milky Way and for inferring its formation history. Since the stellar halo is a dynamically hot system, the velocity distribution of halo stars is well described by the 3-dimensional velocity dispersions $(\sigma_r, \sigma_\theta, \sigma_\phi)$, or by the velocity anisotropy parameter $\beta=1-(\sigma_\theta^2+\sigma_\phi^2)/(2\sigma_r^2)$. Direct measurements of $(\sigma_r, \sigma_\theta, \sigma_\phi)$ consistently suggest $\beta =0.5$-$0.7$ for nearby halo stars. In contrast, the value of $\beta$ at large Galactocentric radius $r$ is still controversial, since reliable proper motion data are available for only a handful of stars. In the last decade, several authors have tried to estimate $\beta$ for distant halo stars by fitting the observed line-of-sight velocities at each radius with simple velocity distribution models (local fitting methods). Some results of local fitting methods imply $\beta<0$ at $r \gtrsim 20 \;\rm{kpc}$, which is inconsistent with recent predictions from cosmological simulations. Here we perform mock-catalogue analyses to show that the estimates of $\beta$ based on local fitting methods are reliable only at $r \leq 15 \;\rm{kpc}$ with the current sample size ($\sim10^3$ stars at a given radius). As $r$ increases, the line-of-sight velocity (corrected for the Solar reflex motion) becomes increasingly closer to the Galactocentric radial velocity, so that it becomes increasingly more difficult to estimate tangential velocity dispersion $(\sigma_\theta, \sigma_\phi)$ from line-of-sight velocity distribution. Our results suggest that the forthcoming Gaia data will be crucial for understanding the velocity distribution of halo stars at $r \geq 20\;\rm{kpc}$.

Equations

$$\beta (r)=1-\frac{{\sigma }_{\theta }^2(r)+{\sigma }_{\varphi }^2(r)}{2{\sigma }_r^2(r)}.$$

$$\displaystyle f_{\rm Gauss}(\mbox{\boldmath$v$}|\mbox{\boldmath$x$})=\frac{1}{(\sqrt{2\pi})^{3/2}\sigma_{r}\sigma_{\theta}\sigma_{\phi}}$$

$$\times \exp \left[-\left(\frac{v_r^2}{2{\sigma }_r^2(r)}+\frac{v_{\theta }^2}{2{\sigma }_{\theta }^2(r)}+\frac{(v_{\varphi }-V_{\mathrm{rot}}(r){)}^2}{2{\sigma }_{\varphi }^2(r)}\right)\right].$$

$$\displaystyle P(v_{\rm los}|\mbox{\boldmath$x$},\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$$

$$\displaystyle=\frac{1}{\sqrt{2\pi}\sigma_{\rm los}(\mbox{\boldmath$x$})}\exp\left[-\frac{(v_{{\rm los},i}-V_{\rm rot}Q_{\phi}(\mbox{\boldmath$x$}))^{2}}{2\sigma_{\rm los}^{2}(\mbox{\boldmath$x$})}\right].$$

$$\displaystyle\sigma_{\rm los}(\mbox{\boldmath$x$})=\sqrt{\sigma_{r}^{2}(r)Q_{r}^{2}(\mbox{\boldmath$x$})+\sigma_{\theta}^{2}(r)Q_{\theta}^{2}(\mbox{\boldmath$x$})+\sigma_{\phi}^{2}(r)Q_{\phi}^{2}(\mbox{\boldmath$x$})}$$

$$\displaystyle\ln L=\ln\left[\prod_{i=1}^{N}P(v_{{\rm los},i}|\mbox{\boldmath$x$}_{i},\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})\right]$$

$$\displaystyle=-\frac{N}{2}\ln(2\pi)-\sum_{i=1}^{N}\left[\ln\sigma_{\rm los}(\mbox{\boldmath$x$}_{i})+\frac{(v_{{\rm los},i}-V_{\rm rot}Q_{\phi}(\mbox{\boldmath$x$}_{i}))^{2}}{2\sigma_{\rm los}^{2}(\mbox{\boldmath$x$}_{i})}\right].$$

$$\displaystyle P(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot}|\{(\mbox{\boldmath$x$}_{i},v_{{\rm los},i})\}_{i=1}^{N})$$

$$\displaystyle=\frac{P(\{v_{{\rm los},i}\}_{i=1}^{N}|\{\mbox{\boldmath$x$}_{i}\}_{i=1}^{N},\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})P(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})}{P(\{\mbox{\boldmath$x$}_{i},v_{{\rm los},i}\}_{i=1}^{N})}.$$

$$P_A({\sigma }_r,{\sigma }_{\theta },{\sigma }_{\varphi },V_{\mathrm{rot}})$$

$$\propto \{\begin{array}{l}1,({\sigma }_{\mathrm{tot}}<v_{\mathrm{esc}},|V_{\mathrm{rot}}|<v_0),\\ 0,(\mathrm{otherwise}).\end{array}$$

$$P_B({\sigma }_r,{\sigma }_{\theta },{\sigma }_{\varphi },V_{\mathrm{rot}})$$

$$\propto \{\begin{array}{l}{\sigma }_r^{-1},(1{\mathrm{km}\mathrm{s}}^{-1}<{\sigma }_{\mathrm{tot}}<v_{\mathrm{esc}},|V_{\mathrm{rot}}|<v_0),\\ 0,(\mathrm{otherwise}),\end{array}$$

$$P_C({\sigma }_r,{\sigma }_{\theta },{\sigma }_{\varphi },V_{\mathrm{rot}})$$

$$\propto \{\begin{array}{l}({\sigma }_r{\sigma }_{\theta }{\sigma }_{\varphi }{)}^{-1},(1{\mathrm{km}\mathrm{s}}^{-1}<{\sigma }_{\mathrm{tot}}<v_{\mathrm{esc}},|V_{\mathrm{rot}}|<v_0),\\ 0,(\mathrm{otherwise}).\end{array}$$

$${\sigma }_{\theta }^2+{\sigma }_{\varphi }^2=(1-\beta ){\sigma }_{r,\mathrm{true}}^2.$$

$$|\Delta {\beta }_{\mathrm{MaxL}}|=2\sqrt{\frac{8}{N}}\left(\frac{r^2}{R_0^2}-{\beta }_{\mathrm{true}}\right).$$

$${\sigma }_{\mathrm{los}}^2=Q_r^2{\sigma }_r^2+(1-Q_r^2){\sigma }_{\theta }^2,$$

$${\sigma }_{\mathrm{los}}^2(QMAX)={\sigma }_r^2,$$

$${\sigma }_{\mathrm{los}}^2(QMIN)={\sigma }_r^2+\frac{R_0^2}{r^2}({\sigma }_{\theta }^2-{\sigma }_r^2)={\sigma }_r^2\left[1-\frac{R_0^2}{r^2}{\beta }_{\mathrm{true}}\right],$$

$$\beta =\frac{r^2}{R_0^2}\left(1-\frac{{\sigma }_{\mathrm{los}}^2(QMIN)}{{\sigma }_{\mathrm{los}}^2(QMAX)}\right).$$

$$\Delta {\sigma }_{\mathrm{los}}^2(QMAX)=\frac{\sqrt{2}}{\sqrt{(N/2)}}{\sigma }_r^2,$$

$$\Delta {\sigma }_{\mathrm{los}}^2(QMIN)=\frac{\sqrt{2}}{\sqrt{(N/2)}}{\sigma }_r^2\left[1-\frac{R_0^2}{r^2}{\beta }_{\mathrm{true}}\right],$$

$$|\Delta \beta {|}_{\mathrm{ideal}}=\left(\frac{r^2}{R_0^2}\right){\left\{{\left(\frac{{\sigma }_{\mathrm{los}}^2(QMIN)}{{\sigma }_{\mathrm{los}}^4(QMAX)}\Delta {\sigma }_{\mathrm{los}}^2(QMAX)\right)}^2+{\left(\frac{\Delta {\sigma }_{\mathrm{los}}^2(QMIN)}{{\sigma }_{\mathrm{los}}^2(QMAX)}\right)}^2\right\}}^{1/2}=\sqrt{\frac{8}{N}}\left(\frac{r^2}{R_0^2}-{\beta }_{\mathrm{true}}\right).$$

$$|\Delta \beta |=2\sqrt{\frac{8}{N}}\left(\frac{r^2}{R_0^2}-{\beta }_{\mathrm{true}}\right)$$

$$f_{\mathrm{const}}(E,L)=A\exp \left[-(\alpha +2)\frac{E}{v_0^2}\right]\left(L_{\max }^2(E)-K{L}^2\right).$$

$$\beta (r)={\beta }_{\mathrm{const}}=\frac{K}{\frac{1}{2}(\alpha +2{)}^{7/2}{\alpha }^{-5/2}\exp (-1)-K}.$$

$$\displaystyle P(v_{\rm los}|\mbox{\boldmath$x$},K,\alpha)$$

$$\times \left\{{\alpha }^{1/2}(\alpha +2{)}^{5/2}-K{\alpha }^{3/2}(\alpha +2{)}^{1/2}\left[1+Q_r^2+(\alpha +2)(1-Q_r^2){\left(\frac{v_{\mathrm{los}}}{v_0}\right)}^2\right]\exp \left[1-{\left(\frac{v_{\mathrm{los}}}{v_0}\right)}^2\right]\right\}.$$

$$f_{\mathrm{OM}}(E,L)=A\exp \left[-\alpha \frac{E+L^2/(2r_a^2)}{v_0^2}\right]$$

$$\rho (r)=(2\pi {)}^{3/2}{v}_0^3{\alpha }^{-3/2}A{\left(\frac{r}{r_0}\right)}^{-\alpha }{\left(1+\frac{r^2}{r_a^2}\right)}^{-1}.$$

$$\displaystyle P(v_{\rm los}|\mbox{\boldmath$x$},r_{a},\alpha)=\frac{1}{\sqrt{2\pi}\sigma_{\rm los}}\exp\left[-\frac{v_{\rm los}^{2}}{2\sigma^{2}_{\rm los}}\right],$$

$${\sigma }_{\mathrm{los}}^2=\frac{(r_a^2+r^2{Q}_r^2)}{(r_a^2+r^2)}\cdot \frac{v_0^2}{\alpha }.$$

Full text

1. Introduction

The velocity distribution of halo stars provides a lot of useful information about the Milky Way. For example, by regarding the halo stars as the dynamical tracers, we can estimate the mass of the Milky Way through Jeans equation, if we can determine or assume their density $\rho$ and their 3-dimensional velocity dispersions $(\sigma_{r},\sigma_{\theta},\sigma_{\phi})$ (Dehnen et al., 2006; Gnedin et al., 2010). Also, since the stellar halo is a collisionless system, the orbital shapes of halo stars are relatively immune to adiabatic change of the gravitational potential. Thus, the current distribution of orbits can provide some insight into how the stellar halo was formed.

A useful quantity to describe the orbital distribution of halo stars is the velocity anisotropy parameter (Binney, 1980),

[TABLE]

By definition, $\beta=0$ corresponds to an isotropic velocity distribution. If radial or circular orbits dominate, $\beta$ is positive ($0<\beta<1$) or negative ($-\infty<\beta<0$), respectively. Since $\beta$ only depends on the velocity distribution at a given position and is independent of the potential, $\beta$ is a potentially powerful tool with which to compare the dynamical state of the Milky Way and that of simulated galaxies.

Interestingly, essentially all the recently published simulations of Milky Way-like galaxies based on $\Lambda$CDM cosmology show a qualitatively similar radial profile of $\beta$. Specifically, the $\beta$ profile of most simulated stellar haloes almost monotonically increases from $\beta\simeq 0$-$0.5$ near the Galactic center, through $\beta\simeq 0.3$-$0.7$ near the Solar circle, to $\beta\simeq 0.5$-$1$ at virial radius (Diemand et al., 2005; Sales et al., 2007). Importantly, most simulated haloes do not show negative value of $\beta$ outside the Solar circle.111 In a companion paper Loebman et al. (2016) discuss some exceptional situations where this general $\beta$ profile is not attained. This characteristic $\beta$ profile was also found in classical simulations of an initially cold stellar system (van Albada, 1982) that collapses and experiences violent relaxation (Lynden-Bell, 1967). The fact that most cosmological simulations predict a qualitatively similar profile of $\beta$ is intriguing, since the radial profile of $\beta$ may be used to test $\Lambda$CDM cosmology.

The observed value of $\beta=0.5$-$0.7$ of the halo stars in the Solar neighborhood is consistent with these simulations (Chiba & Yoshii, 1998; Chiba & Beers, 2000; Smith et al., 2009; Bond et al., 2010). The direct determination of $\beta$ at larger Galactocentric radii has been hampered by difficulty in obtaining accurate proper motion data, and it is only recently that Deason et al. (2013) obtained $\beta=0.0_{-0.4}^{+0.2}$ at $18\,\mathrm{kpc}<r<30\,\mathrm{kpc}$ by using 13 distant halo stars with proper motion data derived with the Hubble Space Telescope (see Cunningham et al. 2016 for their updated results). Their measurement of $\beta\simeq 0$ at large $r$ is intriguing, since it differs significantly from the characteristic $\beta(r)$ profile found in simulated galaxies.

Prior to these direct measurements of $\beta$, there were several attempts to estimate $\beta$ in the distant halo by using 3-dimensional positions and line-of-sight velocities. Broadly speaking, these methods can be classified into two methodologies: global fitting methods and local fitting methods.

The global fitting is based on global models of the stellar halo. In these methods, a functional form of the global model of the stellar halo is assumed. Then the likelihood of the model parameters given the line-of-sight velocity data is evaluated to determine the best-fit parameters. For example, Sommer-Larsen et al. (1994) and Sommer-Larsen et al. (1997) chose the functional forms of the radial and tangential velocity dispersion profiles that satisfy the spherical Jeans equation and fitted the 4-dimensional information of halo stars with these models to claim a declining profile of $\beta(r)$. The results from Sommer-Larsen et al. (1997) suggest that $\beta\simeq 0.5$ at $r=8\,\mathrm{kpc}$, $\beta\simeq 0$ at $r=20\,\mathrm{kpc}$, and $\beta\simeq-1.3$ at $r=50\,\mathrm{kpc}$. Also, Deason et al. (2012) used a distribution function model to fit the 4-dimensional information for halo stars located at $16\,\mathrm{kpc}<r<48\,\mathrm{kpc}$ and claimed that $\beta\simeq 0.5$. More recently, Williams & Evans (2015) analyzed a similar dataset with a more sophisticated action-based distribution function model and derived a radially varying $\beta$ profile with $\beta\simeq 0.4$ at $r=15\,\mathrm{kpc}$ and $\beta\simeq 0.65$ at $r=50\,\mathrm{kpc}$ – a result that is qualitatively similar to cosmological simulations, but rises more slowly with increasing radius.

The local fitting methods, on the other hand, interpret the line-of-sight velocity distribution of halo stars with a simple statistical model (e.g., Gaussian distribution) and estimate the velocity moments of halo stars. The idea behind these studies is simple: if the velocity distribution of halo stars is isotropic ($\beta=0$), we expect the line-of-sight velocity dispersion $\sigma_{\rm los}$ to be independent of the heliocentric line-of-sight direction. On the other hand, if $\beta\neq 0$, we expect the line-of-sight projection of the velocity ellipsoid gradually changes across the sky, due to the off-center location of the Sun in the Milky Way. In order to derive the velocity dispersion profile as a function of Galactocentric radius $r$, authors often divide their sample into several bins according to $r$ and perform the analyses for each radial bin. Therefore, these methods can be considered as a series of local fits of the velocity distribution. For example, Sirko et al. (2004) assumed that their sample of halo stars obeys a Gaussian velocity distribution and claimed $\beta\simeq 0$. Kafle et al. (2012) and King et al. (2015) used a similar formulation as in Sirko et al. (2004) for a larger sample of stars and claimed that $\beta<0$ at $r\gtrsim 20\,\mathrm{kpc}$. Also, Hattori et al. (2013) and Kafle et al. (2013) even claimed a metallicity-dependence of $\beta$ such that $\beta<0$ for relatively metal-poor halo stars and $\beta>0$ for relatively metal-rich halo stars.

It is currently unclear whether these discrepant estimates of $\beta$ are a result of the fitting method or difference in the sample of stars used. However, it is important to note that both global fitting and local fitting methods have their own advantages and disadvantages, and these methods are complementary to each other. In order to illustrate this point, we now compare the advantages and disadvantages associated with the (global) distribution function fitting and the local fitting methods.

When a distribution function model is used to fit a given sample of halo stars, is is generally assumed that the stellar halo is in dynamical equilibrium and the functional form of the distribution function is known. These basic assumptions have some advantages for the distribution function fitting. For example, these models are designed to be physical (e.g., non-zero phase-space density), so the best-fit solution for any given data set is guaranteed to be physical. Also, distribution function (combined with the potential model of the Milky Way) contains all the information needed to calculate any velocity moments. Thus, if we have some external knowledge on the stellar halo (such as the density profile estimated from other surveys) in addition to the kinematical data we are trying to fit, we can naturally incorporate additional information to improve the fit. However, the basic assumptions also have disadvantages: it is unclear if the stellar halo is really in dynamical equilibrium and if the assumed functional form is really adequate to model the stellar halo.

Local fitting methods, on the other hand, do not require the stellar halo to be in dynamical equilibrium, since we just need to fit the current velocity distribution at a given Galactocentric radius. In these methods, we need to assume the functional form of the velocity distribution (or at least some properties of the functions, such as the symmetry of the velocity distribution), but we can assume simple and flexible functions (e.g., Gaussian velocity distribution) to mitigate the arbitrariness of the chosen functions. By design, the best-fit model is not guaranteed to be in dynamical equilibrium (especially when the sample size is small), even if there is some external evidence to believe that the stellar halo is in dynamical equilibrium. However, the local fitting methods are useful if the stellar halo is not in dynamical equilibrium. For example, let us consider an idealized situation where the stellar halo is a sum of smooth component and a substructure, and suppose that the contribution from the substructure is only prominent at a certain Galactocentric radius. In this case, local fitting methods may detect the presence of substructure as an anomaly in the velocity dispersion profile, while global fitting methods which require smooth distribution functions may be adversely affected by this local substructure. Indeed, Kafle et al. (2012) used a local fitting method and claimed that the radial profile of $\beta(r)$ shows a dip-like structure at $r=17\,\mathrm{kpc}$. Although we will show in this paper that this dip might not be a real signature of the $\beta$ profile (due to the small sample size), at least it is safe to say that currently proposed stellar halo distribution functions do not accommodate this kind of dip; and hence these two classes of fitting methods are complementary to each other.

In this paper, we focus on the local fitting methods and evaluate how far out in the stellar halo are these methods can reliably estimate $\beta(r)$ from the currently available 4 dimensional data. To this end, we apply the local fitting methods to a set of mock catalogues and compare the estimated profile of $\beta$ with the input profile of $\beta$. Here we report that the widely-used local fitting methods can reliably recover $\beta$ within $r<15\,\mathrm{kpc}$, but can only weakly constrain $\beta$ at $r>15\,\mathrm{kpc}$, if 1000 sample stars are used for a given radius. This work is an extension of Hattori et al. (2013), who applied a matrix-based local fitting method to Sloan Digital Sky Survey (SDSS) data and pointed out that the estimation of $\beta$ from the line-of-sight velocity distribution can be biased at $r\gtrsim$ (16-18) $\,\mathrm{kpc}$. It is important to note that most of the above-mentioned works of local- and global-fitting methods (except for Sirko et al. 2004 and Hattori et al. 2013) do not present comparison with mock data to validate their methodology. Thus, their results at large $r$ might have been affected by systematic errors that were not accounted for.222 Williams & Evans (2015) confirmed that the observed line-of-sight velocity distribution is well reproduced by that of mock data generated from their best-fit model. However, this procedure is not enough to guarantee that the $\beta(r)$ profile of their best-fit model is similar to the $\beta(r)$ profile of the observed halo population.

The outline of this paper is as follows. In Section 2 we describe the maximum-likelihood and Bayesian formulations for inferring $\beta$ from 4-dimensional data. In Section 3 we describe our mock catalogues. In Section 4, we show the results of our maximum-likelihood analyses of our mock catalogues. In Section 5, we show the results of our Bayesian mock-analyses. Section 6 presents a discussion, and Section 7 sums up.

2. Method

Here we outline our maximum-likelihood and Bayesian formulations of the local fitting method to estimate the 3-dimensional velocity dispersion of halo stars from 4-dimensional phase-space coordinates (3D positions and line-of-sight velocities).

We note that we do not take into account any observational errors in our formulation (nor in our mock catalogues; see Section 3), since the main aim in this paper is to demonstrate how the performance of these widely-used local fitting methods deteriorates as the Galactocentric radius of the sample stars increases, even if we use idealized stellar data.

2.1. Maximum-likelihood method

Following Sirko et al. (2004) and Kafle et al. (2013), we assume that the distribution of velocity $v$ of halo stars at a given location $x$ takes the form333 This functional form includes the classical distribution function, the Osipkov-Merritt model in a singular isothermal potential (see Appendix D.2) as a special case.

[TABLE]

We adopt a spherical coordinate system $(r,\theta,\phi)$ such that r=|\mbox{\boldmathx}| is the Galactocentric radius, $\theta$ is the polar angle ($\theta=0$ corresponds to the Galactic disc plane), and $\phi$ is the azimuthal angle. The principal axes of the velocity ellipsoid are assumed to be aligned with this spherical coordinate system, which can be justified by the recent work of Evans et al. (2016). Also, 3-dimensional velocity dispersions $(\sigma_{r},\sigma_{\theta},\sigma_{\phi})$ and the mean azimuthal velocity $V_{\rm rot}$ are assumed to be functions of $r$ only. The last assumption implies a spherical density distribution of halo stars and a spherical potential of the Milky Way. This contradicts the claimed flattening of the stellar halo (Deason et al., 2011) and the potential (Koposov et al., 2010), although the density distribution may be less flattened in the outer part of the stellar halo (Carollo et al., 2007).

Under this velocity distribution model, the probability density that a star located at $x$ has the line-of-sight velocity $v_{\rm los}$ in the Galactic rest frame for a given set of parameters $(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ is expressed as (see Appendix A of Sirko et al. 2004 for derivation)

[TABLE]

Here,

[TABLE]

is the line-of-sight velocity dispersion of halo stars located at $x$. Also, (Q_{r},Q_{\theta},Q_{\phi})=(\mbox{\boldmathe}_{\rm los}\cdot\mbox{\boldmathe}_{r},\mbox{\boldmathe}_{\rm los}\cdot\mbox{\boldmathe}_{\theta},\mbox{\boldmathe}_{\rm los}\cdot\mbox{\boldmathe}_{\phi}) are dot products of the unit vector along the line-of-sight \mbox{\boldmathe}_{\rm los} and each of the unit vectors (\mbox{\boldmathe}_{r},\mbox{\boldmathe}_{\theta},\mbox{\boldmathe}_{\phi}) of the spherical coordinate system at $x$.

Suppose we have a sample of $N$ halo stars such that the location and line-of-sight velocity of $i$th star, (\mbox{\boldmathx}_{i},v_{{\rm los},i}), are known ($i=1,\cdots,N$) and that $N$ stars have an identical Galactocentric radius r=|\mbox{\boldmathx}_{i}|. Then, the total log-likelihood of the observational data given the parameters is expressed as

[TABLE]

The maximum-likelihood method (in the limit of no observational errors) finds the set of parameters that maximizes $\ln L$ at each radius $r$.

2.2. Bayesian method

Another useful method to estimate $\beta$ from 4-dimensional information is the Bayesian method (Kafle et al., 2012). The Bayesian method explicitly incorporates information about prior knowledge or constraints on parameters. The Bayesian method has the advantage that it makes it clear whether the data are constraining the parameters or the answers (posterior distribution) are driven by the priors. Here we briefly outline the formulation of the Bayesian method.

As in the maximum-likelihood method, let us assume that the distribution function of the stellar halo is given by equation (2). In Bayesian formulation, our aim is to obtain the posterior distribution of the model parameters given the data. In our case, the posterior distribution can be expressed as (via Bayes’ Theorem)

[TABLE]

Here the likelihood P(\{v_{{\rm los},i}\}_{i=1}^{N}|\{\mbox{\boldmathx}_{i}\}_{i=1}^{N},\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot}) is identical to $L$ in equation (5) and the evidence P(\{(\mbox{\boldmathx}_{i},v_{{\rm los},i})\}_{i=1}^{N}) can be regarded as a constant. Thus the only additional task for us is to set a certain prior distribution $P(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ of the model parameters.

In order to be as objective as possible in judging the performance of the Bayesian method, we use three types of relatively uninformative priors A, B and C as described below. We note that in the limit of an infinite number of sample stars with no error, the posterior distribution is expected to be independent of the choice of prior. However, in reality we only have a finite number of stars, so we need to use appropriate prior information in order to make the best use of the available data.

Prior A is a uniform prior for all the parameters $(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ given by

[TABLE]

Here, we define $\sigma_{\rm tot}\equiv\sqrt{\sigma_{r}^{2}+\sigma_{\theta}^{2}+\sigma_{\phi}^{2}}$. Also, $v_{0}(r)=220\,\mathrm{km\ s}^{-1}$ and $v_{\rm esc}(r)$ are the circular velocity and escape velocity at the Galactocentric radius $r$ of a truncated singular isothermal potential, respectively (although the details of the potential model do not affect the results).

From a mathematical point of view, prior A is not purely uninformative, since $(\sigma_{r},\sigma_{\theta},\sigma_{\phi})$ in our model are so-called scale parameters (while $V_{\rm rot}$ in our model is a so-called location parameter). The Jeffreys’ rule (Jeffreys 1961, Section 3.10; Ivezić et al. 2013, Section 5.2.1), suggests that for a scale parameter, a more appropriate choice of a prior is one that is inversely proportional to the scale parameter. However, since the use of Jeffreys’ rule for more than one parameters is controversial (Robert et al. 2009, Section 4.7), we adopt two additional priors:

[TABLE]

and

[TABLE]

We note that the lower limit on $\sigma_{\rm tot}$ is set to be a small but non-zero value ($1\,\mathrm{km\ s}^{-1}$) so that the prior distribution can be normalized. The upper limit on $\sigma_{\rm tot}$ is chosen to be $v_{\rm esc}$ so that most of the stars are bound to the Milky Way, but we have confirmed that our results do not change when a larger value is adopted.

In practice, it is fair to state that these three priors are equally uninformative, so the use of any one of them is equally justified. It is important to note that priors A and B are independent of $(\sigma_{\theta},\sigma_{\phi})$, which implies that using these prior distributions is equivalent to setting a flat prior on velocity anisotropy $\beta$. In contrast, prior C is weighted heavily towards small values of $(\sigma_{\theta},\sigma_{\phi})$, and therefore towards large value of $\beta(\simeq 1)$.

3. Mock catalogues

Here we describe how we generate the mock catalogues with which we test the maximum-likelihood and Bayesian methods.

3.1. Assumptions on our mock catalogues

In generating the mock catalogues we first assume that the sample stars obey the distribution function model in equation (2) with no net rotation $(V_{\rm rot}=0)$. Since we want to quantify the error associated with $\beta$, we simply assume that $(\sigma_{r},\sigma_{\theta},\sigma_{\phi})=(1,\sqrt{1-\beta_{\rm true}},\sqrt{1-\beta_{\rm true}})\times(100\;{\rm{km\;s^{-1}}})$, independent of $r$.

Also, we assume that each mock catalogue contains 1000 stars, and all of them have an identical Galactocentric radius. In reality, most of the previous studies used a few thousand halo stars in total, with stars binned according to their Galactocentric radii. Since such a bin typically contains a few hundred stars, our mock catalogues are better populated than reality. Also, our mock catalogues are much simpler to analyze since we can ignore the radial dependence of the halo density.

Furthermore, we assume a simple model for the spatial selection function that mimics the Sloan Digital Sky Survey (SDSS). Specifically, all the stars in a given mock catalogue are distributed at high Galactic latitude with $|b|>30^{\circ}$ and are distributed more than $5\,\mathrm{kpc}$ away from the Galactic disc plane, but otherwise they are distributed uniformly in $(\sin\theta,\phi)$-space. We have confirmed that our results are essentially unchanged if we do not apply the cuts on Galactic latitude or distance from the disc plane.

Lastly, once we generate mock catalogues, we transform the 3-dimensional velocities in the Galactocentric frame to a line-of-sight velocity in the frame of an observer (Sun) moving on a circular orbit with radius $8\,\mathrm{kpc}$ at a velocity of $220\,\mathrm{km\ s}^{-1}$.

These assumptions (especially the assumption of no observational errors) are rather simplistic and idealistic. However, by using mock catalogs with these assumptions, we can be sure that any systematic errors associated with our mock analyses are no less serious than the systematic errors affecting previous local fitting analyses (Sirko et al., 2004; Kafle et al., 2012, 2013; Hattori et al., 2013; King et al., 2015).

3.2. Parameters of our mock catalogues

In this paper, we generate 1000 mock catalogues for a given set of parameters $(r,\beta_{\rm true})$. The Galactocentric radius of the sample stars are assumed to be either $r=6,7,\cdots,30\,\mathrm{kpc}$ in steps of $1\,\mathrm{kpc}$. Also, we adopt eight values of $\beta_{\rm true}=0.75,0.5,\cdots,-1$ in steps of $0.25$. Thus we generate in total $25\times 8\times 1000=2\times 10^{5}$ mock catalogues, each contains 1000 mock stars.

4. Result 1: maximum-likelihood method

Here we investigate the reliability of the maximum-likelihood method. In this Section we analyze our mock catalogues with the maximum-likelihood method to derive $(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ from 4-dimensional information. Then we calculate the corresponding values of velocity anisotropy $\beta_{\rm MaxL}$ (the maximum-likelihood solution for the velocity anisotropy). These calculations were performed by using GNU Scientific Library (Galassi et al, 2009).

4.1. Illustrative results

Since the total number of mock catalogues employed in this paper is huge, we begin by presenting results that illustrate the performance of the maximum-likelihood method, focusing first on analyses of mock catalogues with $\beta_{\rm true}=0.5$ and $-1$.

Figure 1 shows the distributions of the maximum-likelihood solutions for 1000 mock catalogues with $\beta_{\rm true}=0.5$. From top to bottom the panels show the results for $r/\,\mathrm{kpc}=10,15,20$, and $25$, demonstrating how the results of the maximum-likelihood estimation deteriorate with increasing $r$.

The histograms of $\sigma_{r}$ (5th column from the left), show that $\sigma_{r}$ is well estimated at all radii, and the median value of $\sigma_{r}$ (black dashed line) almost coincides with the true value (red solid line). On the other hand, the maximum-likelihood estimates of both tangential velocity dispersion components become broader as $r$ increases and the median and true values increasingly diverge. For example, the histograms of $\sigma_{\phi}$ (right-most column) at $r=10\,\mathrm{kpc}$ and $r=15\,\mathrm{kpc}$ show that the median values of $\sigma_{\phi}$ coincide with the true values. Note that for $r=15\,\mathrm{kpc}$ a fraction of solutions are clustered at $\sigma_{\phi}=0$, and are unrealistic (the origin of these unrealistic solutions is discussed in the Appendix B.) The fraction of unrealistic solutions increases as $r$ increases. Also, we note that the median value of $\sigma_{\phi}$ begins to deviate from the exact value at $r>15\,\mathrm{kpc}$. These properties are also true for $\sigma_{\theta}$. Since $\sigma_{r}$ is well estimated, the error in $\beta_{\rm MaxL}$ is dominated by the errors in $(\sigma_{\theta},\sigma_{\phi})$. As a result, the median value of $\beta_{\rm MaxL}$ begins to deviate from the true value of $\beta_{\rm true}=0.5$ at $r>15\,\mathrm{kpc}$ (see 4th column).

The performance of the maximum-likelihood method can be well summarized in the $(\sigma_{\theta},\sigma_{\phi})$-space (3rd column). In this space, a curve of constant $\beta$ (when $\sigma_{r}$ is fixed to the true value) is described by an arc defined by

[TABLE]

In the 3rd column of Figure 1, the arc (black dashed curve) corresponds to $\beta={\rm median}(\beta_{\rm MaxL})$, while the blue dots show the distribution of the solutions. At $r=10\,\mathrm{kpc}$, the maximum-likelihood solutions are distributed compactly around the true values (marked by red lines), and the arc of $\beta(\sigma_{\theta},\sigma_{\phi})={\rm median}(\beta_{\rm MaxL})$ goes through the true location of $(\sigma_{\theta,{\rm true}},\sigma_{\phi,{\rm true}})$. At $r=15\,\mathrm{kpc}$, the distribution of $(\sigma_{\theta},\sigma_{\phi})$ is broadened and a fraction of solutions are found to be unrealistic ($\sigma_{\theta}=0$ or $\sigma_{\phi}=0$). At $r>15\,\mathrm{kpc}$, the distribution of the solutions is broadened further and the fraction of unrealistic solution is increased. Also, some fraction of solutions attain a large value of $\sigma_{\theta}^{2}+\sigma_{\phi}^{2}$, which corresponds to a highly negative value of $\beta_{\rm MaxL}$. As a result, the arc of median $\beta_{\rm MaxL}$ begins to deviate from the true location of $(\sigma_{\theta,{\rm true}},\sigma_{\phi,{\rm true}})$ at $r>15\,\mathrm{kpc}$.

Figure 2 shows the same results but with mock catalogues with $\beta_{\rm true}=-1$. Again, the maximum-likelihood solutions deteriorate at $r>15\,\mathrm{kpc}$. In this case, the median value of $\beta_{\rm MaxL}$ happens to stay very close to the exact value of $\beta_{\rm true}=-1$ even at $r=25\,\mathrm{kpc}$. However, this result only suggests that, in the case that $\beta_{\rm true}=-1$, the maximum-likelihood method on average returns the unbiased value of $\beta_{\rm MaxL}$ for a large number of independent datasets.

As we can see from the highly broadened histogram of $\beta_{\rm MaxL}$ (4th column), the maximum-likelihood method hardly ever constrains the true anisotropy at $r>15\,\mathrm{kpc}$ if we only use $\sim 10^{3}$ stars for a given radius. In principle, the quality of the estimate of $\beta$ can be improved by increasing the sample to $\sim 10^{4}$ stars at a given radius (see Appendix C and Section 6.1). However, the prospects for obtaining line-of-sight velocities for such a large sample is observationally infeasible in the near future.

4.2. Detailed analyses of maximum-likelihood method

In Section 4.1, we demonstrated that the maximum-likelihood method becomes unreliable at large $r$, especially at $r>15\,\mathrm{kpc}$. Here we have a closer look at this problem.

For each pair of $(r,\beta_{\rm true})$, we have 1000 mock catalogues, so we have an ensemble of 1000 solutions. In order to evaluate the statistical properties of the solutions, we derive the 2.5, 16, 50, 84 and 97.5 percentiles of $\beta_{\rm MaxL}$ and investigate how the percentile ranges depend on $r$ and $\beta_{\rm true}$.

4.2.1 Results for fixed $\beta_{\rm true}$

Here we investigate how the performance of the maximum-likelihood method depends on $r$.

Figure 3 shows the distribution of $\beta_{\rm MaxL}$ as a function of $r$. In each panel, the value of $\beta_{\rm true}$ is fixed, and it is shown by the horizontal red line. In order to understand the systematic error on $\beta_{\rm MaxL}$, let us first focus on the behavior of the median value of $\beta_{\rm MaxL}$. From Figure 3, we can see that the median value of $\beta_{\rm MaxL}$ matches the value of $\beta_{\rm true}$ within a certain radius $r_{\rm reliable}$ (marked by the vertical black line). For example, when $\beta_{\rm true}=0.5$, the median value of $\beta_{\rm MaxL}$ coincides with $\beta_{\rm true}$ at $r<r_{\rm reliable}=15\,\mathrm{kpc}$. At $r>r_{\rm reliable}$, the median value of $\beta_{\rm MaxL}$ is systematically smaller than $\beta_{\rm true}$, and this deviation from $\beta_{\rm true}$ grows larger with increasing $r$. These properties are also true for other values of $\beta_{\rm true}$. Although the value of $r_{\rm reliable}$ is slightly smaller than $15\,\mathrm{kpc}$ for $\beta_{\rm true}=0.75$, and larger for $\beta_{\rm true}=0$ and $-1$, it is safe to say that the maximum-likelihood solutions are reliable within $r\leq 15\,\mathrm{kpc}$ given the small offset between $\beta_{\rm true}$ and the median value of $\beta_{\rm MaxL}$.

We now focus on the spread in the distributions of $\beta_{\rm MaxL}$. The blue solid curves in each panel of Figure 3 show the 16 and 84 percentiles of the distributions of $\beta_{\rm MaxL}$. Similarly, the blue dashed curves in this figure show the 2.5 and 97.5 percentiles. We hereafter refer to these percentile ranges bracketing 68% and 95 % of the distribution as one- and two-$\sigma$ ranges, respectively. We can see from these panels that the spread of these distributions grows rapidly as a function of $r$ and that its growth is especially prominent at $r\gtrsim 15\,\mathrm{kpc}$. This finding can be understood in a following manner. As $r$ becomes larger, the line-of-sight direction \mbox{\boldmathe}_{\rm los} becomes closer to the radial direction \mbox{\boldmathe}_{r}. This means that the line-of-sight velocity $v_{\rm los}$ is more dominated by the radial velocity $v_{r}$. Therefore, when $r$ is large enough compared to the Galactocentric radius of the Sun ($8\,\mathrm{kpc}$), the contribution of $\sigma_{\theta}$ or $\sigma_{\phi}$ to $\sigma_{\rm los}$ becomes less significant, making it harder to reliably extract information on the tangential velocity components. A geometrical explanation for this result is given in Appendix A.

To summarize, the maximum-likelihood method tends to underestimate the value of $\beta$ beyond a certain radius $r_{\rm reliable}$, and that this systematic error on $\beta$ increases as the Galactocentric radius $r$ of the sample increases, and as the true anisotropy $\beta_{\rm true}$ becomes larger. Also, the random error on $\beta$ increases with increasing $r$ at $r\gtrsim 15\,\mathrm{kpc}$. These systematic and random errors make the estimated values of $\beta$ unreliable at large $r$. Based on these findings, it is fair to conclude that the maximum-likelihood method can in principle reliably estimate $\beta$ at $r\leq 15\,\mathrm{kpc}$, but is unable to estimate $\beta$ at $r>15\,\mathrm{kpc}$ if one uses only 4-dimensional information for $1000$ stars at a given radius.

4.2.2 Results for fixed $r$

Here we shall view our results from a different perspective and investigate how the performance of the maximum-likelihood method depends on $\beta_{\rm true}$.

Figure 4 shows the distribution of $\beta_{\rm MaxL}$ as a function of the input anisotropy $\beta_{\rm true}$. In each panel, the Galactocentric radius $r$ is fixed, and the diagonal red line indicates the line of $\beta_{\rm MaxL}=\beta_{\rm true}$. From Figure 4, we can see the growth of the systematic and random errors as a function of $r$. At $r=10\,\mathrm{kpc}$, we see that the maximum-likelihood method on average returns the correct velocity anisotropy, independent of $\beta_{\rm true}$. At $r=15\,\mathrm{kpc}$, the one-$\sigma$ (68%) range of $\beta_{\rm MaxL}$ becomes wider (larger random error), but the median value of $\beta_{\rm MaxL}$ is still very close to $\beta_{\rm true}$. However, at $r>15\,\mathrm{kpc}$, the systematic error on $\beta_{\rm MaxL}$ becomes prominent. Especially, if $\beta_{\rm true}>0$, the median value of $\beta_{\rm MaxL}$ is systematically smaller than $\beta_{\rm true}$. This systematic offset as well as the larger random error on $\beta_{\rm MaxL}$ indicates that there is a large probability that the maximum-likelihood method mistakenly returns a highly negative value of $\beta_{\rm MaxL}$ at $r>15\,\mathrm{kpc}$ even if the true value of $\beta_{\rm true}$ is positive. Although the median value of $\beta_{\rm MaxL}$ is much closer to $\beta_{\rm true}$ if $\beta_{\rm true}<0$, the one-$\sigma$ range is so large at $r>15\,\mathrm{kpc}$ that it is practically impossible to determine if a measured negative $\beta_{\rm MaxL}$ results from a positive or a negative value of $\beta_{\rm true}$.

5. Result 2: Bayesian method

In Section 4, we have found that the performance of the maximum-likelihood method deteriorates at $r>15\,\mathrm{kpc}$. In this Section, we shall confirm this result by using a Bayesian method.

Due to the relatively large computational cost of Bayesian analyses, we apply this method only to a fraction of our mock catalogues. By using 100 mock catalogs for each pair of $(r,\beta_{\rm true})$, we derive the posterior distributions of $(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ as well as $\beta_{\rm Bayes}$ (hereafter $\beta_{\rm Bayes}$ denotes the velocity anisotropy obtained from Bayesian analyses). We used three types of priors, A, B, and C (see Section 2.2), but it turned out that the use of priors A and B results in almost identical posterior distributions. Since our main aim here is to quantify the systematic error in this method, we combine these 100 posterior distributions for each pair of $(r,\beta_{\rm true})$ for each prior. Then we calculate the 2.5, 16, 50, 84, and 97.5 percentiles of the posterior distribution of $\beta_{\rm Bayes}$. These calculations were performed by using a publicly available Python package emcee  (Foreman-Mackey et al., 2013).

5.1. Results for fixed $\beta_{\rm true}$

Figure 5 shows the posterior distribution of $\beta_{\rm Bayes}$ as a function of $r$. In each panel, the value of $\beta_{\rm true}$ is fixed and shown by the red horizontal line. The black solid line shows the median value of the posterior distribution of $\beta_{\rm MaxL}$, and blue solid and blue dashed lines respectively cover 68% and 95% of the posterior distribution.

The results of the Bayesian analyses for priors A and B are almost identical to each other, while the results for priors B and C look distinctly different (hence we do not show results for prior A in Figure 5). This prior-dependence of the resultant posterior distributions can be explained in a following manner. On the one hand, the uncertainty in $\sigma_{r}$ is relatively small even at large Galactocentric radius $r$ (see Section 4.1). This small uncertainty means that the likelihood function is strongly peaked near the true value of $\sigma_{r}$. Thus, although the priors A and B have very different $\sigma_{r}$-dependence, the resultant posterior distributions of $\sigma_{r}$ are almost the same. On the other hand, errors in $\sigma_{\theta}$ and $\sigma_{\phi}$ are large at large $r$ (see Section 4.1), and the error in $\beta_{\rm Bayes}$ is dominated by such errors. This large uncertainty means that the likelihood function only weakly depends on $(\sigma_{\theta},\sigma_{\phi})$, so the posterior distribution is sensitive to the $(\sigma_{\theta},\sigma_{\phi})$-dependence of the prior. Since prior C is weighted more heavily towards smaller values of $(\sigma_{\theta},\sigma_{\phi})$, adopting prior C is equivalent to adopting a strong prior on $\beta$. As a result, the resultant posterior distribution is strongly peaked near $\beta_{\rm Bayes}\simeq 1$, making the median value of $\beta_{\rm Bayes}$ biased toward large values. Since prior B does not depend on $(\sigma_{\theta},\sigma_{\phi})$, adopting prior B is equivalent to adopting a flat prior on $\beta$. As a result, the posterior distribution traces the $(\sigma_{\theta},\sigma_{\phi})$-dependence of the likelihood function, making the posterior distribution of $\beta_{\rm Bayes}$ for prior B more or less similar to the distribution of $\beta_{\rm MaxL}$.

Next we compare the posterior distribution $\beta_{\rm Bayes}$ for prior B with the distribution of $\beta_{\rm MaxL}$. As seen in Figure 5, the median value of $\beta_{\rm Bayes}$ is close to $\beta_{\rm true}$ at small Galactocentric radii, but it gradually deviates from $\beta_{\rm true}$ at large $r$. When $\beta_{\rm true}>0$ and $r\gtrsim 15\,\mathrm{kpc}$, the median value of $\beta_{\rm Bayes}$ is systematically lower than $\beta_{\rm true}$, similar to the results of maximum-likelihood method. However, when $\beta_{\rm true}<0$ and $r\gtrsim 15\,\mathrm{kpc}$, the median value of $\beta_{\rm Bayes}$ is larger than $\beta_{\rm true}$, unlike in the case of the maximum-likelihood method. In any case, the one- and two-$\sigma$ ranges in the $\beta_{\rm Bayes}$-distribution rapidly increases at $r\gtrsim 15\,\mathrm{kpc}$, making the estimation of velocity anisotropy as difficult and uncertain as with the maximum-likelihood method.

5.2. Results for fixed $r$

Figure 6 is similar to Figure 4 and shows the distribution of $\beta_{\rm Bayes}$ as a function of $\beta_{\rm true}$. In each panel, the Galactocentric radius $r$ of the sample is fixed, and the diagonal red line indicates the line of $\beta_{\rm Bayes}=\beta_{\rm true}$.

From the top row of this figure, we see that the Bayesian analysis on average returns a nearly unbiased estimate of velocity anisotropy independent of $\beta_{\rm true}$, for $r\leq 15\,\mathrm{kpc}$ and for prior B (nearly identical plots were obtained for prior A and are not shown). However, the spread in the posterior distribution becomes quite large at $r>15\,\mathrm{kpc}$, making it hard to infer the true anisotropy. On the other hand, the bottom row of Figure 6 (as well as the bottom row of Figure 5) shows that a prior that is biased towards large $\beta$ like prior C results in an overestimate of $\beta$ $r>10\,\mathrm{kpc}$.

6. Discussion

In Sections 4 and 5, we have explored the systematic and random errors inherent in the local fitting methods with maximum-likelihood and Bayesian formulations. We found that these methods can reliably estimate $\beta$ only at $r\leq 15\,\mathrm{kpc}$. In order to better understand the errors inherent in these methods, we discuss how the performance improves if we increase the sample size in Section 6.1. Then we discuss how our results in this paper can be used in interpreting the measured values of $\beta$ in Section 6.2. We also comment on the use of different (non-Gaussian) functions in the local fitting methods in Section 6.3.

6.1. Effects of the increased sample size

From a mathematical point of view, we expect that we can recover the value of $\beta$ from mock data accurately if we have a large enough sample of stars (and we know the functional form of the distribution function). In order to confirm this expectation, we did the same analyses as in Sections 4.1 and 4.2.2 with sample size of $N=10^{4}$ (at a given radius) instead of $N=10^{3}$. We found that we can recover $\beta$ accurately out to a significantly larger radius of $r=25\,\mathrm{kpc}$ (see Appendix C for details). From these experiments, we conclude that the systematic bias in $\beta_{\rm MaxL}$ seen in Figures 1-4 arises from the small sample size of $N=10^{3}$. Also, from a geometric argument in Appendix A, we found that the random error on $\beta_{\rm{MaxL}}$ at $r>R_{0}$ ($R_{0}$ is the Galactocentric radius of the Sun) is approximately given by

[TABLE]

This expression indicates that the performance of the maximum-likelihood method deteriorates as $r$ increases and improves as $N$ increases. We expect the similar results would be obtained if we use Bayesian method and prior B, based on our results in Sections 4 and 5.

We warn that these experiments do not guarantee that a reliable estimation of $\beta$ can always be obtained with with just 4-dimensional information for an adequately large number of sample stars. Local fitting measurements of $\beta$ presented in this paper make use of the fact that $\sigma_{\rm los}$ depends on the heliocentric direction on the sky when the velocity distribution is anisotropic ($\beta\neq 0$) but not if it is isotropic ($\beta=0$) [see equation (4)]. Our analyses use this direction-dependence of $\sigma_{\rm los}$ across the sky to estimate $\beta$ by assuming that the principal axes of the velocity ellipsoid are perfectly aligned with the Galactocentric spherical coordinates and that $(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ are functions of $r$ only. In the Solar neighborhood, these assumptions are approximately valid (Bond et al., 2010; Evans et al., 2016), but there is no guarantee that they are valid outside the Solar neighborhood. Therefore, we expect that in order to obtain a reliable determination of $\beta$ at $r\gtrsim 20\,\mathrm{kpc}$ proper motion data for halo stars are required.

6.2. Literature values of $\beta$ based on the local fitting methods

In this paper, we found that the use of 4-dimensional information for $\sim 1000$ halo stars at a given radius is not sufficient to reliably estimate $\beta$ for Galactocentric radii $r>15\,\mathrm{kpc}$ (even with error-free, noise-free data). In this subsection we use our results to better understand recent observational determinations of $\beta$ that were based on similar local fitting methods. In particular, we focus on Kafle et al. (2012) and King et al. (2015) as the typical studies in which the Bayesian and the maximum-likelihood methods are applied to a halo sample without any metallicity cuts. Also, we consider the results from Hattori et al. (2013), who used a slightly different method based on solving the matrix equation to 4-dimensional information for a halo sample that was split into two different metallicity ranges. (We note Kafle et al. 2013 did essentially the same analyses as in Hattori et al. 2013, but Kafle et al. 2013 did not discard data points at large $r$ which, as we have shown, are likely to be biased.)

In Figure 7, we show the point-estimate distribution of $\beta_{\rm MaxL}$ taken from Figure 3 as well as the measurements of $\beta$ in the above-mentioned papers. It is worth noting that all published measurements of $\beta$ from line-of-sight velocities rely on the maximum-likelihood, Bayesian, or other similar local fitting methods and that the results of these methods are quite similar (as shown in this paper). Hence a comparison of $\beta$ profiles obtained from maximum-likelihood analyses of our mock data with observationally determined $\beta$ profiles provides a useful way to evaluate the accuracy of previous observational measurements. The distribution of $\beta_{\rm MaxL}$ in Figure 7 represents the expected distribution of the maximum-likelihood solution for a given value of $\beta_{\rm true}$. Thus, if a given data point (not the error bar) in Figure 7 is located outside the two-$\sigma$ range of the distribution of $\beta_{\rm MaxL}$ for a certain value of $\beta_{\rm true}$, then that value of $\beta_{\rm true}$ is disfavored by the data point.

First, we focus on the panels with $\beta_{\rm true}=0$ and $0.5$. These panels suggest that even if the Milky Way stellar halo has a constant profile of, say, $\beta(r)=0.5$, most of the measured values of $\beta$ from observations (data points) shown here lie within two-$\sigma$ of the expected deviation (using the maximum-likelihood method).

Second, let us focus on data points of Kafle et al. (2012) and King et al. (2015) at $r>13\,\mathrm{kpc}$. For these data points, we see a declining profile of $\beta(r)$ as a function of $r$, although the error bars are quite large at large radii. This declining profile is a reminiscent of the declining profile of the median curve of $\beta_{\rm MaxL}$. Therefore, even if the measured $\beta(r)$ profile is mildly declining, it does not necessarily mean that the true $\beta(r)$ is declining. The apparent dip of $\beta(r)$ at $r\simeq$(15-17)$\,\mathrm{kpc}$ is intriguing, and this dip might be a true signal. Indeed, Loebman et al. (2016) demonstrate that dips in $\beta(r)$ can arise due to substructure in the stellar halo. However, it is worthwhile to point out that the measurement of a dip in $\beta$ profile at $r=17\,\mathrm{kpc}$ reported by Kafle et al. (2012) is consistent with $-1\leq\beta_{\rm true}\leq 0.5$. Also, the data points at $r\gtrsim 20\,\mathrm{kpc}$ do not seem to help our understanding of $\beta(r)$ profile, since both the systematic and random errors grows rapidly as a function of $r$. For example, all the data points at $r>20\,\mathrm{kpc}$ of Kafle et al. (2012) and King et al. (2015) are consistent with $-1\leq\beta_{\rm true}\leq 0.75$, which is the full range of $\beta_{\rm true}$ we have explored in this paper.

Third, let us focus on data points of Kafle et al. (2012) and King et al. (2015) at $r\simeq 12\,\mathrm{kpc}$. Although their sample stars are partially overlapping (both of their samples include SDSS blue-horizontal branch stars), their estimated values of $\beta$ are inconsistent with each other if their error bars are correct. Currently we do not know the origin of this discrepancy, but the published error bars might be too small. For example, most of their data points at $r\simeq 12\,\mathrm{kpc}$ are consistent with our results of $\beta_{\rm true}=0.5$ within one-$\sigma$ range. On the other hand, the error bars in Kafle et al. (2012) and King et al. (2015) are $\sim$30% of the one-$\sigma$ range of our mock catalogue analyses, while the error bars in Hattori et al. (2013) are $\sim$ 65% of the one-$\sigma$ range in this paper. Since our mock catalogues do not include observational errors and yet the error bars in our mock catalogues are larger than the above-mentioned papers (Kafle et al., 2012; Hattori et al., 2013; King et al., 2015), the published error bars might not represent the uncertainty in $\beta$. However, it is premature to conclude that the published error bars are incorrect, since our mock catalogues are not realistic enough (e.g., we do not consider the spread in $r$ of sample stars).

Lastly, let us now focus on data points of Hattori et al. (2013). Based on the apparent difference of $\beta$ for the metal-poor and metal-rich samples, they claimed that the kinematics of the halo stars depends on the metallicity (for stars with $r<18\,\mathrm{kpc}$). Although they carefully avoid possible systematic errors on $\beta$ by discarding the results at $r>16\,\mathrm{kpc}$ for metal-rich halo stars and at $r>18\,\mathrm{kpc}$ for metal-poor halo stars, their cut may not be sufficient. Based on the analyses in this paper, we argue that a safer approach is to discard data points at $r>15\,\mathrm{kpc}$. (Even with this spatial cut, a metallicity-dependence in $\beta$ can still be seen in the surviving data points.) Admittedly, there is a possibility that these differences appeared by mere chance. For example, if the Milky Way stellar halo has a constant profile of $\beta(r)=\beta_{\rm true}=0$, then the estimated $\beta(r)$ profiles of both the metal-rich and metal-poor samples are inside the one-$\sigma$ range expected from our model. Although the current 4-dimensional data are not good enough to conclusively assert that $\beta(r)$ depends on stellar metallicity, the current data within $15\,\mathrm{kpc}$ do hint at such a possibility. For example, let us suppose that the $\beta$ profile of metal-rich halo is $\beta(r)=\beta_{\rm true}=0.5$. In this case, the data points for the metal-poor sample at $r\leq 15\,\mathrm{kpc}$ are marginally outside the one-$\sigma$ range of the model prediction. Conversely, let us suppose that the metal-poor halo has $\beta(r)=\beta_{\rm true}=-1$. In this case, most of the data points of the metal-rich sample at $r\leq 15\,\mathrm{kpc}$ are outside the one-$\sigma$ range. A robust determination of a metallicity dependence in $\beta(r)$ awaits confirmation with kinematical data from Gaia and chemical information from ground-based surveys.

6.3. Local fitting methods with non-Gaussian functions

In previous sections, we generated mock catalogues based on a Gaussian velocity distribution and performed local fitting analyses by assuming that the underlying velocity distribution is also a Gaussian function. However, in reality we do not know the correct functional form of the velocity distribution of the halo stars. In order to investigate the reliability of our local fitting methods, we perform some additional tests.

First, we introduce two simple distribution function models that are both functions of energy $E$ and total angular momentum $L$ (see Appendix D for details). One model, $f_{\rm const}(E,L)$, has a constant profile of $\beta(r)=\beta_{\rm const}$, and another model, $f_{\rm OM}(E,L)$, is an Osipkov-Merritt model with a rising $\beta$ profile of $\beta(r)=r^{2}/(r_{a}^{2}+r^{2})$ with $r_{a}$ a constant. Then we assume that the potential of the Milky Way is a spherical singular isothermal potential (with a flat rotation curve) and generate three sets of mock catalogues in the same manner as in Section 3. For two sets of mock catalogues, we adopt $f_{\rm const}(E,L)$ with $\beta_{\rm const}=0.25$ and $-0.42$. For the other set of mock catalogues, we adopt $f_{\rm OM}(E,L)$ with $r_{a}=10\,\mathrm{kpc}$.

We fit each of these three sets of mock catalogues with local fitting methods with the maximum-likelihood formulation. Specifically, at each Galactocentric radius $r$, we fit the data by assuming that the underlying velocity distribution is described by either f_{\rm Gauss}(\mbox{\boldmathv}|\mbox{\boldmathx}), f_{\rm const}(E,L|\mbox{\boldmathx}), or f_{\rm OM}(E,L|\mbox{\boldmathx}). Here, f(E,L|\mbox{\boldmathx}) denotes the velocity distribution at a given location $x$ of a system obeying a distribution function $f(E,L)$, and it is different from $f(E,L)$ itself. In the following, however, we omit the arguments of $f_{\rm const}$ and $f_{\rm OM}$ for brevity.

Figure 8 shows the results of these analyses. As seen in Figure 8(f), when the mock catalogues generated from $f_{\rm OM}$ are locally fitted with $f_{\rm OM}$, the median value of the $\beta_{\rm MaxL}$ for the mock catalogues almost overlaps the true profile of $\beta(r)$ at $6\,\mathrm{kpc}\leq r\leq 30\,\mathrm{kpc}$. However, at $r>15\,\mathrm{kpc}$, the one-$\sigma$ range fills essentially the entire range of $0\leq\beta\leq 1$, which is the entire range of the allowed anisotropy in Osipkov-Merritt models (recall that the Osipkov-Merritt model does not permit models with negative $\beta$). Thus, the recovered value of $\beta$ is informative only at $r\leq 15\,\mathrm{kpc}$, just as in the results in Figure 3. On the other hand, as seen in Figures 8(a) and 8(b), when the mock catalogues generated from $f_{\rm const}$ are locally fitted with $f_{\rm const}$, the median value of the $\beta_{\rm MaxL}$ for the mock catalogues is very close to the correct value $\beta_{\rm const}$ at $6\,\mathrm{kpc}\leq r\leq 30\,\mathrm{kpc}$ for both cases of $\beta_{\rm const}=0.25$ and $-0.42$. Interestingly, the uncertainty in $\beta$ does not change a lot as a function of $r$ for this model. This result suggests that choosing the correct functional form is beneficial in estimating $\beta$.

When an incorrect velocity distribution is assumed, on the other hand, the resultant $\beta$ profiles are sometimes not reliable. For example, as seen in Figure 8(c), when mock catalogues generated from $f_{\rm OM}$ are locally fitted with $f_{\rm const}$, the $\beta(r)$ profile can not be recovered even at $r\leq 10\,\mathrm{kpc}$. What is interesting in these fits is that the formal uncertainties associated with the recovered value of $\beta$ are very small despite the fact that the true profile of $\beta(r)=r^{2}/((10\,\mathrm{kpc})^{2}+r^{2})$ is well outside the two-$\sigma$ range. Also, as seen in Figure 8(e), when mock catalogues generated from $f_{\rm const}$ with $\beta_{\rm const}=-0.42$ are locally fitted with $f_{\rm OM}$, the estimated profile of $\beta$ is far from the true profile at all the radii explored, since $\beta<0$ can not be attained with $f_{\rm OM}$. These examples suggest that when local fitting methods are applied to 4D data, the use of wrong functional forms for the velocity distribution can result in a significant systematic error on the recovered value of $\beta$. On the other hand, as seen in Figures 8(g)-(i), when we use $f_{\rm Gauss}$ to fit the mock catalogues generated from $f_{\rm const}$ and $f_{\rm OM}$, the $\beta$ profile is reliably estimated at $r<15\,\mathrm{kpc}$. This result is intriguing, since the velocity distribution for $f_{\rm const}$ is non-Gaussian.444 We note that the velocity distribution for $f_{\rm OM}$ is always Gaussian in our current example.

These results, combined with the findings in Section 4, can be summarized in the following manner. When local fitting methods are used to estimate $\beta$,

• •

if the correct functional form of the velocity distribution is assumed, the estimated value of $\beta$ and the associated error bars are reliable at least at $r<15\,\mathrm{kpc}$, and even at $r=30\,\mathrm{kpc}$ if the halo obeys $f_{\rm const}$ [see Figures 8(a),(b),(f)];

• •

if a wrong functional form of the velocity distribution is assumed, the estimated value of $\beta$ may be significantly biased (even at $r<10\,\mathrm{kpc}$), and the associated formal error may not be reliable [see Figures 8(c),(e)];

• •

if the functional form of the velocity distribution is unknown, the assumption of a Gaussian function $f_{\rm Gauss}$ is a reasonable choice, since this choice allows unbiased estimation of $\beta$ at $r\leq 15\,\mathrm{kpc}$ for the various kinds of mock catalogues explored in this paper [see Figures 8(g)-(i)].

Recently, many authors have developed sophisticated distribution function models that fit the observed positions and velocities of halo stars in the Milky Way (Deason et al., 2012; Williams & Evans, 2015; Das & Binney, 2016; Das et al., 2016). With these global fitting methods, authors use a sample of halo stars that are not restricted to a single Galactocentric radius $r$ (as with the local fitting methods), but are distributed over a wide range of $r$. However, given that our local fitting methods fail to recover $\beta$ (sometimes even at $r<10\,\mathrm{kpc}$) when a wrong functional form is assumed, it may be worthwhile checking the performance of these global fitting methods when a wrong functional form of the stellar halo distribution function is assumed.

7. Conclusion

In the past 10 years, many authors have tried to infer the velocity distribution of distant halo stars from stellar samples without reliable proper motion measurements (see references in Section 1). A common way of inferring the 3-dimensional velocity dispersion of halo stars from 4-dimensional position and line-of-sight velocity measurements is the local fitting methods. In these methods, they estimate 3-dimensional velocity dispersion by using information on how the line-of-sight velocity dispersion $\sigma_{\rm los}$ varies across the sky, which reflects the different line-of-sight projections of the velocity ellipsoid. However, as stars get farther away from us and from the Galactic center, $v_{\rm los}$ (corrected for the Solar reflex motion) becomes increasingly closer to $v_{r}$. As a result, $\sigma_{\rm los}$ becomes closer to $\sigma_{r}$, and the variation of $\sigma_{\rm los}$ across the sky becomes harder to evaluate, making it difficult to estimate the tangential components of the velocity dispersion (see Appendix A). Thus it is important to explore the random and systematic uncertainties inherent in such methods using mock datasets, in order to build intuition about how far out in the stellar halo we can reliably recover the velocity anisotropy. In this paper, we tackled this problem by performing a series of mock analyses. The main messages of this paper can be summarized as follows.

As shown in Figure 4, the local fitting methods with the maximum-likelihood formulation can in principle reliably estimate the velocity anisotropy $\beta$ of the stellar halo at $r\leq 15\,\mathrm{kpc}$ but is unable to reliably estimate $\beta$ at larger $r$ if 4-dimensional data only (position and line-of-sight velocity) are used for $N=10^{3}$ halo stars at a given radius. 2. 2.

If a flat prior on $\beta$ is adopted, the Bayesian formulation yields similar results to those obtained with the maximum-likelihood formulation [see Figure 4 and the top row (prior B) of Figure 6]. 3. 3.

Previous local fitting analyses to measure $\beta$ from 4-dimensional information used a few hundred halo stars (in each radial bin). Our results suggest that these measurements of $\beta$ are very likely to be biased to low/negative values at $r>15\,\mathrm{kpc}$. 4. 4.

The performance of the local fitting methods to estimate $\beta$ with 4-dimensional information can be improved if we enlarge the sample size (see Section 6.1 and Appendix C). However, we expect that a more direct way of better estimating $\beta$ at large $r$ is to use accurate proper motion data. 5. 5.

If the correct functional form of the velocity distribution of the stellar halo is assumed, the performance of the local fitting methods to estimate $\beta$ can be improved; otherwise assuming a Gaussian velocity distribution is a reasonable choice (see Section 6.3 and Figure 8). 6. 6.

It is important to point out that Deason et al. (2013) and Cunningham et al. (2016) used halo stars at $18\,\mathrm{kpc}<r<30\,\mathrm{kpc}$ with reliable proper motion data and reported $\beta\simeq 0$, although their sample size was small ($N=13$). Their results combined with the results in this paper suggest that the negative values of $\beta$ at $r>20\,\mathrm{kpc}$ reported by Kafle et al. (2012) and King et al. (2015) likely resulted from the large systematic and random errors inherent in the maximum-likelihood method (see Figure 7). 7. 7.

In this paper we have used idealized mock catalogues for which the stellar heliocentric distances and line-of-sight velocities are measured with infinite precision. Also, all the stars in each mock catalogue are assigned the same Galactocentric radius $r$. The errors and noise in real data will only make the task of measuring $\beta$ from 4-dimensional data even harder. However, in the next 3 to 5 years, Gaia (Perryman et al., 2001; Lindegren et al., 2016) will provide proper motion data for a large number of halo stars, opening new avenues for measuring the velocity distribution of halo stars at $r>20\,\mathrm{kpc}$. These new data will yield important insights into our understanding of the structure and the merger history of the Milky Way.

The authors thank the stellar halo group at the Department of Astronomy, University of Michigan for stimulating discussions. MV and KH are supported by NASA-ATP award NNX15AK79G. SRL acknowledges support from the Michigan Society of Fellows.

Appendix A Geometry of local fitting methods

Here we use a geometrical argument to explain why the estimation of $\beta$ deteriorates as $r$ increases and improves as $N$ increases. We assume that the velocity distribution obeys a Gaussian distribution described in equation (2). Also, for brevity, we assume that the sample stars have an identical Galactocentric radius $r>R_{0}$, where $R_{0}$ is the Galactocentric radius of the Sun, and that they are distributed along Galactic longitude $\ell=0^{\circ}$ or $180^{\circ}$ (with any value of Galactic latitude $b$). In this case, the line-of-sight velocity dispersion $\sigma_{\rm los}$ at a given position in the Milky Way is given by

[TABLE]

where Q_{r}=\mbox{\boldmathe}_{\rm los}\cdot\mbox{\boldmathe}_{r} [see equation (4)]. We note that this linear dependence of $\sigma_{\rm los}^{2}$ on $Q_{r}^{2}$ suggests that $\sigma_{\rm los}^{2}=\sigma_{\theta}^{2}$ at $Q_{r}^{2}=0$, while $\sigma_{\rm los}^{2}=\sigma_{\theta}^{2}$ at $Q_{r}^{2}=1$. Also, we note that at a given Galactocentric radius $r(>R_{0})$, the maximum value of $Q_{r}^{2}$ is $Q^{2}_{r,\rm{max}}(r)=1$ ($b=0^{\circ}$) and the minimum value of $Q_{r}^{2}$ is $Q^{2}_{r,\rm{min}}(r)=1-{R_{0}^{2}}/{r^{2}}$ ($b=\pm 90^{\circ}$).

If we have 4D information of $N$ stars, these data can be transformed into $\{(Q_{r},v_{\rm los})_{i}|i=1,\cdots N\}$. The local fitting method fits the distribution of $(Q_{r},v_{\rm los})$ with a model in which $\sigma_{\rm los}^{2}$ varies linearly as a function of $Q_{r}^{2}$ as described in equation (A1). To put it differently, if we bin the data according to $Q_{r}^{2}$ and derive $\sigma_{\rm los}^{2}$ for each bin, the derived $\sigma_{\rm los}^{2}$ profile is fitted with a line. The values of $\sigma_{\rm los}^{2}$ at $Q_{r}^{2}=1$ and $Q_{r}^{2}=0$ correspond to the best-fit values of $\sigma_{r}^{2}$ and $\sigma_{\theta}^{2}$, respectively. The data points are distributed only at $Q^{2}_{r,{\rm min}}(r)\leq Q_{r}^{2}\leq 1$, so the derivation of $\sigma_{\theta}^{2}$ requires an extrapolation of the linear relationship inferred at $Q^{2}_{r,{\rm min}}(r)\leq Q_{r}^{2}\leq 1$ to $Q_{r}^{2}=0$ [see Figure 9(b)]. At larger $r$, $Q_{r,\rm{min}}^{2}$ becomes increasingly closer to $1$, so that the estimation of the slope ${\rm d}\sigma_{\rm los}^{2}/{\rm d}Q_{r}^{2}$ becomes increasingly more difficult. This difficulty results in large uncertainty in $\sigma_{\theta}^{2}$ at large $r$ [see Figure 9(c)], and thus sometimes the maximum-likelihood routine finds unphysical solutions of $\sigma_{\theta}^{2}=0$ (if a constraint of $\sigma_{\theta}^{2}\geq 0$ is imposed), or even $\sigma_{\theta}^{2}<0$ (if such a constraint is not imposed; see Appendix B). On the other hand, the estimation of $\sigma_{r}^{2}$ is not difficult, since $\sigma_{r}^{2}$ is approximately the observed value of $\sigma_{\rm los}^{2}$ at $Q_{r}^{2}\simeq 1$.

A.1. Ideal distribution of sample stars

Let us consider a case where we have $N(\gg 1)$ sample stars at a Galactocentric radius $r$, and half of them ($N/2$ stars) are observed in the direction of $Q_{r}^{2}=1$ (hereafter ‘QMAX direction’) and the other half of them are observed in the direction of $Q_{r}^{2}=Q_{r,\rm{min}}^{2}(r)$ (hereafter ‘QMIN direction’). This spatial distribution of stars is not realistic, but is ideal for inferring the slope ${\rm d}\sigma_{\rm los}^{2}/{\rm d}Q_{r}^{2}$ with local fitting method.

In this case, the true values of $\sigma_{\rm los}^{2}$ in the QMAX and QMIN directions are given by

[TABLE]

respectively. By using observed values of $\sigma_{\rm los}^{2}(QMIN)$ and $\sigma_{\rm los}^{2}(QMAX)$, the value of $\beta$ can be expressed as

[TABLE]

Since the distribution of $v_{\rm los}$ follows a Gaussian distribution [see equation (3)], the observed values of $\sigma_{\rm los}^{2}$ in the QMAX and QMIN directions are associated with uncertainties of

[TABLE]

[TABLE]

respectively. By using equation (A4) and by assuming that the uncertainties $\Delta\sigma_{\rm los}^{2}(QMAX)$ and $\Delta\sigma_{\rm los}^{2}(QMIN)$ are not correlated, we can express the uncertainty in $\beta$ as follows:

[TABLE]

A.2. Realistic distribution of sample stars

In reality, the sample stars are distributed in a wide area in $(\ell,b)$-space, and are not confined around the QMAX and QMIN directions. Therefore, we need to rescale the value of $|\Delta\beta|_{\mathrm{ideal}}$. From our results in Section 4, we find that

[TABLE]

is a good approximation to the random error on $\beta$. This expression clearly illustrates how the performance of the maximum-likelihood method deteriorates when $r$ increases and improves when $N$ increases. For example, with $N=10^{3}$ and $\beta_{\rm true}=0.5$, we see that $\Delta\beta$ is smaller than $(1-\beta_{\rm true})$ only at $r<14.5\,\mathrm{kpc}$ and hence estimation of $\beta$ is not reliable beyond this radius (as discussed in Sections 4 and 5). However, when $N=10^{4}$ and $\beta_{\rm true}=0.5$ are assumed, $\Delta\beta$ is smaller than $(1-\beta_{\rm true})$ at $r<24.4\,\mathrm{kpc}$.

Appendix B Origin of the unrealistic solutions in maximum-likelihood analyses

In Section 4.1, we showed that the maximum-likelihood method returns unrealistic solutions with $\sigma_{\theta}=0$ or $\sigma_{\phi}=0$ for a fraction of mock catalogues at $r\gtrsim 15\,\mathrm{kpc}$. In Appendix A, we explained with an geometrical argument that it becomes increasingly more difficult to estimate the tangential components of the velocity dispersion at larger $r$, and that unphysical solutions arise due to this difficulty. Here we investigate the origin of the unrealistic solutions from a different perspective, by performing maximum-likelihood analyses for 100 mock catalogues with $\beta_{\rm true}=0.5$ and $r=15\,\mathrm{kpc}$ in two ways.

The first set of analyses are done with the same formulation as in Section 2.1. To be specific, we search for a set of parameters $(\sigma_{r},\sigma_{\theta},\sigma_{\phi},V_{\rm rot})$ that maximize the log-likelihood $\ln L$ [see equation (5)] under the condition that $0\leq\sigma_{k}<\infty$ $(k=r,\theta,\phi)$ and $-\infty<V_{\rm rot}<\infty$. Since $0\leq\sigma_{k}$ is a physical requirement, we refer to these solutions as physical solutions.

The second set of analyses are done with a different parametrization. Here, we define new variables $(s_{r},s_{\theta},s_{\phi})=(\sigma^{2}_{r},\sigma^{2}_{\theta},\sigma^{2}_{\phi})$ and we search for a set of parameters $(s_{r},s_{\theta},s_{\phi},V_{\rm rot})$ that maximize $\ln L$ under the conditions of $-\infty<s_{k}<\infty$ $(k=r,\theta,\phi)$ and $-\infty<V_{\rm rot}<\infty$. Obviously, when $s_{k}<0$ for any of $k=r,\theta,\phi$, there is no physical distribution function given by equation (2). However, since $\ln L$ is a function of $(\sigma^{2}_{r},\sigma^{2}_{\theta},\sigma^{2}_{\phi},V_{\rm rot})$, it is mathematically justified to search for the solutions with negative values of $s_{k}$. Hereafter, we refer to these solutions as mathematical solutions.

Figure 10 shows the distribution of the physical and mathematical solutions in $(\sigma^{2}_{\theta},\sigma^{2}_{\phi})$-space [or equivalently, $(s_{\theta},s_{\phi})$-space]. The blue dots and magenta crosses indicate the physical and mathematical solutions, respectively. We see that the mathematical solutions are more or less distributed around the exact location of $(\sigma^{2}_{\theta,{\rm true}},\sigma^{2}_{\phi,{\rm true}})$ with a rather large scatter. Since $(\sigma^{2}_{\theta,{\rm true}},\sigma^{2}_{\phi,{\rm true}})$ is located within the first quadrant, mathematical solutions for a large fraction of mock catalogues are located within the first quadrant. In such cases, the physical and mathematical solutions are identical. However, for a fraction of mock catalogues, the mathematical solutions are located outside the first quadrant. In such cases, the corresponding physical solutions are located along the axes of $\sigma_{\theta}^{2}=0$ or $\sigma_{\phi}^{2}=0$.

The origin of these unrealistic solutions can be explained in a following manner. At $r\gtrsim 15\,\mathrm{kpc}$, the log-likelihood $\ln L$ depends weakly on $(s_{\theta},s_{\phi})$ (as mentioned in Section 5.1). Therefore, for a fraction of mock catalogues, depending on the spatial distribution or velocity distribution of the sample stars, $\ln L$ happens to attain its maximum at $s_{\theta}<0$ or $s_{\phi}<0$. In such situations, $\ln L$ increases as $s_{\theta}$ or $s_{\phi}$ decreases within the first quadrant, so that the physical solutions are distributed along the axes.

Since the scatter in the mathematical solutions of $s_{\theta}$ and $s_{\phi}$ becomes larger with increasing $r$ (due to the enhanced difficulties in extracting the information regarding tangential velocity components), a larger fraction of mathematical solutions are located outside the first quadrant in $(\sigma^{2}_{\theta},\sigma^{2}_{\phi})$-space. This is why the fraction of unrealistic physical solutions with $\sigma_{\theta}=0$ or $\sigma_{\phi}=0$ increases with increasing $r$, as seen in Figure 1.

Appendix C Experiments with larger sample size

In Sections 4 and 5, we use $N=10^{3}$ stars at a given radius. Here we briefly explore whether the results improve if we use $N=10^{4}$ stars instead.

To this end, we generated 1000 mock catalogues containing $N=10^{4}$ stars for each pair of $(r,\beta_{\rm true})$. We adopted four values of $r/\,\mathrm{kpc}=10,15,20,$ and $25$ and $\beta_{\rm true}=0.75,0.5,0,$ and $-1$.

First, we did exactly the same analyses as in Section 4.1 by using these mock catalogues. Figures 11 and 12 show the distributions of $\beta_{\rm MaxL}$ for $\beta_{\rm true}=0.5$ and $\beta_{\rm true}=-1$, respectively. A comparison of Figures 1 and 11 suggests that the performance of the maximum-likelihood method is better when $N=10^{4}$ than when $N=10^{3}$ at all the radii explored here. From Figure 11, we note that the histogram of $\beta_{\rm MaxL}$ is peaked at around $\beta_{\rm true}=0.5$ at $r\leq 20\,\mathrm{kpc}$, and that the median value of $\beta_{\rm MaxL}$ coincides with $\beta_{\rm true}$ even at $r=25\,\mathrm{kpc}$. Also, Figure 12 suggests that the peak of the histogram of $\beta_{\rm MaxL}$ as well as the median value of $\beta_{\rm MaxL}$ coincide with $\beta_{\rm true}=-1$ at all the radii explored here. The peaked histogram of $\beta_{\rm MaxL}$ at $r=25\,\mathrm{kpc}$ seen in Figure 12 is in contrast to the highly flattened histogram at $r=25\,\mathrm{kpc}$ seen in Figure 2.

Secondly, we did the same analyses as in Section 4.2.2 by using the mock catalogues with $N=10^{4}$ stars. Figure 13 shows the distribution of $\beta_{\rm MaxL}$ as a function of $\beta_{\rm true}$ for different Galactocentric radius $r$ of sample stars. We see that the median value of $\beta_{\rm MaxL}$ almost perfectly coincides with $\beta_{\rm true}$ at $r\leq 25\,\mathrm{kpc}$. Also, we found that the one- and two-$\sigma$ ranges of the posterior distribution of $\beta_{\rm MaxL}$ for the case of $N=10^{4}$ seen in Figure 13 are significantly smaller than the corresponding ranges for the case of $N=10^{3}$ seen in Figure 4.

These results indicate that the maximum-likelihood method can in principle reliably estimate $\beta$ at $r\leq 25\,\mathrm{kpc}$ if we have $N=10^{4}$ stars at a given radius.

Appendix D Two distribution function models

Here describe some details on the two distribution function models used in Section 6.3, $f_{\rm const}(E,L)$ and $f_{\rm OM}(E,L)$, which are functions of energy $E$ and total angular momentum $L$. In the following, we assume that the potential of the Milky Way is spherical and is expressed as $\Phi(r)=v_{0}^{2}\ln(r/r_{0})$ with $(r_{0},v_{0})=(220\,\mathrm{km\ s}^{-1},8\,\mathrm{kpc})$.

D.1. Constant $\beta$ model

The model with const $\beta(r)$ is given by

[TABLE]

Here, $L_{\rm max}(E)=r_{0}v_{0}\exp(2E/v_{0}^{2}-1)$ is the angular momentum of a star with energy $E$ moving on a circular orbit. We note that $f_{\rm const}(E,L)\geq 0$ is always satisfied if $K\leq 1$. The density profile of this distribution function is given by $\rho(r)\propto r^{-\alpha}$ ($\alpha>0$). The velocity anisotropy is governed by $(K,\alpha)$ and expressed as

[TABLE]

The probability density that a star at $x$ characterized by Q_{r}=\mbox{\boldmathe}_{\rm los}\cdot\mbox{\boldmathe}_{r} has a line-of-sight velocity $v_{\rm los}$ is expressed as

[TABLE]

In Section 6.3, the mock catalogues with $\beta_{\rm const}=0.25$ and $-0.42$ are generated by assuming $(K,\alpha)=(0.84,2)$ and $(-3,2)$, respectively. Given the mock data, the local fitting method finds the pair $(K,\alpha)$ that maximizes the likelihood.

D.2. Osipkov-Merritt model

Osipkov-Merritt model is a broad class of distribution functions that only depends on $Q=E+L^{2}/(2r_{a}^{2})$ with $r_{a}$ a constant. Here we adopt a family of functions of the form

[TABLE]

with $\alpha>0$. The density profile of this distribution function is given by

[TABLE]

The velocity anisotropy is given by $\beta(r)=r^{2}/(r_{a}^{2}+r^{2})$. The probability density that a star at $x$ characterized by Q_{r}=\mbox{\boldmathe}_{\rm los}\cdot\mbox{\boldmathe}_{r} has a line-of-sight velocity $v_{\rm los}$ is expressed as

[TABLE]

where the line-of-sight velocity dispersion is given by

[TABLE]

In Section 6.3, the mock catalogues are generated by assuming $(r_{a},\alpha)=(10\,\mathrm{kpc},2)$. Given the mock data, the local fitting method finds the pair $(r_{a},\alpha)$ that maximizes the likelihood.