How to break the density-anisotropy degeneracy in spherical stellar systems

TL;DR

This paper introduces GravSphere, a non-parametric code for recovering density and anisotropy profiles in spherical stellar systems, and evaluates methods to break the density-anisotropy degeneracy using mock data.

Contribution

The paper demonstrates the effectiveness of combining multiple populations, higher order VSPs, and proper motions to better constrain density and anisotropy profiles, advancing dynamical modeling techniques.

Findings

Multiple populations improve density estimates between their half-mass radii.

Higher order VSPs enable measurement of density and broad constraints on anisotropy.

Proper motion data significantly enhance the accuracy of density and anisotropy profiles.

Abstract

We present a new non-parametric Jeans code, GravSphere, that recovers the density $\rho(r)$ and velocity anisotropy $\beta(r)$ of spherical stellar systems, assuming only that they are in a steady-state. Using a large suite of mock data, we confirm that with only line-of-sight velocity data, GravSphere provides a good estimate of the density at the projected stellar half mass radius, $\rho(R_{1/2})$, but is not able to measure $\rho(r)$ or $\beta(r)$, even with 10,000 tracer stars. We then test three popular methods for breaking this $\rho-\beta$ degeneracy: using multiple populations with different $R_{1/2}$; using higher order `Virial Shape Parameters' (VSPs); and including proper motion data. We find that two populations provide an excellent recovery of $\rho(r)$ in-between their respective $R_{1/2}$. However, even with a total of $\sim 7,000$ tracers, we are not able to well-constrain $\beta(r)$ for either population. By contrast, using 1000 tracers with higher order VSPs we are able to measure $\rho(r)$ over the range $0.5 < r/R_{1/2} < 2$ and broadly constrain $\beta(r)$. Including proper motion data for all stars gives an even better performance, with $\rho$ and $\beta$ well-measured over the range $0.25 < r/R_{1/2} < 4$. Finally, we test GravSphere on a triaxial mock galaxy that has axis ratios typical of a merger remnant, $[1:0.8:0.6]$. In this case, GravSphere can become slightly biased. However, we find that when this occurs the data are poorly fit, allowing us to detect when such departures from spherical symmetry become problematic.

Tables

Table 1. Table 1: Parameters of the spherical and triaxial mock data that we use in this paper. The columns give, from left to right: the mock data label; the tracer density parameters (see equation 31 ); the dark matter halo density parameters (see equation 31 ); the velocity anisotropy profile; the projected half stellar mass radius; and the long a 𝑎 a , intermediate b 𝑏 b and short c 𝑐 c axis of the figure, where [ 1 , 1 , 1 ] 1 1 1 [1,1,1] corresponds to spherical symmetry. The length units are given in kpc, while ρ 0 subscript 𝜌 0 \rho_{0} has units of 10 7 superscript 10 7 10^{7}  M ⊙  kpc -3 . The mock labelling convention is: < < light profile > > < < Cusp/core > > < < Anisotropy > > , where the < < light profile > > can be Plummer-like (Plum; γ ∗ = 0.1 subscript 𝛾 0.1 \gamma_{*}=0.1 ); cusped (NonPlum; γ ∗ = 1 subscript 𝛾 1 \gamma_{*}=1 ); split-component (Split); or triaxial (Triax); the dark matter halo can be cusped (Cusp; γ D = 1 subscript 𝛾 D 1 \gamma_{\text{D}}=1 ) or cored (Core; γ D = 0 subscript 𝛾 D 0 \gamma_{\text{D}}=0 ); and the < < Anisotropy > > can be isotropic (Iso; β = 0 𝛽 0 \beta=0 ); of Osipkov-Merritt form (Om; see equation 10 ) or tangential (Tan; β = − 0.5 𝛽 0.5 \beta=-0.5 ). For the split-component models, there are two separate tracers moving in the same dark matter halo; these are labelled by the subscript 1 and 2, respectively. For the TriaxCusp mock, the stars and dark matter are assumed to be aligned, with identical axis ratios [ a , b , c ] 𝑎 𝑏 𝑐 [a,b,c] . We project this mock along either the long axis ( X 𝑋 X ) or short axis ( Z 𝑍 Z ), leading to two different R 1 / 2 subscript 𝑅 1 2 R_{1/2} , as marked in the Table.

Label	$[r_{\ast },{\alpha }_{\ast },{\beta }_{\ast },{\gamma }_{\ast }]$	$\left[({\rho }_0,r_{\text{D}},{\alpha }_{\text{D}},{\beta }_{\text{D}},{\gamma }_{\text{D}})\right]$	$\beta \left(r\right)$	$R_{1/2}$	$[a,b,c]$
PlumCuspIso	$[0.25,2,5,0.1]$	$[6.4,1,1,3,1]$	0	0.25	$[1,1,1]$
NonPlumCoreOm	$[0.25,2,5,1]$	$[40,1,1,3,0]$	Eq. 10; $r_a=0.25$	0.20	$[1,1,1]$
NonPlumCuspOm	$[0.1,2,5,1]$	$[6.4,1,1,3,1]$	Eq. 10; $r_a=0.1$	0.079	$[1,1,1]$
PlumCuspTan	$[0.5,0.5,5,0.1]$	$[2.39,2,1,4,1]$	$-0.5$	0.86	$[1,1,1]$
SplitCompCore	${[0.5,2,5,1]}_1$	$[40,1,1,3,0]$	Eq. 10; $r_{a,1}=0.5$	0.39	$[1,1,1]$
	${[1.0,2,5,1]}_2$		Eq. 10; $r_{a,2}=1.0$	0.79
SplitCompCusp	${[0.5,2,5,1]}_1$	$[6.4,1,1,3,1]$	Eq. 10; $r_{a,1}=0.5$	0.39	$[1,1,1]$
	${[1.0,2,5,1]}_2$		Eq. 10; $r_{a,2}=1.0$	0.79
TriaxCusp	$[0.81,0.34,5.92,0.23]$	$[5.5,1.5,1,4,1]$	Eq. 33	0.6 ($X$), 0.8 ($Z$)	$[1,0.8,0.6]$

Equations

$$\frac{\text{d}f}{\text{d}t}=\frac{\partial f}{\partial t}+{\nabla }_{\mathbf{x}}f\cdot \mathbf{v}-{\nabla }_{\mathbf{v}}f\cdot {\nabla }_{\mathbf{x}}\Phi =0,$$

$${\nabla }_{\mathbf{x}}^2\Phi =4\pi G\rho$$

$$\frac{1}{\nu }\frac{\partial }{\partial r}\left(\nu {\sigma }_r^2\right)+\frac{2\beta (r){\sigma }_r^2}{r}=-\frac{GM(<r)}{r^2},$$

$${\sigma }_r^2=⟨v_r^2⟩-⟨v_r{⟩}^2;⟨v_r^n⟩=\int v_r^{n}f{d}^3\mathbf{v};$$

$$\nu =\int f{d}^3\mathbf{v};\beta \equiv 1-\frac{{\sigma }_t^2}{{\sigma }_r^2}$$

$${\sigma }_r^2(r)=\frac{1}{\nu (r)g(r)}{\int }_r^{\infty }\frac{GM(\tilde{r})\nu (\tilde{r})}{{\tilde{r}}^2}g(\tilde{r})\text{d}\tilde{r}$$

$$g(r)=\exp \left(2\int \frac{\beta (r)}{r}\text{d}r\right)$$

$${\sigma }_{\text{LOS}}^2(R)=\frac{2}{{\Sigma }_{\ast }(R)}{\int }_R^{\infty }\left(1-\beta \frac{R^2}{r^2}\right)\frac{\nu (r){\sigma }_r^2(r)r}{\sqrt{r^2-R^2}}\text{d}r,$$

$$\beta (r)={\beta }_0+\left({\beta }_{\infty }-{\beta }_0\right)\frac{1}{1+{\left(\frac{r_0}{r}\right)}^n}$$

$$\beta =r^2/(r^2+r_0^2)$$

$$g(r)=r^{2{\beta }_{\infty }}{\left({\left(\frac{r_0}{r}\right)}^n+1\right)}^{\frac{2}{n}({\beta }_{\infty }-{\beta }_0)}$$

$$\beta (r)={\beta }_0+({\beta }_{\infty }-{\beta }_0)\left[\sum _{i=0}^N{\left(1+{\left(\frac{r_i}{r}\right)}^{n_i}\right)}^{-1}\right]$$

$$g(r)$$

$$\tilde{\beta }=\frac{{\sigma }_r-{\sigma }_t}{{\sigma }_r+{\sigma }_t}=\frac{\beta }{2-\beta }$$

$$\nu =\sum _j^{N_P}\frac{3M_j}{4\pi a_j^3}{\left(1+\frac{r^2}{a_j^2}\right)}^{-5/2}$$

$${\Sigma }_{\ast }(R)=\frac{R^2}{\rho (R)}\sum _j^{N_P}\frac{15M_j}{4\pi a_j^5}{\left(1+\frac{R^2}{a_j^2}\right)}^{-7/2}$$

$$M(r)=M_{\ast }(r)+M_{\mathrm{dm}}(r)$$

$$\rho_{\rm dm}(r)=\left\{\begin{array}[]{ll}\rho_{0}\left(\frac{r}{r_{0}}\right)^{-\gamma_{0}}&r<r_{0}\\ \rho_{0}\displaystyle\prod_{n=0}^{j<N_{\rm dm}}\left(\frac{r_{n+1}}{r_{n}}\right)^{-\gamma_{n}}\left(\frac{r}{r_{j+1}}\right)^{-\gamma_{j+1}}&r_{j}<r<r_{j+1}\\ \end{array}\right.$$

$$r_{j=\mathrm{0...}(N_{\mathrm{dm}}-1)}=[0.25,0.5,1.0,2.0,4.0]R_{1/2}$$

$$v_{s1}$$

$$=$$

$$v_{s2}$$

$$=$$

$${\zeta }_A$$

$$\simeq$$

$${\zeta }_B$$

$$\simeq$$

$${\sigma }_{\text{pmr}}(R)=\frac{2}{{\Sigma }_{\ast }(R)}{\int }_R^{\infty }\left(1-\beta +\beta \frac{R^2}{r^2}\right)\frac{\nu (r){\sigma }_r^2(r)r}{\sqrt{r^2-R^2}}\text{d}r$$

$${\sigma }_{\text{pmt}}(R)=\frac{2}{{\Sigma }_{\ast }(R)}{\int }_R^{\infty }\left(1-\beta \right)\frac{\nu (r){\sigma }_r^2(r)r}{\sqrt{r^2-R^2}}\text{d}r$$

$${\chi }^2$$

$${\rho }_X(r)={\rho }_0{\left(\frac{r}{r_X}\right)}^{-{\gamma }_X}{\left[1+{\left(\frac{r}{r_X}\right)}^{{\alpha }_X}\right]}^{({\gamma }_X-{\beta }_X)/{\alpha }_X}$$

$$r^2\to m^2={\left(\frac{x}{a}\right)}^2+{\left(\frac{y}{b}\right)}^2+{\left(\frac{z}{c}\right)}^2$$

Full text

We find that two populations provide an excellent recovery of $\rho(r)$ in-between their respective $R_{1/2}$. However, even with a total of $\sim 7,000$ tracers, we are not able to well-constrain $\beta(r)$ for either population. By contrast, using 1000 tracers with higher order VSPs we are able to measure $\rho(r)$ over the range $0.5<r/R_{1/2}<2$ and broadly constrain $\beta(r)$. Including proper motion data for all stars gives an even better performance, with $\rho$ and $\beta$ well-measured over the range $0.25<r/R_{1/2}<4$.

Finally, we test GravSphere on a triaxial mock galaxy that has axis ratios typical of a merger remnant, $[1:0.8:0.6]$. In this case, GravSphere can become slightly biased. However, we find that when this occurs the data are poorly fit, allowing us to detect when such departures from spherical symmetry become problematic.

1 Introduction

The mass modelling of stellar systems has a long and rich history dating back to Eddington (1915)’s work on modelling star clusters and Kapteyn (1922)’s seminal model of the Milky Way. It has been used to probe the distribution and nature of dark matter in galaxies and galaxy clusters (e.g. Zwicky, 1937; Kleyna et al., 2001; Łokas & Mamon, 2003; Walker et al., 2009a; Cappellari, 2016); the initial mass function of stars (e.g. Cappellari et al., 2012; Sollima et al., 2012; Napolitano et al., 2014; Shanahan & Gieles, 2015); the structure of star clusters (e.g. Pryor & Meylan, 1993; Watkins et al., 2013; Peuten et al., 2016); the mass of central ‘supermassive’ black holes (e.g. Magorrian et al., 1998; van den Bosch & de Zeeuw, 2010); and even the nature of weak-field gravity (e.g. Famaey & McGaugh, 2012; Lüghausen et al., 2014; McGaugh, 2016).

To a very good approximation, stars in galaxies obey the collisionless Boltzmann equation (e.g. Binney & Tremaine, 2008):

[TABLE]

where $f(\mathbf{x},\mathbf{v})$ is the distribution function of stars in phase space; $\mathbf{x}$ and $\mathbf{v}$ are the position and velocity of the stars, respectively; and $\Phi$ is the gravitational potential.

For galaxies in a ‘steady state’, the $\frac{\partial f}{\partial t}$ term is small and can be discarded. Thus, in principle, given the positions and velocities of many stars, we can directly solve equation 1 for the force field $\nabla_{\mathbf{x}}\Phi$. From this we can derive the mass distribution from the Poisson equation:

[TABLE]

where $G=6.67398\cdot 10^{-11}\text{m}^{3}/\text{kg}\,\text{s}^{2}$ is Newton’s gravitational constant, and the matter density $\rho$ includes all contributions from stars, gas and dark matter. However, in practice solving equation 1 is hard because it involves derivatives of $f$ that is six dimensional and poorly sampled. Several solutions have been proposed for this, including assuming a well-motivated functional form for $f$ (e.g. Eddington, 1915; King, 1966; Hernquist, 1990; Cuddeford, 1991; Heggie & Ramamani, 1995; Gieles & Zocchi, 2015); modelling $f$ as a sum of many individual stellar orbits (e.g. Schwarzschild, 1979; van den Bosch et al., 2008; Breddels & Helmi, 2013); or dynamically evolving an $N$-body simulation of $f$ towards a target distribution (e.g. Syer & Tremaine, 1996; de Lorenzi et al., 2007; Dehnen, 2009; Hunt & Kawata, 2013). In this paper, we focus on a fourth approach: moment methods. These were first proposed by Jeans (1922) and involve integrating equation 1 over each component of the velocity to create a set of three moment equations. There are several well-known disadvantages to this approach (e.g. Binney & Tremaine, 2008): the hierarchy of moment equations is not closed; there is no guarantee that their solution will correspond to a physical distribution function; and the method typically requires the data to be binned, leading to the loss of potentially useful information. However, the key advantages are that the form of $f$ need not be specified and that the Jeans equations can be solved very rapidly. This latter is particularly important for the method we present in this work, where we employ a large number of free parameters that we marginalise over when fitting data.

In this paper, we concern ourselves with non-rotating near-spherical stellar systems that are a reasonable approximation to star clusters, spheroidal and elliptical galaxies and galaxy clusters. Integrating equation 1 over the radial velocity $v_{r}$ and assuming no rotation, we arrive at the ‘spherical Jeans equation’ (Binney, 1980; Binney & Tremaine, 2008):

[TABLE]

where:

[TABLE]

$M(<r)$ is the total cumulative mass as a function of radius $r$; $\nu$ and $\beta$ are the ‘tracer density’ and ‘velocity anisotropy’, respectively:

[TABLE]

and $\sigma_{r}$ and $\sigma_{t}$ are the radial and tangential velocity dispersion, respectively.

The tracer density $\nu(r)$ describes the radial density profile of a population of massless tracers moving in the gravitational potential of a spherical mass distribution $M(r)$. These could be massless ‘stars’ moving in a star cluster or galaxy, or massless galaxies moving in a galaxy cluster. If the tracers correspond to stars/galaxies that do have significant gravitating mass, then this is subsumed into the definition of $M(r)$. The velocity anisotropy $\beta(r)$ describes the orbital structure of the stellar system, where $\beta=0$ corresponds to an isotropic distribution; $\beta=1$ is a fully radial distribution; and $\beta=-\infty$ is a fully tangential distribution.

Integrating both sides of equation 3 gives the radial velocity dispersion as function of radius $r$ (e.g. van der Marel, 1994; Mamon & Łokas, 2005):

[TABLE]

where:

[TABLE]

However, typically, only line-of-sight velocities are observable. Projecting equation 6 along the line of sight, we obtain (Binney & Mamon, 1982):

[TABLE]

where $\Sigma_{*}(R)$ denotes the tracer surface mass density at projected radius $R$; and $\sigma_{\text{LOS}}(R)$ is the line-of-sight velocity dispersion of these tracer stars.

From equation 8, it is clear that the velocity anisotropy $\beta(r)$ trivially degenerates with $\sigma_{r}(r)$ and, therefore, with the cumulative mass distribution $M(<r)$ that we would like to measure. This is the well-known $M-\beta$ degeneracy (e.g. Merrifield & Kent, 1990; Wilkinson et al., 2002; Łokas & Mamon, 2003; de Lorenzi et al., 2009). Since the radial density profile $\rho(r)$ uniquely maps to the enclosed mass $M(<r)$, this can also be cast as an equivalent density-velocity anisotropy degeneracy that we refer to from here on as the ‘$\rho-\beta$ degeneracy’.

Several methods have been proposed in the literature to break the $\rho-\beta$ degeneracy. Saglia et al. (2000) suggest using multiple tracer populations, typically split by metallicity, colour, or stellar type (and see also Battaglia et al. 2008; Walker & Peñarrubia 2011; Amorisco & Evans 2011; Napolitano et al. 2014). Wojtak et al. (2009) suggest performing a joint analysis on similar galaxy clusters to recover a statistical estimate of $\beta(r)$ for the whole population. van der Marel & Franx (1993) suggest measuring the shape of the line-of-sight velocity distribution function, using deviations from Gaussianity to constrain $\beta(r)$ (and see also Wilkinson et al. 2002; Breddels et al. 2013; Mamon et al. 2013). Łokas & Mamon (2003) propose using higher order moments of $f$ – e.g. the kurtosis – to break the degeneracy. However, this can only work if assumptions are made about the form of $\beta(r)$ and its higher order counterparts (Merrifield & Kent, 1990; Richardson & Fairbairn, 2013). For this reason, Merrifield & Kent (1990) propose instead using fourth order ‘Virial Shape Parameters’ (VSPs; and see Richardson & Fairbairn 2014). These are two fourth order moments integrated over radius $r$ to produce two scalar numbers (see §2). The advantage of these VSPs is that they depend only on $\beta(r)$. Finally, proper motion data can be used to break the degeneracy (e.g. Strigari et al., 2007; van der Marel & Anderson, 2010; Watkins et al., 2013).

While there is no shortage of methods for breaking the $\rho-\beta$ degeneracy, there is often a wide dispersion of conflicting results in the literature. For example, for the Sculptor dwarf spheroidal galaxy, split population methods favour a constant density dark matter core (Battaglia et al., 2008; Walker & Peñarrubia, 2011; Amorisco & Evans, 2011; Agnello & Evans, 2012); a VSP method favours a cusp (Richardson & Fairbairn, 2014); and Breddels et al. (2013) have argued, using a Schwarzschild distribution function method, that the data are not sufficiently constraining to be able to tell (and see Breddels & Helmi 2014). Similarly, there is an ongoing debate about the presence or absence of an ‘intermediate’ mass black hole in the globular clusters Omega-Centauri (e.g. Noyola et al., 2008; van der Marel & Anderson, 2010; Noyola et al., 2010; Watkins et al., 2013; Zocchi et al., 2016) and G1 (e.g. Baumgardt et al., 2003; Gebhardt et al., 2005; Miller-Jones et al., 2012), while conflicting results are also found in the elliptical galaxy modelling community (e.g. Romanowsky et al., 2003; de Lorenzi et al., 2009).

To shed light on the above discrepancies, in this paper we develop a new non-parametric spherical Jeans code, GravSphere, that is designed to recover $\beta(r)$ and $\rho(r)$ under the assumptions only of spherical symmetry and a steady-state distribution function. We test GravSphere on a large suite of spherical and triaxial mock data, probing different tracer distributions, $\beta(r)$ profiles, and $\rho(r)$ profiles. Our specific goal is to assess – for ideal data – how well different methods for breaking the $\rho-\beta$ degeneracy work. We focus in particular on three methods: split populations; VSPs; and proper motions. In a companion paper, we will describe these mock data in more detail and compare and contrast GravSphere with over methods in the literature (Read et al. in prep.).

This paper is organised as follows. In §2, we describe the GravSphere method. In §3, we briefly describe the mock data that we use in this paper that is a subset of the ‘default’ spherical and triaxial suite from the Gaia Challenge wiki site (Walker & Peñarrubia, 2011; Dehnen, 2009)111http://astrowiki.ph.surrey.ac.uk/dokuwiki/.. These mock data are publicly available and described in detail in Read et al. in prep. In §4, we test GravSphere on these mock data, assuming only line-of-sight data (§4.1); including split populations (§4.2); testing the utility of VSPs (§4.3); and adding in proper motion data (§4.4). In §4.5, we assess how well GravSphere performs on triaxial data for which it should become biased. Finally, in §5 we present our conclusions.

2 Method

In this section, we describe our new GravSphere ‘non-parametric’ Jeans modelling code. GravSphere primarily solves equation 8 to calculate the projected velocity dispersion profile, $\sigma_{\text{LOS}}(R)$, for a set of tracers moving within a spherical mass distribution $M(r)$. These ‘tracers’ can be stars in star clusters or galaxies; globular clusters or planetary nebulae, or even galaxies moving in a galaxy cluster; they do not need to self consistently generate the full gravitational potential. At minimum, we assume that we have measurements of $\sigma_{\text{LOS}}(R)$ and the surface brightness profile of the tracers, $\Sigma_{*}(R)$. We also explore having $N_{t}$ distinct tracers that each have their own $\sigma_{\text{LOS},i}(R),\Sigma_{*,i}(R)$ and velocity anisotropy $\beta_{i}(r)$, while moving in the same mass distribution $M(r)$. And we consider the utility of higher order ‘Virial Shape Parameters’ (VSPs) and proper motion data. Otherwise, our goal is to assume only spherical symmetry and a steady-state distribution function. We describe our method as ‘non-parametric’ because we have sufficient freedom in our functional forms for $\nu(r)$ and $M(r)$ that we can have many more parameters than data constraints. These are then marginalised over to explore model degeneracies (see e.g. the discussion in Coles et al. (2014) on this point).

In the following subsections, we describe our functional forms for $\beta(r)$ (§2.1), $\nu(r)$ (§2.2) and $M(r)$ (§2.3); how we wrap in additional data from multiple populations (§2.4), VSPs (§2.5) and proper motions (§2.6); and how we compare our models to data (§2.7).

2.1 The velocity anisotropy

We assume an expression for the velocity anisotropy $\beta(r)$ that ensures broad generality while making the function $g(r)$ (equation 7) analytic:

[TABLE]

where $\beta_{0}$ is the inner asymptotic anisotropy (in the limit $r\rightarrow 0$); $\beta_{\infty}$ is the outer asymptotic anisotropy (in the limit $r\rightarrow\infty$); $r_{0}$ is a transition radius; and $n$ controls the sharpness of the transition. Notice that equation 9 is identical to the Osipkov-Merritt anisotropy profile (Osipkov, 1979; Merritt, 1985):

[TABLE]

in the limit $n=2$, $\beta_{0}=0$ and $\beta_{\infty}=1$.

Putting equation 9 into equation 7, we derive:

[TABLE]

More general forms for $\beta(r)$ can be readily implemented if the data quality warrant it222In an early version of GravSphere, we employed a fully non-parametric model for $\beta(r)$ that was linearly interpolated over a number of discrete bins. This was problematic for two reasons. Firstly, it required us to numerically solve the integral in equation 7. This must be solved with care since it is then used as the exponent of an exponential. Secondly, without regularisation priors, this led to a wildly fluctuating $\beta(r)$ that varied discretely from bin to bin. For these reasons, we moved instead to a smooth analytic form for $\beta(r)$.. Indeed, the above form can be generalised to have $N$ radii at which $\beta(r)$ transitions to a different value:

[TABLE]

which also gives an analytic solution to equation 7:

[TABLE]

In this paper, we explore only the simpler form for $\beta$ in equation 9.

Equation 9 is a useful form for the velocity anisotropy because it appears directly in equation 8 and it avoids us needing to numerically integrate equation 7. However, it is problematic for model fitting because $\beta_{0}$ and $\beta_{\infty}$ are defined over an infinite range ($-\infty<\beta_{0}<1$). For this reason, throughout this paper we work instead with a symmetrised anisotropy parameter (Read et al., 2006):

[TABLE]

where $\tilde{\beta}=-1$ corresponds to full tangential anisotropy; $\tilde{\beta}=1$ to full radial anisotropy; and $\tilde{\beta}=0$ to isotropy. By using $\tilde{\beta}$ instead of $\beta$, we can set flat priors on $-1<\tilde{\beta}_{0,\infty}<1$ that give equal weight to tangentially and radially anisotropic models. (The transformation from $\tilde{\beta}_{0,\infty}$ to $\beta_{0,\infty}$ is given by equation 14.)

2.2 The tracer density profile

For the tracer density, we use a sum of $N_{P}$ Plummer spheres, similarly to Rojas-Niño et al. (2016):

[TABLE]

where $M_{j}$ and $a_{j}$ are the mass and scale length of each individual component. (These components should not be thought of as physical, but rather as terms in a density expansion.)

This has the advantage that the tracer surface density (that we will compare with observational data (§2.7) and that appears in equation 8) is then analytic:

[TABLE]

We fit the above expansion in advance of running our full fit using the Python package LMFIT (Newville et al., 2014)333http://zenodo.org/record/11813#.V2P8khV97eQ.. We then allow $M_{j}$ and $a_{j}$ to vary in the final fits by up to 50% of their best-fit LMFIT values. This makes the method substantially more efficient since $\Sigma_{*}(R)$ is usually much better constrained than the stellar kinematics. It also avoids sampling trivial degeneracies between the individual Plummer components in equation 15. By default, we use ${N_{P}}=3$ Plummer components, corresponding to $6$ free parameters, but arbitrarily many can be employed if the data require it. (Note that the above expansion is similar to the Multi-Gaussian Expansion proposed by Emsellem et al. (1994). However, since both dwarf galaxies and globular star clusters are well-approximated by Plummer spheres, for many applications the above expansion is particularly efficient; e.g. Plummer 1911; Irwin & Hatzidimitriou 1995; Martin et al. 2008.)

2.3 The mass profile

We split the enclosed mass profile into a contribution from stars, $M_{*}$, and a contribution from dark matter $M_{\rm dm}$:

[TABLE]

The stellar component is assumed to have the same form as the photometric light profile (fit similarly to the tracer density, as described in §2.2) with a mass to light ratio parameter $\Upsilon$ that we fit for. For the mock data in this paper, $\Upsilon=0$ and the stars do not contribute to the gravitational potential (see §3).

The dark matter component is described by a sequence of power laws defined in $N_{\rm dm}$ radial bins with bin edges $r_{j}$:

[TABLE]

where $\rho_{0}$ is a normalisation parameter that controls the mass of the dark matter halo, and $\gamma_{j}$ sets the logarithmic density slope of each bin $j$. The cumulative mass profile $M_{\rm dm}(r)$ then follows analytically from integrating equation LABEL:eqn:rhodm.

The above form allows a wide range of density profiles to be explored while ensuring that these are positive definite444This is important, for example, in satellite dwarf galaxies that may have their inner dark matter ‘heated’ by bursty star formation (e.g. Navarro et al., 1996a; Read & Gilmore, 2005; Pontzen & Governato, 2012) and their outer dark matter profile sculpted by tides (e.g. Read et al., 2006). These two processes introduce two new scales to the dark matter halo on which we expect it to deviate from the ‘NFW’ profile found in $\Lambda$ Cold Dark Matter ‘dark matter only’ structure formation simulations (Navarro et al., 1996b).. It also defines the dark matter density in terms of its local logarithmic slope, $\gamma_{j}$. This is a quantity of key interest for probing the nature of dark matter, but often difficult to extract because it requires a double differential of $M(<r)$ that can become very noisy. With the above functional form, we fit directly for $\gamma_{j}$, thereby avoiding this problem. Finally, we are able to place useful theory priors on $\gamma_{j}$ that it must be positive and monotonically decrease outwards. The former ensures that the density monotonically decreases outwards; the latter prevents the dark matter profile from having any inflection points. If such inflection points are expected in a given dark matter theory, then this prior can be relaxed (e.g. Young et al., 2016; Gonzáles-Morales et al., 2016).

For a single tracer with only line-of-sight velocities, $M(<r)$ is maximally constrained at $\sim R_{1/2}$ (e.g. Walker et al., 2009b; Wolf et al., 2010). For this reason, we use as default $N_{\rm dm}=5$ dark matter mass bins spaced logarithmically around $R_{1/2}$:

[TABLE]

More bins can be added if the data quality warrant it; we find that, for the mocks presented in this paper, we are not sensitive to this choice.

By default, for a single tracer component, our model has 6 parameters to describe the light profile (that are typically very well constrained by the data); 4 parameters to describe $\beta(r)$; 6 parameters to describe the mass profile ($\rho_{0}$ and 5 $\gamma_{j}$ for each radial bin); and one parameter to describe the mass to light ratio of the stars $\Upsilon$. This gives 17 parameters in total.

2.4 Multiple tracer populations

Our method can cope naturally with $N_{t}$ sets of tracer stars moving in the same potential. As discussed in §1, if the $N_{t}$ populations have distinct $R_{1/2,i...N_{t}}$, it has been argued in the literature that they can break the $\rho-\beta$ degeneracy and allow for a measurement of $\rho(r)$.

Each tracer $i$ has a distinct light profile $\nu_{i}(r)$ and velocity anisotropy profile $\beta_{i}(r)$. These lead to distinct projected light ($\Sigma_{*,i}(R)$) and velocity dispersion ($\sigma_{\text{LOS},i}(R)$) profiles that can be compared with data. Thus, for each tracer we must solve an additional equation 8.

Since there are now multiple $R_{1/2,i}$, we also update our binning strategy for $M_{\rm dm}$, placing one bin at each $R_{1/2,i}$, two at $[0.5,1]R_{1/2,{\rm min}}$ and two at $[2.0,4.0]R_{1/2,{\rm max}}$. Thus, our default method for $N_{t}$ tracers requires a total of $10N_{t}$ parameters to describe $\nu_{i}(r)$ and $\beta_{i}(r)$; and $6+N_{t}$ parameters to describe $M(r)$ (one for $\rho_{0}$; $4+N_{t}$ for the $\gamma_{j}$; and one for $\Upsilon$). This gives $11N_{t}+6$ parameters in total, which can become substantially more expensive than the single tracer case.

2.5 Including ‘Virial Shape Parameters’ (VSPs)

Following Merrifield & Kent (1990) and Richardson & Fairbairn (2014), we also include the ability to use higher order moments of the velocity distribution via the fourth order ‘Virial Shape Parameters’ (VSPs):

[TABLE]

and

[TABLE]

The key advantage of these VSPs is that they involve fourth-order moments of the line of sight velocities $\langle v_{\rm LOS}^{4}\rangle$ (equations 21 and 23), but by integrating over projected radius $R$, $v_{s1}$ and $v_{s2}$ depend only on $\beta$ and not on its fourth-order counterparts (see equations 21 and 23 and Merrifield & Kent, 1990; Richardson & Fairbairn, 2014). Thus, $v_{s1}$ and $v_{s2}$ allow us to obtain constraints on $\beta$ via line of sight velocities alone, breaking the $\rho-\beta$ degeneracy.

However, we take a slightly different approach to Richardson & Fairbairn (2014) in comparing $v_{s1}$ and $v_{s2}$ to data. They suggest recasting $v_{s1}$ and $v_{s2}$ as dimensionless estimators:

[TABLE]

[TABLE]

where $N_{\rm tot}=\int_{0}^{\infty}\Sigma_{*}RdR$ is $1/2\pi$ times the total number of stars, and $N_{s}$ is the actual number of observed stars.

The problem with $\zeta_{A}$ and $\zeta_{B}$ is twofold. Firstly, when calculating 25 and 27 from the data, we must take a noisy fourth order estimator and divide it by a noisy second order estimator, amplifying noise in the data. Secondly, the approximations in equations 25 and 27 allow for the inclusion only of errors in the velocity measurements. This could lead to artificially small errors on $\zeta_{A}$ and $\zeta_{B}$. For these reasons, in GravSphere we work directly with $v_{s1}$ and $v_{s2}$, solving equations 21 and 23 by numerically integrating the observed $\langle v_{\rm LOS}^{4}\rangle$ over our best-fit model $\Sigma_{*}$. In this way, we can estimate the uncertainties on $v_{s1}$ and $v_{s2}$ by Monte-Carlo sampling the errors on both $\langle v_{\rm LOS}^{4}\rangle$ and, if significant, $\Sigma_{*}$. We describe this procedure in more detail in 2.7. We compare our estimators for $v_{s1}$ and $v_{s2}$ with those used in Richardson & Fairbairn (2014) in Appendix B.

2.6 Including proper motions

Where proper motion data are available, we add two further constraint equations for the radially and tangentially averaged projected proper motions, $\sigma_{\rm pmr}$ and $\sigma_{\rm pmt}$ (e.g. Strigari et al., 2007; van der Marel & Anderson, 2010):

[TABLE]

[TABLE]

Equations 8, 28 and 29 each have a different dependence on $\beta$. For an isotropic distribution ($\beta=0$) all three become equal, while any anisotropy will cause them to differ. Thus, by fitting all three simultaneously to data for $\sigma_{\text{LOS}}$, $\sigma_{\text{pmr}}$ and $\sigma_{\text{pmt}}$, we can obtain tight constraints on both $\beta$ and $\rho$, thereby breaking the $\rho-\beta$ degeneracy (see §1).

2.7 Fitting the model to data

To compare our model to data, we bin the data in projected radius $R$. Following Walker et al. (2009b), we use $N_{\rm bin}=\sqrt{N_{*}}$ bins where $N_{*}$ is the total number of stellar tracers. This gives us $\sqrt{N_{*}}$ stars in each bin, ensuring good sampling and equal Poisson error for each bin, while not compromising too much on the radial coverage. (In Appendix A, we show that we are not sensitive to this choice.) The error on each bin includes the effect of Poisson sampling and measurement error, calculated as in Pryor & Meylan (1993). From here on, we will denote binned data quantities with an line over the top. For example, the binned projected surface density of a tracer population $i$ will be denoted $\overline{\Sigma}_{*,i}$ to differentiate it from the GravSphere model surface density $\Sigma_{*,i}$ in equation 16.

Our binned data comprise the projected surface density of each set of tracers $\overline{\Sigma}_{*,i}$; the projected velocity dispersion of these tracers $\overline{\sigma}_{\text{LOS},i}$ and, if proper motions are available, the radial and tangential projected dispersions $\overline{\sigma}_{\text{pmr},i}$ and $\overline{\sigma}_{\text{pmt},i}$. (In general, each of $\overline{\Sigma}_{*,i}$, $\overline{\sigma}_{\text{LOS},i}$, $\overline{\sigma}_{\text{pmr},i}$ and $\overline{\sigma}_{\text{pmt},i}$ will have its own bins.)

To estimate the VSPs for a given tracer population, we use a novel technique. First, we calculate the binned fourth order line-of-sight velocity moment $\overline{\langle v_{\rm LOS}^{4}\rangle}$, using the same bins as for $\overline{\sigma}_{\text{LOS}}$. We estimate the errors on $\overline{\langle v_{\rm LOS}^{4}\rangle}$ by Monte-Carlo sampling the velocity error on each star 1000 times and then adding the measurement variance on $\overline{\langle v_{\rm LOS}^{4}\rangle}$ in quadrature with the Poisson error. We then numerically integrate equations 21 and 23 to calculate $\overline{v}_{s1}$ and $\overline{v}_{s2}$. For these integrals, we use our best-fit model $\Sigma_{*}(R)$, assuming that the error on this is negligible as compared with the error on $\overline{\langle v_{\rm LOS}^{4}\rangle}$. To obtain an error on $\overline{v}_{s1}$ and $\overline{v}_{s2}$, we then perform a Monte-Carlo sampling of the errors on $\overline{\langle v_{\rm LOS}^{4}\rangle}$, repeating our numerical integration of equations 21 and 23 1000 times. The variance on these integrals is then taken as a Gaussian error on $\overline{v}_{s1}$ and $\overline{v}_{s2}$. (We verified that the probability distribution functions of $\overline{v}_{s1}$ and $\overline{v}_{s2}$ that result from these integrals are indeed close to Gaussian and that our error bars are not sensitive to the choice of 1000 Monte-Carlo samples.)

We then fit our model to these binned data to obtain a joint $\chi^{2}$:

[TABLE]

where $\overline{\Sigma}_{*,ij}$ is the observed surface density of a given tracer $i$ in bin $j$; $\Sigma_{*,ij}$ is the model surface density of tracer $i$ evaluated at bin $j$; $\sigma_{\Sigma_{*},ij}$ is the error on that bin; and similar $\chi^{2}$ terms are added for $\overline{\sigma}_{\text{LOS},i}$, $\overline{\sigma}_{\text{pmr},i}$ and $\overline{\sigma}_{\text{pmt},i}$ (if available), and $\overline{v}_{s1,i}$ and $\overline{v}_{s2,i}$ (if available).

We fit our model to the data using the emcee affine invariant Markov Chain Monte Carlo (MCMC) sampler from Foreman-Mackey et al. (2013). We assume uncorrelated Gaussian errors such that the Likelihood function is given by $\mathcal{L}=\exp(-\chi^{2}/2)$. We use as default 1000 walkers, each generating 5000 models and we throw out the first half of these as a conservative ‘burn in’ criteria. (These numbers are substantially larger than usually employed due to the high dimension of our parameter space.) We explicitly checked that our results are converged by running more models and examining walker convergence.

In setting our priors, our basic philosophy is to use flat priors on parameters that are well constrained and/or dimensionless, and logarithmic priors on everything else. Specifically, we use flat priors on the dark matter logarithmic slope parameters $0<\gamma_{j}<3$; a flat logarithmic prior on the dark matter density normalisation $6<\log_{10}[\rho_{0}/({\rm M}_{\odot}\,{\rm kpc}^{-3})]<10$; a flat prior on the stellar mass to light ratio $\Upsilon$ over a range $\pm 25$% of the stellar population synthesis model value (for the mock data we use in this paper $\Upsilon=0$); flat priors on the scale length parameters describing the tracer surface density $a_{j}$ within 50% of their pre-determined best-fit values; a flat logarithmic prior on the mass parameters describing the tracer surface density $\log_{10}[M_{j}]$, also within 50% of their pre-determined best-fit values; a flat logarithmic prior on the velocity anisotropy scale length $-0.3<\log_{10}[r_{0}/R_{1/2}]<0.3$; and a flat prior on the transition rate parameter $1<n<3$. We use priors on the inner and outer asymptotic velocity anisotropy parameters $\beta_{0}$ and $\beta_{\infty}$, as described in §2.1. Our results are not sensitive to these prior choices. Finally, we do not include any density slope-anisotropy prior. Such a prior could be used to ensure that the distribution function is everywhere positive (e.g. An & Evans, 2006; Ciotti & Morganti, 2010). We will explore this in future work.

3 The mock data

We take mock data from the ‘default spherical-triaxial suite’ of the Gaia Challenge wiki site555http://astrowiki.ph.surrey.ac.uk/dokuwiki/.. The full spherical-triaxial suite is described in detail in Read et al. in prep. We have run all of the mock data in the default suite using GravSphere. However, for brevity, we show results here for a representative subset of mocks as detailed in Table 1. These are chosen to highlight both the best and worst performance of GravSphere over all of the mock data.

The spherical mocks assume that both the tracer distribution and their host dark matter halo are given by a split power law form (Zhao, 1996):

[TABLE]

where $\rho_{0}$ is a normalisation parameter; $r_{X}$ sets the scale length; $\gamma_{X}$ is the inner asymptotic logarithmic slope; $\beta_{X}$ is the outer asymptotic logarithmic slope; and $\alpha_{X}$ controls how sharp the transition is at $r_{X}$ between $\gamma_{X}$ and $\beta_{X}$. For the tracers, we use the notation: $\nu=\rho_{*}(r,\nu_{0},r_{*},\alpha_{*},\beta_{*},\gamma_{*})$; for the dark matter we write similarly: $\rho_{\text{D}}=\rho_{\text{D}}(r,\rho_{0},r_{\text{D}},\alpha_{\text{D}},\beta_{\text{D}},\gamma_{\text{D}})$.

For the radially anisotropic mocks, the tracers are set up in equilibrium inside their host dark matter halo assuming an Osipkov-Merritt distribution function (Osipkov, 1979; Merritt, 1985), as in Walker & Peñarrubia (2011). This has an anisotropy profile $\beta(r)$ given by equation 10, parameterised by the anisotropy radius $r_{a}$. The tangentially anisotropic mocks are set up using the ‘made to measure’ code from Dehnen (2009). These also assume the form given in equation 31 for the tracer and dark matter density profiles. The velocity anisotropy is, however, assumed to be a constant $\beta=-0.5$.

In addition to the above ‘single population’ mocks, we consider two ‘split population’ mocks that have two sets of tracers stars with distinct $R_{1/2}$ moving in the same dark matter halo. Like the single population mocks, these models assume an Osipkov-Merritt distribution function (Walker & Peñarrubia, 2011), where both tracers have a density profile of the form given in equation 31, but with different parameters and different numbers of tracer stars. We report the model parameters for these mocks in Table 1. For this paper, we assume that we can perfectly split these two populations into their constituent stars.

Finally, we consider a triaxial mock galaxy for which GravSphere (that assumes spherical symmetry) should become biased. This was set up, similarly to the tangential mocks, using the Dehnen (2009) made to measure code. As with the other mocks, the tracer stars and dark matter halo are assumed to be of the form given in equation 31. However, now the radius $r$ is replaced by an ellipsoidal radius parameter $m$:

[TABLE]

where $a>b>c$ set the long, intermediate and short axis of the triaxial figure. The mock we use here has $b/a=0.8$ and $c/a=0.6$, consistent with the typical triaxiality that occurs following the merger of two spherical galaxies (e.g. Gualandris et al., 2017). The tracer stars and their host dark matter halo were assumed to have exactly the same triaxiality. The velocity anisotropy profile for the triaxial mocks is given by:

[TABLE]

where $r_{t}=0.81$ kpc, $\beta_{0}=0$ and $\beta_{\infty}=0.5$, such that these models are isotropic at $r=0$ and become radially anisotropic at radii $r>r_{t}$.

The parameters and labels that we use for all of our mock data are reported in Table 1. In all mocks, the stars are massless tracers that do not contribute to the gravitational potential666For this reason, we do not report the value of $\nu_{0}$ in Table 1 since this parameter is arbitrary.. For the single component mocks, we consider 1000 and 10,000 tracer stars; for the split component mocks, we use $[1877,4703]$ stars for the inner/outer tracer population in SplitCompCore, and $[3075,4030]$ for SplitCompCusp. In this way, the total number of stars for the split component mocks is comparable to the 10,000 star case for our single component tests. In all cases, we assume excellent data quality with constant errors of $2$ km s*-1* on the measured velocities (line of sight or proper motion) for each star.

4 Results

4.1 Line of sight velocities only

We first consider how well GravSphere performs when fitting only line-of-sight velocities, $\overline{\sigma}_{\text{LOS}}(R)$, and the surface density of tracer stars $\overline{\Sigma}_{*}(R)$. The results are shown in Figures 1 (top panels), 2 (using 1000 tracer stars) and 3 (using 10,000 tracer stars).

In Figure 1 (top panels), we show how well $\overline{\sigma}_{\text{LOS}}(R)$ and $\overline{\Sigma}_{*}(R)$ are fit by GravSphere for the NonPlumCoreOm mock with 1000 tracers (GravSphere performs similarly well for the other mocks and so we omit these plots for brevity). In Figure 2, we show the radial density profile $\rho(r)$ (left), logarithmic density slope $d\ln\rho/d\ln r$ (middle) and symmetrised velocity anisotropy $\tilde{\beta}$ (right) for the PlumCuspIso (top), NonPlumCoreOm (middle) and PlumCuspTan (bottom) mocks (see §3 and Table 1 for details of these mock data sets) that result from this fit. In Figure 3, we show the same but for 10,000 tracers.

As can be seen, in all three mocks, we obtain a good estimate of the density at $R_{1/2}$, but a poor measure of both the radial profile of $\rho(r)$ and the velocity anisotropy profile $\tilde{\beta}(r)$. This is true even when using a rather optimistic 10,000 tracer stars (Figure 3), highlighting the familiar $\rho-\beta$ degeneracy (see §1 and Walker et al. 2009b; Wolf et al. 2010). In particular, the logarithmic derivative of the density, $d\ln\rho/d\ln r$ (middle panels), is poorly recovered, with only very weak constraints inside $r<R_{1/2}$ at 95% confidence (light grey shaded regions). (Notice that this parameter moves in discrete steps. This owes to our parameterised form for $\rho_{\text{D}}$ that we described in detail in §2.3.)

Finally, notice that our recovery for the NonPlumCoreOm mock (compare the shaded regions with the blue dashed lines in Figure 2, middle row) is biased, with the bias becoming even slightly worse as we move from 1000 to 10,000 tracers (Figure 3, middle row). This is our worst mock from the whole default spherical-triaxial suite (§3) which is why we have included it here. It is also the only mock from the full suite for which we become biased in this way. The problem occurs because this particular mock is both very cored (with $0<\gamma(r)<0.5$ for $r<R_{1/2}$) and very radially biased. Furthermore, the radial bias sets in only far beyond $R_{1/2}$ where there are few tracer stars (see Figure 2, middle row, middle and rightmost panels). In this case, the hypervolume of solution space remains large due to the $\rho-\beta$ degeneracy, but the correct solution lies right on the edge of this space. Correct models are found within our chains, but they are rare. This is not because these models are a poor representation of the data, but rather because there are many more models that are equally good representations of the data and that have $\tilde{\beta}<0$. This leads to the tangential bias on $\tilde{\beta}$ that can be seen in Figures 2 and 3, with an associated bias on $\rho$ that it becomes too steep.

The above presents us with a good opportunity to remind the reader that the shaded regions in Figure 2 are marginalised histograms of the model parameters in radial bins. They do not represent any single model and, indeed, it can be the case that no model actually passes through the black median line (though in practice, this is not the case for the mock data we consider in this paper). These marginalised histograms are a convolution of both the data and our prior. Where the data constraints are poor, as is the case when only line-of-sight velocities are available, then the shaded regions that result depend as much on our choice of priors as the data. (A prior, for example, that $\beta_{\infty}=1$ would remove the bias from the NonPlumCoreOm mock. But it would introduce a bias into the isotropic and tangential models.)

Ultimately, the only solution to the above bias problem is improved data. In the following subsections, we explore three methods that have been proposed in the literature to break the $\rho-\beta$ degeneracy by introducing additional data constraints: multiple tracer populations with distinct $R_{1/2}$ (§4.2); VSPs (§4.3) and proper motions (§4.4).

4.2 Split populations

As discussed in §1, a popular method for breaking $\rho-\beta$ degeneracy is to simultaneously model multiple tracers with distinct $R_{1/2}$ (e.g. Saglia et al., 2000; Battaglia et al., 2008; Walker & Peñarrubia, 2011; Amorisco & Evans, 2011; Napolitano et al., 2014). Each tracer $i$ provides a robust measurement of $\rho(R_{1/2,i})$ (c.f. §4.1), allowing the density profile to be pinned down at multiple radii.

While the idea of using split populations to break the $\rho-\beta$ degeneracy is clear, there has been some debate in the literature over its efficacy (e.g. Battaglia et al. 2008; Walker & Peñarrubia 2011; Breddels et al. 2013; Breddels & Helmi 2014 and see the discussion in §1). In this context, GravSphere is particularly useful since it assumes only that the system is close to spherical, that the tracers are in Virial equilibrium, and that the tracers have been perfectly separated (see §2).

In Figure 4, we show results for our two split population mocks: SplitCompCusp and SplitCompCore (see Table 1). The top row shows an example of the data fit for the SplitCompCusp case (the SplitCompCore mock is fit similarly well and so we have omitted these plots for brevity). The lines and shaded regions are as in Figure 1, except that now there are two vertical lines to mark the $R_{1/2}$ of the inner (blue) and outer (red) population, while the grey and magenta shaded bands show the 68% and 95% confidence intervals of our fits for the inner and outer population, respectively. The middle and bottom rows show the radial density profile $\rho(r)$; logarithmic slope of the density profile $d\ln\rho/d\ln r$; and symmetrised velocity anisotropy parameter $\tilde{\beta}$ for the SplitCompCusp mock (middle) and SplitCompCore mock (bottom). The lines and shaded regions are as in Figure 2, but with the inner (vertical blue) and outer (vertical red) $R_{1/2}$ of the two populations marked. In the rightmost panels, we show results for the inner $\tilde{\beta}_{1}$ (grey shaded region) and outer $\tilde{\beta}_{2}$ (magenta shaded region); the correct solutions are marked by the blue and red dashed lines, respectively.

As can be seen from Figure 4, using two populations performs as advertised, dramatically improving the recovery777Since the total number of tracers in these split mocks is 6,580 for SplitCompCore and 7,105 for SplitCompCusp (see Table 1 and §3), we should compare the performance of GravSphere on these with the 10k single component models in Figure 3. of $\rho(r)$. We obtain an excellent recovery of $\rho(r)$ in-between $R_{1/2,1}$ and $R_{1/2,2}$, with tight constraints also over the range $0.5R_{1/2,1}<r<2R_{1/2,2}$. Inside $r<0.5R_{1/2,1}$ the constraints become poorer again and we are not able to well-measure the very inner logarithmic cusp slope. (Note, however, that this does depend on our choice of mass model and prior. If we assume a single power law inside $R_{1/2,1}$ then this will be very well constrained, as pointed out by Walker & Peñarrubia 2011. Similarly, our constraints are improved by our monotonicity prior on $d\ln\rho_{\rm dm}/d\ln r$ (§2.3) since this forbids an inner density slope that is steeper than the outer density slope.)

Interestingly, however, while the split component data provide an excellent recovery of $\rho(r)$, $\tilde{\beta}_{1}$ and $\tilde{\beta}_{2}$ are poorly recovered (see Figure 4, middle and bottom rows, rightmost panels). The inner population is reasonably well recovered inside $r<R_{1/2,1}$ within our 95% confidence intervals, but for $r>R_{1/2,1}$ we favour only weakly radially anisotropic models. Similarly, the outer population is poorly recovered, with a weak tangential bias on $\tilde{\beta}_{2}$ for $r<R_{1/2,2}$ and only a hint of radial anisotropy detected at larger radii. It is interesting that despite this $\rho(r)$ is not biased, suggesting that it benefits from the torsioning of one population against the other. We defer a discussion of how GravSphere performs when populations are imperfectly split to future work.

4.3 Virial Shape Parameters

Split populations appear to be a promising avenue for breaking the $\rho-\beta$ degeneracy. However, unambiguously split populations are not always readily available. Furthermore, as we showed in §4.2, they do not yield strong constraints on $\beta$. In this subsection, we consider instead the utility of the fourth order ‘Virial Shape Parameters’ (VSPs) that we introduced in §2.5, $v_{s1}$ and $v_{s2}$ (see equations 21 and 23).

In Figure 5, we show results for the single component mocks PlumCuspIso, NonPlumCoreOm and PlumCuspTan, but now including constraints from $v_{s1}$ and $v_{s2}$. The top three rows show results using 1000 tracer stars; the bottom row shows what happens when we use 10k tracer stars for the PlumCuspTan mock. The lines and shaded regions are as in Figure 2. We show an example fit for the NonPlumCoreOm mock in Figure 1 (middle row), where the middle and right panels show marginalised histograms of $v_{s1}$ and $v_{s2}$, respectively. The two additional VSP data points that we fit are marked on these histograms, with errors calculated as in §2.7. (For brevity, we omit the similar plots for the PlumCuspIso and PlumCuspTan mocks since GravSphere performs similarly well on these.)

As can be seen from Figure 5, even with just 1000 tracers, including the VSPs produces a dramatic improvement in the performance of GravSphere. The density $\rho(r)$ is well recovered over the range $0.5R_{1/2}<r<2R_{1/2}$, with poorer constraints for $r<0.5R_{1/2}$ and a small bias creeping in at larger radii. This bias is not a particular cause for concern since there are few tracers this far out. For the PlumCuspIso and PlumCuspTan mocks with 1000 tracers, the $\rho-\beta$ degeneracy is only weakly broken. However, for the NonPlumCoreOm mock, the improvement is dramatic. The VSPs are able to detect this mock’s strong radial anisotropy, reproducing the inner isotropic $\beta(0)\sim 0$ and its sharp radially anisotropic rise to large $r$ (see Figure 5, middle row, rightmost panel).

With 10,000 tracers, the VSPs perform even better. Figure 5, bottom row, shows results for the PlumCuspTan mock using 10k tracers (similar results are obtained with 10k tracers for the other mocks and so we omit these for brevity). We now obtain tight constraints on $\rho$ over the range $0.25R_{1/2}<r<4R_{1/2}$ at better than 95% confidence. In particular, this leads to an excellent recovery of $d\ln\rho/d\ln r(r)$ (middle panel, bottom row of Figure 5) and $\tilde{\beta}(r)$ (rightmost panel, bottom row of Figure 5).

Finally, we note that combining VSPs with split populations may yield even tighter constraints. We will explore this idea further in future work.

4.4 Including proper motion data

In this subsection, we consider the utility of adding proper motion data. In projection, this means adding the tangential ($\sigma_{\text{pmt}}$) and radial ($\sigma_{\text{pmr}}$) velocity dispersion constraints, as described in §2.6. Our results for the PlumCuspIso, NonPlumCoreOm and PlumCuspTan mocks are shown in Figure 6. The lines and shaded regions are as in Figure 2. An example model fit to these data is shown in Figure 1 (bottom row) for the NonPlumCoreOm mock (GravSphere performed similarly well on the other mocks and so we omit their similar plots for brevity.)

As can be seen from Figure 6, even with 1000 tracers stars GravSphere gives an excellent recovery of $\rho(r)$, $d\ln\rho/d\ln r$ and $\tilde{\beta}$, once proper motion data are added. As discussed in §2.6, the addition of $\sigma_{\text{pmt}}$ and $\sigma_{\text{pmr}}$ should provide a strong breaking of the $\rho-\beta$ degeneracy and this is reflected in the performance of GravSphere in Figure 6. As was the case when adding VSPs (Figure 5), the NonPlumCoreOm mock becomes slightly biased at large radii. Otherwise, the results for all three mocks are excellent, recovering the input models within our 95% confidence intervals.

While proper motion data clearly gives the best performance, collecting data of the quality assumed in Figure 6 will require a new and dedicated satellite astrometric mission (e.g. Strigari et al., 2007). Even then, we would only be able to study to the very nearest stellar systems. By contrast, both split populations and VSPs can be used on existing data and will continue to improve as spectra for more and more tracers are collected.

4.5 Triaxial Mock Data

Finally, we consider in this subsection how GravSphere performs on triaxial mock data for which we expect it to become biased (because the method fundamentally assumes spherical symmetry; see §2). In Figure 7, we show results for the TriaxCusp mock (see Table 1) looking down the long axis ($X$) and in Figure 8, we show the same looking down the short axis ($Z$). The lines and shaded regions are as in Figure 2. From top to bottom, the panels show results for 1000 tracers with line-of-sight velocities only; including VSPs; and including proper motions. We show the model fits to the data for the short axis case in Figure 9 (the long axis behaves similarly).

Firstly, notice that when staring down the long axis, $\rho(r)$ is recovered very well with no noticeable bias (compare the grey shaded regions and blue dashed lines in the left and middle columns of Figure 7). Only when adding full proper motion data do the models become biased at greater than 95% confidence, and then only in $\tilde{\beta}$ (see the bottom right panel of Figure 7).

The performance of GravSphere is slightly worse when looking along the short axis, however (Figure 8). Now when adding VSPs. the density becomes biased low at all radii, while a similar albeit weaker effect is also seen when adding the proper motion data. Interestingly, however, in both cases we would be able to detect that something has gone wrong. As can be seen in Figure 9, the actual fits to the data are excellent when only line-of-sight velocities are available (top row). However, when adding either VSPs (middle row) or proper motions (bottom row), the fits become visibly poor. In the case of VSPs, the parameter $v_{s1}$ can no longer be fit, while for the proper motions, $\sigma_{\text{LOS}}$ becomes visibly biased high as compared to the data. A corollary of this is that there appears to be power in the data to determine that the system is not spherical. We will explore this further in future work.

Finally, notice that although GravSphere becomes slightly biased when applied to triaxial mock data, in almost all cases the density at $R_{1/2}$ remains well-recovered. This is consistent with recent work by Sanders & Evans (2016) who show that Jeans estimators remain valid for stellar systems that are flattened along constant ellipsoids.

5 Conclusions

We have presented a new non-parametric Jeans method – GravSphere – and tested it on a wide range of mock data. Our goal was to assess how well three popular methods for breaking the $\rho-\beta$ degeneracy work: split populations; ‘Virial Shape Parameters’ (VSPs); and proper motions. Our key results are as follows:

• •

We confirm that with only line-of-sight velocity data, GravSphere provides a good estimate of the density at the projected stellar half mass radius, $\rho(R_{1/2})$, but is not able to measure $\rho(r)$ or $\beta(r)$, even with 10,000 tracer stars.

• •

We find that two populations provide an excellent recovery of $\rho(r)$ in-between their respective $R_{1/2}$. However, even with a total of $\sim 7,000$ tracers, we are not able to well-constrain $\beta(r)$ for either population.

• •

Using 1000 tracers with higher order VSPs, we are able to measure $\rho(r)$ over the range $0.5<r/R_{1/2}<2$ and broadly constrain $\beta(r)$. With 10,000 tracers, we obtain an exquisite measure of $\rho$ and $\beta$ over the range $0.25R_{1/2}<r<4R_{1/2}$ at better than 95% confidence.

• •

Including proper motion data for all stars gives an even better performance than using VSPs. With just 1000 tracers, we can measure $\rho$ and $\beta$ over the range $0.25<r/R_{1/2}<4$. However, obtaining these data will require a specialised astrometric mission. For the foreseeable future, this will only be possible for the very nearest stellar systems.

• •

We tested GravSphere on a triaxial mock galaxy that has axis ratios typical of a merger remnant, $[1:0.8:0.6]$. In this case, with 1000 tracers and using VSPs or proper motions, GravSphere can become slightly biased. However, we find that when this occurs the data are poorly fit, allowing us to detect when such departures from spherical symmetry become problematic. A corollary of this is that there appears to be power in the data to determine that the system is not spherical. We will explore this further in future work.

6 Acknowledgments

JIR would like to acknowledge support from SNF grant PP00P2_128540/1; STFC consolidated grant ST/M000990/1; and the MERAC foundation. We would like to thank the Gaia Challenge initiative for many stimulating workshops that led to the mock data compilation that we have used here. In particular, JIR would like to thank the Nordita institute for hosting the last Gaia Challenge meeting at which these ideas came to fruition. We would like to thank Matt Walker, Jorge Peñarrubia, Mark Gieles, Michelle Collins, Payel Das, Laura Watkins, Gary Mamon, Eugene Vasiliev and John Magorrian for useful feedback and stimulating discussions.

Appendix A The choice of binning

In this Appendix, we explore the sensitivity of our results to our choice of data binning. By default, we use $N_{\rm bin}=\sqrt{N_{*}}$ bins where $N_{*}$ is the total number of stellar tracers in the data (see §2.7). For 1000 tracers, this gives $N_{\rm bin}=32$. In Figure 10, we show how our results for the PlumCuspIso mock change when we use instead $N_{\rm bin}=16$ (top row) and $N_{\rm bin}=64$ (bottom row), assuming only line-of-sight velocities for 1000 tracer stars. As can be seen, these changes in binning produce no noticeable effect on the results.

Appendix B Virial Shape Parameter Estimators

In this work, we estimate $v_{s1}$ and $v_{s2}$ from the data by numerically integrating the observed $\langle v_{\rm LOS}^{4}\rangle$ over our model-fitted $\Sigma_{*}$ (see equations 21 and 23). We then estimate the errors on $v_{s1}$ and $v_{s2}$ using a Monte-Carlo sampling method described in 2.7. In this appendix, we compare this procedure with the estimators given in Richardson & Fairbairn (2014):

[TABLE]

where we recall that $\Sigma_{*}$ is the surface number density of stars and $dN=2\pi\Sigma_{*}RdR$. Similarly:

[TABLE]

We estimate errors on the Richardson & Fairbairn (2014) estimators by Monte-Carlo sampling over the errors on the individual $v_{{\rm LOS},i}$ measurements and calculating the variance in the solutions to equations 34 and 35.

In Figure 11, we compare our estimators for $v_{s1}$ and $v_{s2}$ with those in equations 34 and 35. We plot the estimated over the true VSP (as calculated from the mock density, tracer density and velocity anisotropy profiles), VSP/VSPtrue, as a function of the number of tracer stars, $N_{s}=100,1000$ and $10,000$. The data points are slightly offset from one another horizontally to aid clarity. The horizontal dashed line marks VSP/VSP${}_{\rm true}=1$, corresponding to zero bias. The blue data points show results for $v_{s1}$, the green for $v_{s2}$. The squares show results from the estimators presented in this work; the circles show those from Richardson & Fairbairn (2014) (equations 34 and 35). We show results for two different Gaia Challenge mock data sets: the NonPlumCoreOm mock discussed already in this paper, and a similar mock but with a dark matter cusp, NonPlumCuspOm, as marked on the Figure (see Table 1).

As can be seen in Figure 11, increasing the number of tracer stars, $N_{s}$, reduces the bias for both estimators. As might be expected, there is less bias for $v_{s1}$ (blue) than $v_{s2}$ (green) since $v_{s2}$ is much more sensitive to the fall-off of the light profile at large radii where there is little data. Our new estimators for $v_{s1}$ and $v_{s2}$ perform similarly well to those from Richardson & Fairbairn (2014). As anticipated in §2.5, the main advantage of our new estimators is a more realistic error estimate at low $N_{s}$. Finally, notice that the NonPlumCoreOm mock is less biased, particularly for $v_{s2}$ than the NonPlumCuspOm mock. This is because the NonPlumCuspOm mock has the same dark matter scale length as NonPlumCoreOm but the stars are much more concentrated, with a scale length of $r_{*}=0.1$ kpc as compared to $r_{*}=0.25$ kpc. This means that at fixed $N_{s}$, there is correspondingly less information about the large radius behaviour of the mass distribution that $v_{s2}$ is particularly sensitive to. We will discuss this mock in more detail in a companion paper where we describe the full Gaia Challenge mock data suite (Read et al. in prep.).