Gravitational Potential from small-scale clustering in action space: Application to Gaia DR2

TL;DR

This paper introduces a novel method using clustering in action space to measure the Milky Way's gravitational potential, avoiding common assumptions and providing new insights from Gaia DR2 data.

Contribution

It develops a likelihood-based approach from first principles that leverages small-scale clustering in action space to infer gravitational potential parameters.

Findings

Measured dark matter halo profile parameters from Gaia DR2 data.

Confirmed that maximum likelihood corresponds to maximum clustering in action space.

Derived a Milky Way circular velocity curve consistent with previous measurements.

Abstract

Most measurements of mass in Astronomy that use kinematics of stars or gas rely on assumptions of equilibrium that are often hard to verify. Instead, we develop a novel idea that uses the clustering in action space, as a probe of underlying gravitational potential: the correct potential should maximize small-scale clustering in the action space. We provide a first-principle derivation of likelihood using the two-point correlation function in action space, and test it against simulations of stellar streams. We then apply this method to the 2nd data release of Gaia, and use it to measure the radial force fraction $f_h$ and logarithmic slope $\alpha$ of dark matter halo profile. We investigate stars within 9-11 kpc and 11.5-15 kpc from Galactic centre, and find $(f_h,\alpha)= (0.391\pm 0.009, 1.835\pm 0.092) $ and $(0.351\pm 0.012,1.687\pm 0.079)$, respectively. We also confirm that the set of parameters that maximize the likelihood function do correspond to the most clustering in the action space. The best-fit circular velocity curve for Milky Way potential is consistent with past measurements (although it is $\sim$ 5-10\% lower than previous methods that use masers or globular clusters). Our work provides a clear demonstration of the full statistical power that lies in the full phase space information, relieving the need for {\it ad hoc} assumptions such as virial equilibrium, circular motion, or steam-finding algorithms.

Tables

Table 1. Table 1: Initial conditions of progenitor for the generation of stream stars (where galactocentric radius and velocity are normalized by solar radius value.

Stream Number	Initialization
	R	$\varphi$	z	$V_R$	$V_T$	$V_Z$	${\sigma }_v\left(\mathrm{km}/s\right)$	$t_{\mathrm{disrupt}}$ (Gyr)
Stream 1	1.56	0.12	0.89	0.35	$-1.15$	$-0.48$	0.3	2
Stream 2	1.00	$-0.05$	0.001	$-0.60$	0.51	0.0086	0.3	2
Stream 3	1.20	$-0.05$	$-1$	$-0.30$	0.51	0.16	0.3	2

Table 2. Table 2: Constraints on normalization and logarithmic slope of dark matter profile ( f h , α ) subscript 𝑓 ℎ 𝛼 (f_{h},\alpha) for both radial samples with selection cuts. Here, we summarize all sources of measurement errors: the stochastic errors estimated from the posterior distribution of parameters, the systematic errors from simulation, and the total error given by the root of stochastic errors squared plus systematic errors squared.

Error Type ($\left|z\right|$ > 1   kpc)	[9.0  kpc <R<11.0  kpc]		[11.5  kpc <R<15.0  kpc]
	$f_h$ = 0.375	$\alpha$ = 1.967	$f_h$ = 0.387	$\alpha$ = 1.786
Stochastic	${}_{-0.038}{}^{}{}_{}{}^{+0.029}$	${}_{-0.070}{}^{}{}_{}{}^{+0.090}$	${}_{-0.027}{}^{}{}_{}{}^{+0.024}$	${}_{-0.181}{}^{}{}_{}{}^{+0.070}$
Systematic	$\pm$ 0.004	$\pm$ 0.079	$\pm$ 0.004	$\pm$ 0.071
Total	$\pm$ 0.034	$\pm$ 0.112	$\pm$ 0.026	$\pm$ 0.144

Table 3. Table 3: Constraints on normalization and logarithmic slope of dark matter profile ( f h , α ) subscript 𝑓 ℎ 𝛼 (f_{h},\alpha) for both radial samples without selection cuts. Here, we summarize all source of measurement errors: the stochastic errors estimated from the posterior distribution of parameters, the systematic errors from simulation, and the total error given by the root of stochastic errors squared plus systematic errors squared.

Error Type (no $\left|z\right|$ cut)	[9.0  kpc <R<11.0  kpc]		[11.5  kpc <R<15.0  kpc]
	$f_h$ = 0.392	$\alpha$ = 1.756	$f_h$ = 0.345	$\alpha$ = 1.656
Stochastic	${}_{-0.006}{}^{}{}_{}{}^{+0.004}$	${}_{-0.040}{}^{}{}_{}{}^{+0.030}$	$\pm$ 0.006	$\pm$ 0.040
Systematic	$\pm$ 0.004	$\pm$ 0.070	$\pm$ 0.003	$\pm$ 0.066
Total	$\pm$ 0.006	$\pm$ 0.078	$\pm$ 0.007	$\pm$ 0.077

Table 4. Table 4: Constraints on dark matter halo parameters for both samples with/without selection cuts. The top four rows compare results for different cuts (including only, mostly independent, stochastic errors), while the final two rows is an attempt to combine these results, including the systematic errors introduced due to selection, as detailed in Appendix D .

	[9.0  kpc <R< 11.0  kpc]	[11.5  kpc <R< 15.0  kpc]
$f_h$	${0.375}_{-0.038}^{+0.029}$	${0.387}_{-0.027}^{+0.024}$
$\alpha$ ($\left|z\right|$>1 kpc, stochastic error only)	${1.967}_{-0.070}^{+0.090}$	${1.786}_{-0.181}^{+0.070}$
$f_h$	${0.392}_{-0.006}^{+0.004}$	0.345$\pm$ 0.006
$\alpha$ (no $\left|z\right|$ cut, stochastic error only)	${1.756}_{-0.040}^{+0.030}$	1.656$\pm$ 0.040
$f_h$	0.391$\pm$ 0.009	0.351$\pm$ 0.012
$\alpha$ (combined, stochastic+systematic error)	1.835$\pm$ 0.092	1.687$\pm$ 0.079

Equations

$$\rho_{b}(r)\propto r^{-1.8}\exp\bigg{[}-\bigg{(}\frac{r}{1.9~{}{\rm kpc}}\bigg{)}^{2}\bigg{]}.$$

$${\Phi }_d(r,z)\propto -\frac{1}{\sqrt{R^2+{\left[3\mathrm{kpc}+\sqrt{z^2+(0.28\mathrm{kpc}{)}^2}\right]}^2}},$$

$${\rho }_{dm}(r)\propto r^{-\alpha },$$

$$J_i=\frac{1}{2\pi }{\oint }_{{\gamma }_i}\mathbf{p}\cdot d\mathbf{q},$$

$$R=\Delta \sinh u\sin v;$$

$$\Phi (u,v)=\frac{U(u)-V(v)}{{\sinh }^{2}u+{\sin }^{2}v}.$$

$$J_r$$

$$J_z$$

$$\frac{p_u^2}{2{\Delta }^2}$$

$$\frac{p_v^2}{2{\Delta }^2}$$

$$\frac{\partial^{2}}{\partial u\partial v}\bigg{[}(\sinh^{2}u+\sin^{2}v)\Phi_{\textrm{model}}\bigg{]}\approx 0.$$

$$\Delta^{2}=z^{2}-R^{2}+\bigg{(}3R\frac{\partial\Phi}{\partial z}-3z\frac{\partial\Phi}{\partial R}+Rz\bigg{(}\frac{\partial^{2}\Phi}{\partial R^{2}}-\frac{\partial^{2}\Phi}{\partial z^{2}}\bigg{)}\bigg{)}/\frac{\partial^{2}\Phi}{\partial R\partial z}.$$

$$\frac{\Delta p}{p},\frac{\Delta {\mu }_{\text{RA}}}{|{\mu }_{\text{RA}}|},\frac{\Delta {\mu }_{\text{dec}}}{|{\mu }_{\text{dec}}|},\frac{\Delta V_{\text{radial}}}{|V_{\text{radial}}|}<0.2,$$

$$D=\sqrt{(\Delta J_R/{\sigma }_{J_R}{)}^2+(\Delta J_{\varphi }/{\sigma }_{J_{\varphi }}{)}^2+(\Delta J_z/{\sigma }_{J_z}{)}^2},$$

$$\mathcal{P}(D){|}_{\mathrm{uniform}}dD$$

$$⟹\mathcal{P}(\ln D){|}_{\mathrm{uniform}}d\ln D$$

$$1+\xi (\ln D)\equiv \frac{\mathcal{P}(\ln D)}{\mathcal{P}(\ln D){|}_{\mathrm{uniform}}}=\frac{D_{\text{max}}^3\mathcal{P}(\ln D)}{3D^3{\int }_{-\infty }^{\ln D_{\text{max}}}\mathcal{P}(\ln D^{\prime })d\ln D^{\prime }},$$

$$\ln \mathcal{L}(\mathrm{data}|f_h,\alpha )\simeq {⟨\sum _{\mathrm{pairs}}\ln \left[1+\xi (\ln D_{\mathrm{pair}})\right]⟩}_{\mathrm{pairings}}=$$

$$N_{\mathrm{pairs}}{\int }_{-\infty }^{\ln D_{\text{max}}}\mathcal{P}(\ln D)\ln \left[1+\xi (\ln D)\right]d\ln D.$$

$$S_{\mathrm{relative}}\equiv -\int \mathcal{P}(\ln D)\ln \left[\frac{\mathcal{P}}{\mathcal{P}{|}_{\mathrm{uniform}}}\right]d\ln D=-\frac{\ln \mathcal{L}(\mathrm{data}|f_h,\alpha )}{N_{\mathrm{pairs}}},$$

$${\chi }^2={\chi }_{min}^2+(X_i-\overline{X_i})F_{ij}(X_i-\overline{X_i}{)}^T,$$

$$0.5<$$

$$0.25<$$

$$\mathcal{P}(\Theta |\mathcal{D})=\frac{\mathcal{P}(\mathcal{D}|\Theta )\mathcal{P}(\Theta )}{\mathcal{P}(\mathcal{D})}$$

$$\mathcal{P}(\{{\nu }_a\},\{{\chi }_a\})=\frac{\exp \left(-\frac{1}{2}{\sum }_{a,b}{\chi }_a{\xi }_{ab}^{-1}{\chi }_b\right)}{\sqrt{\text{det}(\xi )}(2\pi {)}^{M/2}}\prod _a\frac{(\bar{n}+\bar{n}{\chi }_a{)}^{{\nu }_a}\exp (-\bar{n}-\bar{n}{\chi }_a)}{{\nu }_a!}.$$

$$⟨{\chi }_a{\chi }_b⟩={\xi }_{ab},$$

$$\mathcal{P}(f_h,\alpha |\{a_k\})\propto {\mathcal{P}}_{\mathrm{prior}}(f_h,\alpha )\int \frac{\exp \left(-\frac{1}{2}{\sum }_{a,b}{\chi }_a{\xi }_{ab}^{-1}{\chi }_b\right)}{\sqrt{\text{det}(\xi )}(2\pi {)}^{M/2}}\left\{\prod _k\left[1+{\chi }_{a_k}\right]\right\}\prod _{a}d{\chi }_a,$$

$$\mathcal{L}(\mathrm{data}|f_h,\alpha )\equiv \frac{\mathcal{P}(f_h,\alpha |\{a_k\})}{{\mathcal{P}}_{\mathrm{prior}}(f_h,\alpha )}\propto ⟨\prod _k\left[1+{\chi }_{a_k}\right]⟩\propto \sum _{\mathrm{pairings}}\prod _{\mathrm{pairs}}\left[1+{\xi }_{\mathrm{pair}}\right]\propto {⟨\exp \left[\sum _{\mathrm{pairs}}\ln \left(1+{\xi }_{\mathrm{pair}}\right)\right]⟩}_{\mathrm{pairings}}.$$

$$\ln \mathcal{L}(\mathrm{data}|f_h,\alpha )\approx {⟨\sum _{\mathrm{pairs}}\ln \left(1+{\xi }_{\mathrm{pair}}\right)⟩}_{\mathrm{pairings}}$$

$$\mathcal{L}(x,{\sigma }_{\mathrm{sys}}|\{x_i,{\sigma }_i\})=\prod _i\frac{\exp \left[-\frac{(x-x_i{)}^2}{2({\sigma }_i^2+{\sigma }_{\mathrm{sys}}^2)}\right]}{\sqrt{2\pi ({\sigma }_i^2+{\sigma }_{\mathrm{sys}}^2)}}.$$

$$P(x)\propto {\int }_0^{{\sigma }_{\mathrm{sys},\max }}d{\sigma }_{\mathrm{sys}}\prod _i\frac{\exp \left[-\frac{(x-x_i{)}^2}{2({\sigma }_i^2+{\sigma }_{\mathrm{sys}}^2)}\right]}{\sqrt{2\pi ({\sigma }_i^2+{\sigma }_{\mathrm{sys}}^2)}}.$$

Full text

††pubyear: 2019††pagerange: Gravitational Potential from small-scale clustering in action space: Application to Gaia DR2–LABEL:lastpage

1 Introduction

Understanding the nature of dark matter is one of the most significant challenges in the 21st century for both physicists and astrophysicists. While the cold dark matter (CDM) paradigm the most popular model of dark matter, observational tensions on the galactic scale (Bullock & Boylan-Kolchin, 2017; Del Popolo & Le Delliou, 2017) and non-detection of dark matter particles in the ground-based experiments have led scientists to examine alternative possibilities. Part of the reason why it is difficult to probe the nature of dark matter is that (as far as we can tell) it only interacts with ordinary matter via gravity; only affecting astrophysical observations on very large scales ($\gtrsim$ few kpc). On these scales, extracting the full 6-dimensional phase space information (that is necessary to infer dark matter mass unambiguously), has been a difficult task.

Nonetheless, European Space Agency (ESA)’s Gaia mission has recently started probing the kinematics of the Milky Way stars with unprecedented precision (Gaia Collaboration et al., 2016). Gaia is a space-based observatory launched by European Space Agency (ESA) in December 2013, which aims at constructing the largest catalogue of 3d positions and velocities of the Milky Way stars, using Astrometric techniques. The first data set (Gaia DR1), released in September 2016, did not have the measurement of radial velocities. However, the Data Release 2 (Gaia DR2), which was released in 2018, with more complete magnitude measurement and longer span compared to Gaia DR1, contains the proper motion, parallax as well as the radial velocity information of more than 7 million stars (Gaia Collaboration et al., 2018). The proliferation of data on the kinematics of Milky Way’s stars has therefore opened a new avenue to probe the structure of the Milky Way potential and the nature of dark matter.

In this work, we propose a new method to constrain the potential of the Milky Way and apply it to two six dimensional sub-samples of stars in Gaia DR2. Our method is based on maximizing the statistical clustering of the stars in the space of actions. The current theories and observational evidence suggest that the growth of structure in our universe is hierarchical, where smaller structures merge to form bigger ones. During the formation of galaxies, however, the smaller structures are tidally disrupted and due to various relaxation mechanisms at play, the memory of their common origin in configuration space is erased. This makes identifying stars with common origin nearly impossible. Nevertheless, the information regarding their common origin may still be present in the phase space of action variables. The action variables will remain conserved as long as the host potential evolves adiabatically. When the smaller structures are tidally disrupted in the Milky Way potential, their spread in the action space is much smaller than the rest of the stars in the Milky Way. Therefore, we expect the small scale structure of the action space to contain the hierarchical tidal disruption/assembly history of the Milky Way (Afshordi et al., 2009).

Since action variables are conserved, due to their common origin, the action variables of the various structures would be clustered on small scales (in action space). This principle has previously been proposed to infer the potential by maximizing Kullback-Liebler divergence (KLD) or "relative entropy" in the action space (Kullback & Leibler, 1951; Sanderson et al., 2015, 2017). In these two studies, the viability of this idea was tested in simulated stellar distribution, where they successfully recover a spherical isochrone potential (Sanderson et al., 2015) and a spherical NFW profile (Sanderson et al., 2017) using this method. The same method has also been applied to constrain the parameters of the globular cluster. By minimizing the KL entropy in the phase-space, Buckley et al. (2019) successfully constrain the mass and the King radius of the simulated M4 globular cluster, which provides another proof that the original parameters of a system can be recovered if the true phase-space information are found. Maximizing the statistical clustering in the action space does not require identification of the membership of any star, which is one of the major merits for this method. Related methods have been proposed by Peñarrubia et al. (2012) and Magorrian (2014). However, neither have applied 2-point correlation function as a measure of clustering in the action space. Therefore, the significance of our work is to demonstrate the viability of using 2-point correlation function to measure the mass of a system.

This principle can be used to infer the potential of the Milky way. If the action variables are estimated using the incorrect potential, the resulting quantity will not be conserved with the dynamical evolution. Therefore, the clustering of the stars, in action space, on small scales will be destroyed if we use the wrong potential. Conversely, using the correct potential will maximize the small scale clustering in this space. We provide a first-principle derivation that the likelihood for the potential can be expressed as an integral (or KLD) over the 2-point correlation function in the action space, and test it using simulations of mock streams. Then, as an example, we fit a power-law dark matter profile (assuming a fixed form for bulge and disk component) to Gaia DR2, and compare our results to those that use other methods.

The paper is structured as follows: In Section 2, we introduce the required theoretical background for our methodology, including our parametrized models of the Milky Way potential and the computation of action variables. In Section 3, we briefly discuss the data sets we used and the selection cuts imposed on the raw data. Details related to the two-point correlation function and the likelihood test we used in the action space are presented in Section 4. To check the viability of our method, we first apply it to simulations with only stream stars, where the streams are simulated using the Python package galpy (Version 1.3.0. See Bovy, 2015, for more details). Then we proceed to apply our method to the real observations taken from Gaia DR2 with the selection cuts listed in Section 3 from two different radial bins. The results of our analysis is presented in Section 5. In Section 6, we discuss the shortcomings and future possibilities of our method before concluding in Section 7.

2 Theory

2.1 Modelling the Milky Way Potential

We use a parametrized model for the Milky Way potential which can be regarded as a combination of the bulge, the disk and the dark matter halo. The specific model we are using is a slight modification to the MWPotential2014 potential in galpy, which is assumed to be a good approximation to the Milky Way potential.

The bulge is modelled with a power-law density with an exponential cut-off:

[TABLE]

The contribution from the central bulge is negligible at radius greater than 9 kpc, but we still include this component in the model for completeness.

The disk is modelled as a Miyamoto-Nagai Potential profile, but with fixed parameters (Bovy, 2015):

[TABLE]

where $R$ and $z$ are the radial and vertical galactocentric cylindrical coordinates, respectively ($z=0$ is the plane of the galaxy).

Finally, we model the dark-matter halo profile as a spherical power-law profile

[TABLE]

which is also a built-in potential expression categorized as PowerSphericalPotential in galpy package. Note in MWPotential2014, the halo potential is characterized as NFW profile (NFW: Navarro et al., 1996). Instead, in this work, we use a power law potential for simplicity. Also, as we shall see later, the data set for constraining the potential does not span over a large range, so a localized power law potential should be a good approximation for NFW profile.

Finally, we set the normalizations of bulge, disk, and dark matter components, $\rho_{b},\Phi_{d},\rho_{dm}$, so that the fraction of radial force due to dark matter is $f_{h}$ at the position of the Sun, while the ratio of force due to stellar bulge to disk is fixed to be 1:12. This gives the normalization of halo component as 0.65$f_{h}$/(1-$f_{h}$), and the normalization of the bulge and disk component are fixed to be 0.05 and 0.60 respectively. Therefore, in the end we are left with two free parameters, $f_{h}$, and the power-law index of the density profile, $\alpha$, which we aim to constrain using our method.

2.2 Action Variables

Regular (i.e. non-chaotic) orbits in the galactic potential should admit 3 integrals of motion (e.g., Mo et al., 2010). However, finding these integrals of motions in terms of the phase space coordinates could be a difficult task. Nevertheless, it may be possible to find canonical transformations so that, finding integrals of motions in these coordinates are easy. One particularly convenient system of canonical variables is the so-called action-angle variables [denoted by (${\bm{\theta}}$, ${\bm{J}}$)], where the canonical momenta $\bm{J}$, or actions, are also the integrals of motion. The angle variables $\bm{\theta}$ are periodic in orbital torus, where an increase of $2\pi$ in the angle would be associated with the same point in phase space. The action conjugate to this angle is then defined as:

[TABLE]

where $\gamma_{i}$ is the orbit section where the $i$-th angle, $\theta_{i}$ increases from [math] to $2\pi$.

Finding a canonical transformation to transform to the action-angle variables provide a convenient convention to define integrals of motion for an integrable potential.

2.2.1 Calculating the action variables

The action-angle variables provide a convenient way to find the integrals of motion. However, finding closed analytic forms for the actions is only possible for a few special potentials. Therefore, in galactic dynamics, we often have to rely on approximate methods which involve integration of the orbits, e.g., adiabatic approximation (Binney, 2010) or the torus construction method (Binney & McMillan, 2016). Here, we use the ‘Stäckel approximation’ (Binney, 2012), which is implemented in galpy.

Stäckel potentials are a special class of potentials where the Hamiltonian can be written a separable form, using a canonical transformation. The Stäckel potentials are expressed in the spheroidal coordinates, $(u,v)$ which are related to cylindrical coordinates through

[TABLE]

In these coordinates, the Stäckel potential takes the form

[TABLE]

The radial and the azimuthal actions can be expressed in closed analytic forms

[TABLE]

where,

[TABLE]

In the above relations, $E$ is the energy of the orbit and $I_{3}$ is a third integral of motion (apart from the energy and the azimuthal angular momentum) which can be expressed analytically in terms of the conjugate variables and the momenta.

The Stäckel approximation involves approximating any nominal potential of the Milky Way as a Stäckel potential, and use this potential to find approximate action variables. To do this, we need to find the effective confocal length, $\Delta$ in Equation 5.

A prescription for this was given by Sanders (2012). To obtain the value of $\Delta$, we make use of the fact that $\Phi(\sinh^{2}u+\sin^{2}v)$ is a separable function of $u$ and $v$. Therefore, any mixed derivative of this quantity must vanish. Using the model potential, $\Phi_{\textrm{model}}$, we obtain

[TABLE]

This equation can then be solved for $\Delta$ in terms of the derivatives of the potential, $\partial\Phi/\partial R$, $\partial\Phi/\partial z$, $\partial^{2}\Phi/\partial R^{2}$, $\partial^{2}\Phi/\partial z^{2}$ and $\partial^{2}\Phi/\partial R\partial z$. The relation between $\Delta$ and the derivatives were obtained in Sanders (2012):

[TABLE]

This algorithm is implemented using the python package galpy.

3 Data Set

The kinematics of stars in this paper are obtained from the Data Release 2 (DR2) of the Gaia mission, which was released in April 2018. All data can be accessed though Gaia Archive111Gaia Archive website: http://gea.esac.esa.int/archive/. With the aid of the Gaia radial velocity spectrometer (Cropper et al., 2018), we can directly obtain the full six-dimensional phase space information of approximately 7 million stars. For the sample we used for our analysis, we make a few simple quality cuts in order to avoid the stars with large errors on the parallax or proper motion measurement. We shall impose the cut that the relative error on the parallax, and the proper motion are less than $20\%$:

[TABLE]

where, $p,{\mu}_{RA},{\mu}_{\textrm{dec}}$, $V_{\textrm{radial}}$ are the parallax, proper motion in the right ascension direction, the proper motion in the declination direction, and radial velocity respectively. Here, $\Delta$’ denotes the measurement error in each of these quantities. After applying the selection cuts mentioned above to the raw data, we are left with around 5.6 million stars in the data sample for further analysis222Note that for the actual data analysis, there are additional cuts for halo vs. all stars, and radial distribution, which will be discussed in Section 5.2.. Even though the $20\%$ cut (while common as in Sanderson et al. (2015, 2017), also in Astraatmadja & Bailer-Jones (2016); Schönrich & Aumer (2017); Schönrich et al. (2019)) is ad hoc, we further verify that this choice has little effect on our results, as stars with large uncertainties in their phase space coordinates are unlikely to form close pairs in action space333We did change the choice of relative error cut to another value, for instance, to 10%. However, this does not affect our final estimations..

This data allow us to constrain the two-parameter power law potential in Equation 3 by using the likelihood test, which we will discuss in the next section. As calculations are conducted by using the built-in functions in galpy, the inputs for most of the functions we used are in cylindrical coordinates. For the convenience of computation, all calculations are done in galactocentric coordinate system, and the coordinate transformation are as well handled by the built-in functions in galpy library (Bovy, 2015).

4 Method

We are interested in the small scale clustering of the stars in the action space. There are many different measures of clustering which are useful for different purposes. Here, we use two-point correlation function as our measure of clustering, which should be one of the most straightforward ones. However, to define the correlation function, we need to have a measure of the distance. While this choice is not unique, we shall use the following measure to find the distance of two stars in action space

[TABLE]

where $\Delta J_{i}$ denotes the difference in the action coordinates of the two stars, while $\sigma_{J_{i}}$’s are standard deviations of $J_{i}$’s over all stars. The reason why we normalize the difference in action by the standard deviation in the action space is that stream stars, due to their common origin, should be significantly more clustered than the background, i.e. $\Delta J\ll\sigma_{J}$. Since the background could be anisotropic in the action space, this normalization provides a more appropriate distance measure. We further discuss this choice in Section 6 below.

The calculation of $\sigma_{J_{i}}$ can be affected by the outliers in the raw data or numerical artifacts in galpy action approximation. Therefore, another constraint is added to effectively exclude outliers out of the sample with $\frac{|J_{i}-\bar{J_{i}}|}{\sigma_{J_{i}}}>3$. This choice of 3$\sigma$ seems relatively arbitrary. However, we confirm that the constraints are not affected by this choice, as long as those outliers are safely removed. Also, there is only a small fraction of stars being cut off by the “3$\sigma$-cut” from the original catalogue. Therefore, we do not believe our results are significantly biased by this choice. We then use the remaining action variables that satisfy the above criterion to re-calculate the standard deviation.

For points distributed randomly with a uniform distribution in a three dimensional action space, the probability of finding pairs at a separation between $D$ and $D+dD$ is given by

[TABLE]

However, the actual probability distribution $\mathcal{P}(\ln D)$ will be different from $\mathcal{P}(\ln D)|_{\rm uniform}\propto D^{3}$ due to clustering in the action space. This clustering can be quantified using the 2-point correlation function $\xi(\ln D)$:

[TABLE]

where we used the fact that both $\mathcal{P}$ and $\mathcal{P}|_{\rm uniform}$ should integrate to unity over the range $\ln D\in(-\infty,\ln D_{\textrm{max}})$.

The blue solid line in Figure 1 shows $\ln\left[\frac{P(\ln D)}{D^{3}}\right]$ as a function of $\ln(D)$ calculated by using the Gaia DR2 real data from galactocentric radius 11.5-15 kpc with [$f_{h}$ = 0.34, $\alpha$ = 1.66]444As shall be seen in Section 5.2, this corresponds to the best-fit halo potential recovered by our method, where $D$ is the normalized distance of pairs of stars in the action space (as shown in Equation 14), while $\mathcal{P}(\ln D)$ is its probability density over all pairs of Gaia DR 2 within our sample. According to Equation 4, if stars are uniformly distributed, a plateau is expected at small values of $\ln D$, and deviations from this plateau would correspond to clustering. As can be seen from this figure, at small distance, $\mathcal{P}(\ln D)$ roughly obeys uniform distribution as stated in Equation 4, and the probability drops significantly due to the lack of stellar pairs at large values of $\ln D$.

As it turns out, with certain assumptions, the statistical likelihood of any action-space distribution can be expressed in terms of $\xi(\ln D)$. The key idea here is to assume the star distribution in the action space is the Poisson sampling of a near-uniform background plus a random gaussian field. The correlation function of this random gaussian field encodes all the clustering information at small scale in the action space. This model is agnostic about the distribution function, $f({\bf J})$ and instead relates the likelihood to the correlation function in the action-space $\xi(\ln D)$, after marginalizing over all possible $f({\bf J})$’s. More explicitly, we find that the log-likelihood for a potential is given by

[TABLE]

A detailed derivation of this expression is presented in Appendix A. In this equation, $N_{\rm pairs}$ is the total number of pairs, i.e. half of the number of stars in the sample. Also, as shown in Figure 1, $\mathcal{P}(\ln D)$ obeys the scaling of uniform distribution only at small values of $\ln D$. Therefore, when evaluating the value of likelihood function, integration is terminated at a chosen value of $D_{\textrm{max}}$. This $D_{\rm max}$ characterizes the scale of homogeneity in the action space background, and we shall discuss the choice of $D_{\rm max}$ in later section.

We further note that relative entropy (Kullback & Leibler, 1951) of the distribution $\mathcal{P}(\ln D)$, with respect to the uniform distribution, is defined as

[TABLE]

i.e. the maximization of the likelihood function corresponds to minimizing the entropy relative to the uniform pair distribution. In other words, the best-fit values for the dark matter halo density produce the most non-uniform distribution of pairs in the action space. We should note that while this is similar to the criterion proposed by Sanderson et al. (2015), their relative entropy is based on phase space density in the action space $f({\bf J})$, while our derivation in Appendix A shows that likelihood depends on the relative entropy of the pair distance probability distribution $\mathcal{P}(\ln D)$.

5 Results

5.1 Simulations

To validate our method, we simulate the orbits of a few stars in a known parametrized potential of the same form and then applied the above mentioned analysis to check if we can recover the true parameters of the potential. The simulation includes three groups of tidal stream stars with different initial conditions of progenitors. This is achieved by using the built-in modelling method in galpy package (Bovy, 2014). One can specify the gravitational potential that stars evolve in, the method for action variables calculation, the initial conditions of progenitors, the velocity distribution of progenitors and the time when the disruption began. Initial conditions of the progenitors’ orbit for three streams are tabulated in Table 1. The header of the table is organized in the order of R, $\phi$, z, $v_{R}$, $v_{T}$, $v_{z}$, velocity dispersion ($\sigma_{v}$) and the disruption time ($t_{\rm disrupt}$). Each of the three streams consist of 3000 stream stars. With the simulated stellar trajectories, the action variables of stars can be calculated based on the Stäckel approximation as explained in Section 2. We now wish to test whether our proposed likelihood function (Equation 4) leads to constraints that are consistent with parameters that are used in our simulated host potential.

There are two free parameters in the expression of the dark matter halo density profile (Equation 3), $f_{h}$ which fixes the normalization, and the logarithmic slope $\alpha$. We choose the mass fraction of the halo $f_{h}=0.35$ based on Table 1 in Bovy (2015), and we produce two sets of stream simulations with different choices of $\alpha=$ 1.70 and 2.00 respectively. The progenitor stars are evolved in these two host gravitational potentials respectively. We then compute the action variables on a grid in the ($f_{h}$, $\alpha$) space, and compute the corresponding likelihood function using Equation 4. The likelihood functions evaluated with 9000 simulated stream stars for both potentials are shown in Figure 2 (assuming $\ln D_{\textrm{max}}\sim-1$). For each case, we find clear constraints on both parameters as expected, in reasonable agreement with input parameters of the simulations, subject to caveats that we discuss next.

To determine the location and the uncertainties of the measurements at each $D_{\textrm{max}}$, we fit the log-likelihood distribution with a quadratic function around its maximum. The assumption made by this procedure is that the likelihood only has a single peak that can be approximated by a gaussian distribution. In order to check the validity of this assumption, we plot the posterior of parameter at different values of $D_{\textrm{max}}$, which is shown in Appendix B. As can be seen from Figure 10, both of the posterior distributions for $f_{h}$ and $\alpha$ have a single peak that can be reasonably approximated by gaussian. Also, the probability distributions do not drastically vary with the choices of $D_{\textrm{max}}$. Therefore, we conclude that, at least for our simulated streams, our likelihood distribution is well approximated by gaussian statistics:

[TABLE]

where $\chi_{min}^{2}$ is given by the likelihood peak value by assuming $\mathcal{L}\propto e^{-\frac{\chi^{2}}{2}}$, and $F_{ij}$ represents the Fisher information matrix for $X_{1}=f_{h}$ and $X_{2}=\alpha$. The covariance of the parameters is then given by the inverse of the Fisher matrix, $F_{ij}^{-1}$. There are six parameters to be determined in Equation 19. In practice, we fit for these parameters using a $3\times 3$ grid around the peak of the likelihood.

As can be seen from Figures 2, the final measurements of mass fraction and index do not significantly vary with the choices of $D_{\textrm{max}}$, with small error bars that further shrink by increasing the number of stars or $D_{\textrm{max}}$. Based on these plots, we conclude that the measurement error for the simulations are dominated by systematic error, which is at the level of $1\%$ for mass fraction $f_{h}$ and $4\%$ for the power law index $\alpha$. This systematic error, while small, might arise due to the use of the Stäckel approximation to compute our action variables. We need to keep this in mind when we apply our method to real data, highlighting where further improvements may be needed when other sources of error are small.

As a further sanity check, we can verify that the parameters of potential found by the maximum likelihood test do correspond to the most clustering in the action space. In other words, the two-point correlation function defined in Equation 4 should be maximized for the correct potential. To verify this, we show how a 2D projection of stellar distribution in the action space varies with different choices of potential for both simulations. Figures are shown in Appendix B. As expected, the most compact distributions occur when parameters approach the correct values for the simulation, which is also reflected in the behavior of the two-point correlation function.

To summarize, we have confirmed that our method to maximize likelihood (Equation 4) based on clustering in action space can yield reasonable constraints on simulated potentials, subject to small systematic errors of $1\%$ ($4\%$) on normalization and logarithmic slope. There also does not seem to be any significant dependence on the maximum separation of included pairs in action space $D_{\rm max}$. A more exciting step is to apply our method to real Gaia DR2 data to see how well it can constrain the Milky Way potential, which we shall do next.

5.2 Real data

After confirming the reliability of the method, we proceed with our analysis using real data from Gaia DR2. The criteria for data selection were already discussed in Section 3. Let us now introduce some additional selection cuts. Recall that in the derivation of likelihood function (Appendix A), we assume the stellar distribution in the action space is a uniform background plus fluctuations. This assumption is more appropriate for halo stars in our galaxy, as disk stars have $J_{z}\simeq 0$. Additionally, as we are trying to constrain dark matter profile, which mostly occupies the Milky Way halo, halo stars should be better candidates compared to disk stars. Due to these considerations, we only select stars that have vertical distance to the galactic plane $>$ 1 kpc555We will discuss the effect of this cut, as well as the measurements error cuts on the final results later in Section 6..

Figure 3 shows the galactocentric distance and the tangential velocity distribution (in cylindrical coordinate) for all the data with relative measurement error smaller than 20$\%$ and $|z|$ > 1 kpc. There are around 337,022 stars in total. As expected, the peak of radial distribution is around solar radius and the peak of $V_{T}$ distribution is around the value of circular velocity at solar radius. Here, we assume $R_{\odot}$ = 8.122 kpc, the vertical distance to the galactic plane $z_{\odot}$ = 0.025 kpc, and the galactocentric velocity of the Sun $V_{x,\odot}$ = -11.1 km/s, $V_{y,\odot}$ = 245.8 km/s, $V_{z,\odot}$ = 7.8 km/s as taken in Eilers et al. (2019), but we shall discuss this choice further in Sec. 6. Due to the limitation of the computational time, another galactocentric radius cut is also applied to the data: We choose two different radial ranges 9 kpc < R < 11 kpc (hereafter real-data-9-11) and 11.5 kpc < R < 15 kpc (hereafter real-data-115-15). After applying all of these cuts, there are approximately 61,000 and 16,000 stars in each sample, respectively.

Now, taking the NFW profile as reference, the expected value of the $\alpha$ in the power law density profile should be within 1 to 3: For $r\gg r_{s}$, the density is proportional to $r^{-3}$, while for $r\ll r_{s}$, it goes to $r^{-1}$. Furthermore, Bovy & Rix (2013) used the assumption of Jean’s equilibrium for G-dwarfs from SEGUE survey to constrain $\alpha<1.53$ (at 95% confidence) between R=4 kpc and 9 kpc. Therefore, to allow for a conservative prior, we consider the range:

[TABLE]

for our dataset within 9 kpc to 15 kpc.

The likelihood plots showing the constraints on the mass fraction and the index for both radial samples are shown in the left panels of Figure 4. Although calculated within the same $\ln D_{\textrm{max}}\sim-1$, the log-likelihood values for “real-data-9-11” is larger than those of “real-data-115-15” as there are more stellar pairs included in the more nearby sample. Furthermore, the error bar plots in Figure 4 show how the likelihood peaks (black solid points) and median constraints on parameters (blue hollow circles with error bars) vary with different choices of $D_{\textrm{max}}$. We see that for $\ln D_{\textrm{max}}\lesssim-1$, the $f_{h}$ and $\alpha$ constraints are stable and robust to the choice of free parameter $D_{\textrm{max}}$.

For the choice of $\ln D_{\rm max}$, we notice that for the more distant sample “real-data-115-15”, there is a jump in the error bar plots for both parameters when ln $D_{\textrm{max}}$ is smaller than -1.5. Taking this into account, we treat -1.14 as our final choice of ln $D_{\textrm{max}}$. Results estimated at this point have the smallest uncertainties and the constraints are consistent (within error bars) for all ln $D_{\textrm{max}}\lesssim$ -1.14. We shall refer to this value as ln $D_{\textrm{max, optimum}}$. To be consistent, we use the same value of $D_{\rm max}$ for both Gaia samples. More discussions about choosing an appropriate $\ln D_{\textrm{max}}$ are presented in Section 6 and Appendix C.

Having seen statistical constraints on both parameters from the likelihood plots, we would like to evaluate the true uncertainties of the measurements. However, the determination of uncertainties is more subtle compared with the simulations. Unlike simulations, where we found the posterior distribution of parameters had a sharp gaussian peak, we notice that the likelihood 2D plots have multiple peaks for real Gaia data. This can be seen more clearly in the 1D posterior distributions in Figure 5, where (depending on the choice of $D_{\rm max}$) there can be multiple peaks. As a result, it is no longer appropriate to simply assume the likelihood distribution is approximated by a gaussian. In particular, the jump in $f_{h}$ around $\ln D_{\rm max}\simeq$ -1.5 in “real-data-115-15” sample is due to the change in relative heights of the two main peaks in posteriors shown in the top panel of Figure 5. Therefore, we calculated the median of the parameters using the full posterior distribution within our prior range (Equation 5.2), as it is a more robust statistical estimator than average whenever multiple peaks or outliers are presented in the distribution. The $68\%$ confidence interval (68% CI) around the median, which can be also computed from the posterior distribution, is treated as the error on the parameter.

As can be seen in the error bar plots in Figure 4, the maxima of likelihood (black points) are all consistent with the median values within 68% CI (blue hollow points with blue error bars), and the error bars become larger for smaller $D_{\textrm{max}}$, where fewer stellar pairs are included. When choosing ln($D_{\textrm{max}}$) as -1.14, the constraints we get under both situations are summarized in Table 2.

Furthermore, recall that from simulated measurements in Section 5.1, we do expect an additional $4\%$ percent systematic discrepancy for index measurements and an $1\%$ off in mass fraction measurement. We do include these estimates in Table 2 as systematic errors, which can be combined with our stochastic errors to obtain the total expected uncertainties.

Let us now perform the same consistency checks we did in Section 5.1 for simulated data, and see how stars from real data are distributed in the action space. Figures 6 and 7 shows how (2D projections of) the stellar distribution in the action space, as well as its two-point correlation function change as we vary $f_{h}$ or $\alpha$ in the Milky Way halo potential. As expected, the correlation function values for both radial ranges are maximized while the potential approaches the parameters that maximize the likelihood. However, unlike in simulations, instead of a compact cluster, stars in the action space present a more “lath-shaped” distribution, where most stars concentrate around $J_{R}\sim$ 0 for the best-fit potential. In this situation, stars are extended along $J_{\phi}$ axis with no distinct $J_{R},J_{z}$ contributions, which indicates the property of circular (or disk) motion for most stars. This is not surprising as can be seen from the tangential velocity distribution in Figure 3; the disk stars are still the dominant component in our real data samples even though we performed a $|z|$ > 1 kpc cut666Note that the inclusion of disk stars does not contradict with our assumption made in the likelihood. Disk is a mixture of correlated structures and uniform background, so the disk component can also contribute to the small-scale clustering signal. Also, our distance metric defined above has already included the effect of disk component (Equation 14).. Just as in the simulations, stars within the potential that maximizes the small-scale clustering statistics (Figure 7) present the most compact distribution in the action space. The fact that stellar distribution reduces to circular motion for the best-fit potential is in practice consistent with the traditional assumption of circular motion for disk stars, in order to estimate the mass of Milky Way galaxy. However, our method does not explicitly make this assumption, and thus can account for deviations from circular motion, effectively combining (thin+thick) disk+halo stars.

6 Discussion

In the previous analysis, we used some measurement error cuts and a vertical distance cut to real data. However, selection cuts to the raw data could cause unexpected biases in the measured parameters. To investigate the degree to which our results are sensitive to an arbitrary choice of $z$-cut, we randomly choose 90,000 stars from 9-11 kpc sample777There are 607,257 stars in total, but we only choose a subset of the catalogue due to the limitation of computational time. and take all data from 11.5-15 kpc ($\sim$ 90,000 stars in total ), without imposing any of the previous error or distance cuts 888However, we did apply some minimal cuts to the raw data in order to get rid of the unreliable observations, including $|z|$ ¡ 10 kpc and the absolute values of all three components of velocity in cylindrical coordinates are smaller than 500 km/s. The error bar plots are shown in Figure 8. For comparison, we also overplot the results obtained before using the sample with selection cuts. Generally, at same value of $D_{\textrm{max}}$, the uncertainties on parameters are significantly reduced when using the data samples without selection cuts. To be consistent, for both radial ranges, we still take ln $D_{\textrm{max}}$ $\simeq$ -1.14 and check the corresponding constraints on $f_{h}$ and $\alpha$. The results are tabulated in Table 3.

Compared the constraints at the same $D_{\textrm{max}}$ obtained previously but with error selection cuts with the results from full data set, we find a consistency in $f_{h}$ constraint for both radial ranges. However, for $\alpha$, we notice a 11$\%$ systematic discrepancy within 9-11 kpc, and an 7$\%$ discrepancy for the 11.5-15 kpc. For both radial ranges, the index estimates for the uncut sample are lower than those of the cut sample.

This systematic shift is primarily due to the selection cut to the vertical distance, $z$. Although some measurement error cuts are also imposed on the raw data, $z$ distance cut seems to be the most severe: 80$\%$ of raw data survives the measurement error cuts, while only around 6$\%$ remain after the $|z|$>1 kpc cut is imposed. One possible reason for this systematic difference could be the inaccuracy of the simple analytic model for the disk potential used in Equation 2. It remains to be seen whether a more realistic model (e.g. using other datasets), or including the disk parameters in the likelihood marginalization, could lead to more consistent (and realistic) estimates. In order to account for additional potential systematic errors due selection cuts, we use the probability function defined in Appendix D, which yields our final constraints on the mass fraction of dark matter and the index in the localized density profile in Table 4.

After obtaining the constraints on both parameters in the dark matter halo density profile, we can translate them to less model-dependent constraints by computing the rotation curve (circular Keplerian velocity) of the Milky Way, as a function of distance from the centre. This result can then be compared to other studies that use different parametrizations and methods. To obtain a more robust estimation to the circular rotation curve, we evaluate the average and standard deviation of $v_{\rm circ.}(R)\equiv\sqrt{\frac{\partial\Phi_{\rm tot}(R,z=0)}{\partial\ln R}}$ given the likelihoods found from our different Gaia samples (Equation 4) over our prior range of $f_{h}$ and $\alpha$ (Equation 5.2).

As the expression of disk potential is fixed (where we also fix $z$ = 0 in Equation 2), its contribution to the total rotation curve (as well as that of the bulge) can be simply added to the halo part in quadrature. Right panel of Figure 9 displays the circular velocity curve obtained from this work (purple and black solid line) and its corresponding uncertainty shaded area. For comparison, the results obtained from Bovy (2015), McMillan (2017), Vasiliev (2019), and Eilers et al. (2019) are also shown in the same figure. McMillan (2017) used kinematic data from maser observations with (expected) near-circular motion to fit a Milky Way model with an NFW spherical halo, a stellar and gas disk plus a central bulge. Using a nearly identical model, Vasiliev (2019) assumed Jeans equilibrium of Milky Way globular clusters in Gaia data to constrain the gravitational potential. Eilers et al. (2019) also used Jeans equilibrium for Gaia luminous red-giant stars to determine the circular velocity of the Milky way over radial range 5 kpc < R < 25 kpc. Although we approximate the localized halo density profile as a simple power law, our result is consistent with the estimation from Eilers et al. (2019) and relatively close to (but around 4% higher than) the best-fit NFW dark matter potential found in Bovy (2015) (MWPotential2014). However, the circular velocity (radial force) is about 5-10% (10-17%) smaller than the other two studies. Compared with these studies, our method might be more robust as it does not rely on assumptions of circular motion or Jeans equilibrium, and can be equally applied to halo or disk stars.

An important consideration for comparison to other measurements of circular velocity is our choices of the solar coordinates in the coordinate transformation. If we make the choice of solar coordinates consistent with the analysis of Bovy (2015) ($V_{\phi,\odot}=220$ km/s and $R_{\odot}=8$ kpc), as shown in the left panel of Figure 9, the recovered circular velocity curve is comparable to Bovy (2015)’s estimation but $\sim$ 9-17% lower than other three curves. However, the measurements of ${V_{\odot}}$ and $R_{\odot}$ have been progressively improving. For example, if we take the GRAVITY results (Gravity Collaboration et al., 2018) used in Eilers et al. (2019) to convert the coordinates, it does bring our curve close to their measurements. As shown in the right panel of the figure, measurements using GRAVITY solar coordinates are higher than Bovy (2015)’s and comparable to those of Eilers et al. (2019) from Gaia DR2, but are still significantly lower than maser and globular cluster measurements by $\sim$ 5-10%.

Let us now comment on our choice of distance (or metric) in the action space (Equation 14). The reason why we normalize action variables by their standard deviation to compute distance is partly due to the assumption we made in the likelihood derivation in Appendix A. Our derivation starts from a uniformly distributed background plus gaussian fluctuations which model clustering in action space. Therefore, the structures we consider should be on smaller scale than the background distribution in the action space. Since the extent of the background could be different in different directions in the action space, the normalization has the effective role of making the distribution homogeneous and isotropic, at least for $D\ll 1$, i.e. close pairs.

Note that here we ignore the covariance between different action variables. As a sanity check, we modified our distance definition to $D^{2}=\sum_{i,j=1}^{3}\Delta J_{i}\Delta J_{j}\sigma_{ij}^{-2}$ (where $\sigma_{ij}^{-2}$ is the inverse covariance matrix of $J_{i}$’s over the entire sample) accounting for the correlation between action variables. We checked this using both our real data and simulations, and found no significant change in our results (e,g., for 11.5$<$R$<$15 kpc and same $\ln D_{\rm max}$, the relative changes in $f_{h}$ and $\alpha$ are less than 1%).

For the choice of the free parameter $\ln D_{\textrm{max}}$, the main criterion is that we do not expect the constraints on the parameters to significantly vary with $\ln D_{\textrm{max}}$. Therefore, when $\ln D_{\rm max}\lesssim\ln D_{\rm max,optimum}$, the constraints on both parameters should not be a strong function of $\ln D_{\textrm{max}}$ and also be self-consistent within error bars. Meanwhile, the measurements are better to be the least uncertain at the $\ln D_{\rm max,optimum}$. One of the reasons why we plot the median values (with error bars) and two-point correlation function, $\ln\left[\frac{P(\ln D)}{D^{3}}\right]$, as a function of $\ln D_{\rm max}$ (or generally, $\ln D$) is to see which $\ln D_{\rm max}$ value can give us stable and reliable constraints. Therefore, we do not expect the optimum choice of $\ln D_{\rm max}$ to be necessarily the same for different systems. This criterion is further explored in Appendix C, where we include a background in our simulations, leading to a different $\ln D_{\rm max,optimum}$.

Here, for real data analysis, we choose ln $D_{\textrm{max, optimum}}=$ -1.14 as mentioned in previous section, and we also choose the same range for the stream-only simulation for consistency (as under this specific case, neither parameters drastically change with ln $D_{\textrm{max}}$). To give more intuition about what the chosen $D_{\rm max}$ physically means, we calculate the values of two-point correlation function using the Gaia-Enceladus globular clusters data (Myeong et al., 2018b), which is shown as orange dashed line in Figure 1. By doing this comparison, we can see that maximum separations of the pairs we considered (within the considered $D_{\rm max}$) are slightly smaller, but comparable to the size of Gaia-Enceladus globular clusters distribution (also known as Gaia Sausage) in the action space, which indicate that the largest size of the structure contributing to the likelihood estimate is close to the characteristic size of the Gaia-Enceladus structure. We suggest that the choice of $\ln D_{\rm max}$ needs to be inferred from the behavior of two-point correlation function (as a point of transition from clustered streams to background) and $f_{h,\rm fit}$($\alpha_{\rm fit}$) as a function of $\ln D$ or $\ln D_{\rm max}$, which does help to determine the point where our method and the estimations of parameters are still reliable.

As the assumption made in our likelihood derivation (Appendix A) is a uniform background plus a random gaussian field, one might be skeptical about the validity of systematic error estimations with a stream-only simulation as presented in Section 5.1. To improve this, we conducted another simulation with the inclusion of a background. Background stars are directly taken from Gaia DR2 and their action variables are calculated in a simulated host potential with [$f_{h}$ = 0.35, $\alpha$ = 1.70] using the ‘Stäckel approximation’. Based on this, we subsequently randomize the stellar distribution in action space by adding a random gaussian scatter to each $J_{i}$ respectively, then stars with randomized action distribution are transformed back to ($\vec{X}$,$\vec{V}$) in cylindrical coordinate using the TorusMapper code (Binney & McMillan, 2016). We combine three original streams with the simulated background and evolve them with different choices of ($f_{h}$, $\alpha$). We would like to test whether the parameters recovered by the likelihood function is consistent with the initial input of the simulated host potential (which is [$f_{h}$ = 0.35, $\alpha$ = 1.70]). A detailed analysis is presented in Appendix C. As can be seen from the error bar plots (Figure 14), wherever the constraints do not significantly depend on $\ln D_{\rm max}$, the inclusion of background actually improves the measurements (at $\ln D_{\rm max,optimum}$, the constraint we get is $f_{h}=0.352\pm 0.003$, $\alpha=1.678\pm 0.058$ from quadratic fit, while for a stream-only simulation, we have $f_{h}=0.352\pm 0.001$, $\alpha=1.634\pm 0.014$). Therefore, we conclude that the systematic error estimated from a stream-only simulation should be conservative and can be propagated to further analysis.

As we noted in Sec. 1, there are other proposals to use the action-angle (or similar) variables to constrain the potential. Sanders & Binney (2013) use the correlations in the angle-frequency999For an integrable system, frequencies can be simply thought as another coordinate system in the action space. space for stars of a single stream to constrain gravitational potential. For a true potential, the angle and frequency differences of stars in a long narrow stream should lie along a straight line. An incorrect potential could cause a misalignment between the stream orbit and the underlying progenitor orbit. By minimizing this misalignment, which is potential-dependent, they manage to recover the expected constraints to a spherical logarithm potential using a simulated tidal stream. While this method uses more information (i.e. angle variables) than ours, and thus can be potentially more precise, it requires identifying only stream stars and relies on the assumption of a cold stream, which does limit its precision and accuracy. Magorrian (2014) relates the clumpiness of the stellar action-space distribution to the potential and define a likelihood, where the stellar action-space distribution is drawn from a Dirichlet process. This study justifies the viability of constraining the potential using stellar action-space distribution. However, instead of assuming a specific functional form of stellar distribution in the action space, our method is independent of f($\bf J$) as we marginalized over all possible distribution of f($\bf J$) in the derivation of likelihood function (or more precisely, what we assumed here is the probability functional P{f($\bf J$)} is gaussian with an arbitrary 2-point function that only depends on $D$.). Methods introduced in Peñarrubia et al. (2012), Sanderson et al. (2015), Sanderson et al. (2017), and Buckley et al. (2019) are the closest compared with our methodology, which minimize relative entropy (or KLD) of a system in the space of action variables (or more generally, integrals of motion). However, these studies do not directly connect their statistical representations to the two-point correlation function in the action space as proposed in this study. Indeed, our derivation in Appendix A suggests that relative entropy of the distance distribution $P(\ln D)$ [rather than $f({\bf J})$] is more directly related to the likelihood. On a more practical note, given that the density of stars (or pairs of stars) is discrete, the answer does depend on the coarse-graining procedure. However, since there are many more stellar pairs than stars ($N(N-1)/2$ vs $N$), our likelihood computation is much more robust to coarse-graining. Furthermore, to our knowledge, none of these methods have yet been applied to real data.

Finally, to be fair, we should also highlight some of the caveats in our study. Several assumptions are made in the derivation of our likelihood test in Appendix A, most importantly that of a uniform background with statistically uniform gaussian fluctuations in the action space. How much do gravitationally bound structures or non-uniformity of the background can bias our finding? While the latter effect is partially captured by the dependence on $D_{\rm max}$, a more systematic test using numerical simulations of galaxy formation may be more satisfactory. Another point of concern is the dependence of the best-fit parameter on the selection cuts. While, this could signal the inadequacy of our current potential model (either for stellar disk or dark matter halo, which may need more free parameters), it could also signal deeper problems such as errors in computing the action variables for Milky Way potential (using methods in Sec. 2), or their non-adiabatic evolution.

7 Conclusion

In this work, we develop a novel method to constrain gravitational potentials from small-scale clustering in the space of action variables, and use it to provide precise constraints on the Milky Way dark matter halo potential within 9-15 kpc from Gaia DR2. We first derive the likelihood function for different host potentials, with the assumption that the stellar distribution in the action space is a uniformly distributed background with correlated gaussian fluctuations on small scales, showing that it can be written as an integral over the two-point correlation function evaluated in the action space. The main advantage of our method is that, contrary to past studies, it does not require identification of any kind of compact structures or streams beforehand, assume the circular orbits, or any equilibrium state of the distribution. We first check the viability of our method in simulations of streams with different host potentials, showing that it recovers the normalization (slope) of the host potential with less than 1% (4%) systematic error (while stochastic errors shrink with the number of stars in the sample). We then apply our method to analyze two samples from Gaia DR2 over radial ranges of 9-11 kpc and 11.5-15 kpc, and studied the effect of selection cuts on the final results. Including all the known systematic errors, we find the parameters $(f_{h},\alpha)=(0.391\pm 0.009,1.835\pm 0.092)$ and $(0.351\pm 0.012,1.687\pm 0.079)$, for 9-11 kpc and 11.5-15 kpc respectively, for the median and 68% CI uncertainty from the posterior distribution. For both simulations and real data, we can visually confirm that the potential that maximizes the likelihood function does indeed correspond to the largest two-point correlation function and most compact distribution in the action space, which again, demonstrates the reliability of our method.

We would like to clarify that the fraction of DM, $f_{h}$, characterizes the fraction of halo component contribution to the radial force extrapolated to the position of the Sun. Given the uncertainty, our final index constraints obtained at two radial bins are consistent within $\sim 1.2\sigma$. Based on the NFW prediction, the absolute value of index gets larger at outer radii. Therefore, a power law extrapolation back to the solar radius would lead to a larger local halo density. Therefore, the fact that we find a larger value of $f_{h}$ in outer radii is consistent with the expectations from NFW (or any profile that gets steeper at larger radii).

To our knowledge, this is the first study that constrains the halo potential of the Milky Way using the action space clustering with real data. While more work is needed to fully understand the systematic error of this method (as discussed in Section 6), its sheer statistical power is formidable as it scales with the number of all the stars in the sample, and with proper calibrations can provide exquisite constraints on dark matter potential. Further improvements (or checks) may come from identification of streams beforehand or other criteria to separate disk and halo components (Bonaca et al., 2017; Helmi et al., 2017; Myeong et al., 2018a; Necib et al., 2018). Additionally, in this study, we only varied the parameters in the local dark matter halo density profile but kept the stellar disk potential fixed. More robust constraints, left for future work, requires varying the parameters in the disk potential as well, and possibly include other probes of stellar density.

Acknowledgements

We would like to thank Ana Bonaca, Jo Bovy, Charlie Conroy, and David Hogg for useful discussion.

This work was supported by the University of Waterloo, Natural Sciences and Engineering Research Council of Canada (NSERC), and the Perimeter Institute for Theoretical Physics. Research at the Perimeter Institute is supported by the Government of Canada through Industry Canada, and by the Province of Ontario through the Ministry of Research and Innovation.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Derivation of the Likelihood Function

In this section, we discuss the derivation of the likelihood which we use to constrain the parameters of the Milky Way potential. Using Bayes’ Theorem, we have,

[TABLE]

where $\Theta$ are the parameters we want to constrain, while ${\cal D}$ is the data we have available. In particular, in the context of constraining the Milky Way potential which is our goal here, $\Theta$ represents the parameters of the potential, $(f_{h},\alpha)$.

Let us divide the action space into $M$ small bins such that the average number of stars per bin (over all of action space) is $\bar{n}\ll 1$. However, because of hierarchical structure formation, the expected number counts of stars in different bins will not be independent of each other. In order to capture this, we assume that the star count in bin $a$ is a Poisson sampling of a mean, $\bar{n}+\bar{n}\chi_{a}$, where $\chi_{a}$’s are correlated random gaussian variables. Therefore, the probability of measuring star counts $\{\nu_{a}\}$ is given by:

[TABLE]

From this definition, it follows that

[TABLE]

is the covariance matrix of the random gaussian variables $\{\chi_{a}\}$. Since we do not directly observe $\chi_{a}$’s, we should marginalize over them. Therefore, in the limit of $M\to\infty$ (i.e. when $\nu_{a}=0$ or $1$) the posterior (Equation 21) is given by:

[TABLE]

where $a_{k}$ is the action-space bin in which the k-th star lies, and we assume $\sum_{a}\chi_{a}=0$ over the entire action space. Using Wick’s theorem, the above Gaussian integral can be expressed in terms of a sum of the product of 2-point functions over all possible pairings of stars. This yields the likelihood (defined as the ratio of posterior to prior):

[TABLE]

Now, for a large number of pairs, we can expect that the exponent in this expression have smaller and smaller relative fluctuations around its mean for different possible pairings. This is often known as the mean-field approximation in statistical mechanics (where the sum represents the partition function and the exponent is proportional to the energy), and allows us to move the average inside the exponent:

[TABLE]

This equation defines our log-likelihood formula adopted in Equation 4 in the main text (and subsequent statistical analyses), where we further assume that the 2-point function $\xi$ only depends on the normalized distance $D$ (Equation 14) in the action space.

Appendix B Sanity check for simulations

To check whether the maximum likelihood test does correspond to the most clustering in the action space for either simulations, we plot the 2D projection of stellar distribution and two-point correlation function varying the choice of potential parameters $f_{h}$ and $\alpha$ used in the action computation. Here, we summarize the results. Figure 11 shows how a 2D projection of stellar distribution in the action space varies with different choices of potential for one of the simulations. In these figures, $f_{h}$ ($\alpha$) is fixed, while $\alpha$ ($f_{h}$) is varying across its correct value. As expected, for both simulations, the most compact distributions occur when parameters approach the correct values for the simulation (middle panel in both figures). This is also verified in the behavior of the two-point correlation function in Figure 12, where one of the parameters is fixed and the other one is varying. For both simulations, we see that the two-point correlation function is indeed maximized around the expected value, which proves the viability of our method.

Appendix C Simulation with the inclusion of a background

In Section 5.1, we conducted a simulation using a system that is entirely composed of stream stars. The discrepancies between the initial simulated host potentials and the parameters recovered by simulation are regarded as the systematic errors, and the systematic errors are further propagated to the real data analysis. However, as the assumption made in the likelihood derivation (Appendix A) is that the stellar distribution in action space is a uniform background plus gaussian fluctuations, in this section, we present the results obtained from another set of simulations with a realistic background.

To generate a uniform background in the action space, we take real observations from Gaia DR2 and then calculate their actions using a power law potential with [$f_{h}$ = 0.35, $\alpha$ = 1.70]. Then, we randomize the action distribution by adding a random gaussian scatter to each $J_{i}$. The scatter is generated from a gaussian distribution $N(0,\sigma^{2}_{\rm scatter})$ , where $\sigma^{2}_{\rm scatter}=\frac{1}{3}<J_{i}^{2}>D_{\rm max}^{2}$ with $\ln D_{\rm max}$ = -1 (equivalent to coarse-graining the action distribution using a gaussian filter of width $\ln D_{\rm max}$ = -1). Then, the randomized action variables are transformed back to the position and velocity in the cylindrical coordinates using the TorusMapper code developed by Binney & McMillan (2016), which is implemented using galpy package. Following these steps, a set of background stars, which are generated by randomizing the action distribution calculated from Gaia real data, can be produced. Therefore, we can then combine this background with three stream stars that are evolved in the same host potential for further analysis. Using the likelihood function defined in Equation 4, we will then test whether the recovered potential that corresponds to the most clustering distribution in action space is the same as the initial input.

Here, we include a wider prior (Equation 5.2) to capture whether there might be any double peak features that are presented in real data analysis. Figure 13 shows the posterior distribution of two parameters. Although the posterior distribution of each parameter is dominated by a single peak around the correct value, at some values of $\ln D_{\rm max}$, the inclusion of the background stars seems to cause the distribution present a double peak features. In Figure 14, we present the correlation function as a function of $\ln D$ and the error bar plots, where the error bars are determined by the quadratic fit (as introduced in Section 5.1) and from the posterior distribution (as introduced in Section 5.2), respectively. From the top panel, we can see that the two-point correlation is maximized when approaching the correct values. To be noted that there is a “dip” presented in the correlation function, and this is also reflected in the error bar plots: both of the $f_{h}$ and $\alpha$ estimation do not drastically change with $\ln D_{\rm max}$ when $\ln D_{\rm max}\lesssim$ -2. This justifies our criteria for choosing the free parameter $\ln D_{\rm max}$. We want to choose a value of $\ln D_{\rm max}$ which gives the least uncertain measurements, and at $\ln D_{\rm max}\lesssim\ln D_{\rm max,optimum}$, the estimations of both parameters should be stabilized (not a strong function of $\ln D_{\rm max}$) and the constraints need to be all consistent with each other (within error bars). For this specific set of simulation, it should be $\ln D_{\rm max,optimum}$ $\sim$ -2 (while in real data, we chose $\ln D_{\rm max,optimum}$ $\sim$ -1). That is the point where our method is still valid and the constraints can be safely obtained with confidence. Also, compare the systematic errors obtained from this simulation with those got from Section 5.1, it seems like the inclusion of background stars improve the constraints, “bringing” the constraints recovered from likelihood function closer to the correct values (for example, at $\ln{D_{\rm max}}\sim$ -2, the constraint we get here is $f_{h}=0.352\pm 0.003$ $\alpha=1.678\pm 0.058$ from quadratic fit and $f_{h}=0.352\pm 0.003$ $\alpha=1.657^{+0.066}_{-0.051}$ from median and posterior distribution, while for a stream-only simulation, we have $f_{h}=0.352\pm 0.001$ $\alpha=1.634\pm 0.014$ at same $\ln{D_{\rm max}}$ from quadratic fit only). Therefore, we here conclude that the systematic errors estimated from a stream-only simulation should be conservative, and we thus propagate those systematic errors to the real data analysis.

Appendix D Combining measurements with unknown systemics

Here, we discuss how to combine measurements $x_{i}$ (of a single quantity $x$) that have independent known stochastic gaussian errors $\sigma_{i}$, as well as an unknown (but independent) systematic gaussian error $\sigma_{\rm sys}$. The joint likelihood is given by:

[TABLE]

Now, assuming a flat prior on $\sigma_{\rm sys}$, up to some maximum $\sigma_{\rm sys,max}$, we can find the posterior on the parameter $x$:

[TABLE]

In practice, based on the difference in results we find with and without the cuts, we make conservative choices of $\sigma_{\rm sys,max}=$ 0.02 and 0.2 for $f_{h}$ and $\alpha$ determinations, respectively.