Self-consistent Modelling of the Milky Way using Gaia data

TL;DR

This paper discusses using angle/action distribution function models, optimized with Gaia and other survey data, to accurately map the Milky Way's matter distribution.

Contribution

It introduces a self-consistent modeling approach that leverages Gaia data to improve the accuracy of Milky Way mass distribution estimates.

Findings

Enhanced accuracy in Milky Way modeling with Gaia data.

Demonstrated the importance of full selection function in modeling.

Improved understanding of the Galaxy's matter distribution.

Abstract

Angle/action based distribution function (DF) models can be optimised based on how well they reproduce observations thus revealing the current matter distribution in the Milky Way. Gaia data combined with data from other surveys, e.g. the RAVE/TGAS sample, and its full selection function will greatly improve their accuracy.

Equations

$$\rho (x)=\int d^{3}vf(\text{J}(\text{x},\text{v}))$$

$$f(J_r,J_z,J_{\varphi })=f_{{\sigma }_r}(J_r,J_{\varphi })f_{{\sigma }_z}(J_z,J_{\varphi }),$$

$$f_{{\sigma }_r}(J_r,J_{\varphi })\equiv \frac{\Omega \Sigma }{\pi {\sigma }_r^2\kappa }[1+\tanh (J_{\varphi }/L_0)]{\mathrm{e}}^{-\kappa J_r/{\sigma }_r^2},f_{{\sigma }_z}(J_z,J_{\varphi })\equiv \frac{\nu }{2\pi {\sigma }_z^2}{\mathrm{e}}^{-\nu J_z/{\sigma }_z^2}.$$

$$f_{\mathrm{P}}(\text{J})=\frac{N}{J_0^3}\frac{(1+J_0/h{)}^{5/3}}{(1+h/J_0{)}^{2.9}}$$

$$g(h)={\left[\frac{h_0^2}{h^2}-\beta \frac{h_0}{h}+1\right]}^{-5/6},$$

Full text

††volume: 330††journal: Astrometry and Astrophysics in the Gaia sky††editors: eds. Recio-Blanco, de Laverny & Brown.

1 Introduction

Our knowledge of how baryons were accreted by galaxies such as the Milky Way is limited however but unprecedented amounts of data are becoming available from large scale surveys. It is vital that we use these data to improve our understanding of the local group and an excellent starting point is to find the current distribution of matter in the Galaxy. We cannot observe directly the distribution of dark matter so we must use the observable components as tracers to build dynamical models from which we can discover the distribution of dark matter. If we assume the Galaxy is in statistical equilibrium, we can exploit Jeans theorem ([Jeans(1916), Jeans 1916]) and presume that the distribution function (DF) $f(x,v)$ is a function of integrals of motion $I(x,v)$ only. There is an infinite choice of integrals because any function of the integrals is also an integral but the best choices are the actions $J_{i}$, which can be uniquely complemented by canonically conjugate variables to make a complete set ($\theta$,J) for phase-space coordinates.

2 Modelling Process

The Milky Way can be built up from its components; the dark halo, thin and thick discs, stellar halo, bulge and gas disc. Each of these components is modelled either by a distribution function, $f(\textbf{J})$ (first three components above), or a fixed potential. By making a sensible estimate for the parameters of these components we can make an estimate for the total potential. We then use the Stäckel Fudge ([Binney (2012), Binney (2014), Binney 2012,2014]) to find the actions and by integrating over velocity:

[TABLE]

we estimate the density. We then solve Poisson’s equation to find the new potential and iterate this procedure until after a few iterations the total potential converges.

2.1 Disc DF

The DF of the discs in our models is superposition of the “quasi-isothermal” components introduced by [Binney & McMillan (2011), Binney & McMillan (2011)]. It has the form

[TABLE]

where $f_{\sigma_{r}}$ and $f_{\sigma_{z}}$ are

[TABLE]

The thick disc is represented by a single quasi-isothermal DF, while the thin disc’s DF is built up with a quasi-isothermal for each coeval cohort of stars. The velocity-dispersion parameters depend on $J_{\phi}$ and the age of the cohort. The star-formation rate in the thin disc decreases exponentially with time, with characteristic time scale $t_{0}=8\,{\rm Gyr}$. In addition a parameter $F_{\rm thk}$ controls the fraction of mass contributed by the thick disc.

2.2 Dark halo DF

Our dark halo DF is based on the form introduced by [Posti et al.(2015)Posti, Binney, Nipoti, & Ciotti, Posti et al. (2015)] which in isolation self-consistently generates a density distribution which has a NFW ([Navarro, Frenk & White(1997)Navarro, Frenk, & White, Navarro, Frenk & White 1997]) profile (equation 4). This DF is a function of $h(\textbf{J})$ which is an almost linear function (a homogeneous function of degree unity) of the actions $J_{i}$. The haloes generated are isotropic centrally and mildly radial when $r>r_{s}$. The anisotropy can be changed by varying the linear function of the $J_{i}$. These haloes closely resemble the haloes formed in dark-matter-only simulations. Specifically

[TABLE]

The scale action $J_{0}$ encodes the scale radius around which the slope of the radial density profile shifts from $-1$ at small radii to $-3$ far out. From equation (4) it follows that $f_{P}(\mathrm{\textbf{J}})\sim|\mathrm{\textbf{J}}|^{-5/3}$ as $|\mathrm{\textbf{J}}|\rightarrow 0$.

2.3 Observational constraints

We use several sets of observations to constrain the parameters of the DFs at various stages of the modelling. The observations include the astrometry of H2O and SiO masing stars ([Reid et al.(2014)Reid, Menten, Brunthaler, Zheng, Dame, Xu, Wu, Zhang, Sanna, Sato, Hachisuka, Choi, Immer, Moscadelli, Rygl, & Bartkiewicz, Reid et al. 2014]), the distribution of radio-frequency lines of Hi and CO emission in the longitude-velocity plane, the stellar parameters and distance estimates in the fourth RAVE data release ([Kordopatis et al.(2013)Kordopatis, Gilmore, Steinmetz, Boeche, Seabroke, Siebert, Zwitter, Binney, de Laverny, Recio-Blanco, Williams, Piffl, Enke, Roeser, Bijaoui, Wyse, Freeman, Munari, Carrillo, Anguiano, Burton, Campbell, Cass, Fiegert, Hartley, Parker, Reid, Ritter, Russell, Stupar, Watson, Bienaymé, Bland-Hawthorn, Gerhard, Gibson, Grebel, Helmi, Navarro, Conrad, Famaey, Faure, Just, Kos, Matijevič, McMillan, Minchev, Scholz, Sharma, Siviero, de Boer, & Žerjal, Kordopatis et al. 2013]) (see Section 4) and finally the vertical density profile from SDSS. We assume that the population from which the RAVE sample is drawn is identical to that studied by [Juric2008, Jurić et al. (2008)].

3 Results

[Piffl et al.(2014)Piffl, Binney, McMillan, Steinmetz, Helmi, Wyse, Bienaymé, Bland-Hawthorn, Freeman, Gibson, Gilmore, Grebel, Kordopatis, Navarro, Parker, Reid, Seabroke, Siebert, Watson, & Zwitter, Piffl et al. (2014)] used RAVE data and SDSS Juric (2008) data to constrain the mass of DM within solar radius, R0. The dark halo was included as a potential not a DF. [Binney & Piffl(2015), Binney & Piffl (2015)] used the DF in equation 4 for the dark halo in their self-consistent model of the Milky Way. In this model the halo was in self-consistent equilibrium with the other components so it had been adiabatically compressed by the baryons from its original NFW form. Hence it mimicked a scenario where baryons accumulated quiescently in the Galaxy’s dark halo. An NFW halo becomes more centrally concentrated, and with so much dark matter at low radii (Fig. 1) the matching disc contains too few stars to satisfy the microlensing data ([Popowski et al.(2005)Popowski, Griest, Thomas, Cook, Bennett, Becker, Alves, Minniti, Drake, Alcock, Allsman, Axelrod, Freeman, Geha, Lehner, Marshall, Nelson, Peterson, Quinn, Stubbs, Sutherland, Vandehei, Welch, & MACHO Collaboration, Popowski et al. 2005]). This implies that the infinite phase-space density of particles at $J=0$ characteristic of an NFW DF does not survive the accretion of baryons. The baryons cannot have accumulated entirely adiabatically but the most tightly bound dark matter particles were upscattered. This scattering of DM particles by baryons reduces their phase-space density and had greatest impact near $J=0$.

In order to model this scattering process [Cole & Binney(2017), Cole & Binney (2017)] modified the NFW DF by setting $f(\textbf{J})=g(h)f_{NFW}(h)$ with $g\sim{h}^{5/3}$ for small $h$ and $g\sim 1$ at large $h$. This shifts particles from very low $J$ to higher $J$. The functional form for $g$ is

[TABLE]

where $h_{0}$ is an arbitrary constant with the dimensions of action that sets the scale of the almost constant-density core of the DF and $\beta$ is determined by requiring the total mass of dark matter is conserved. The free parameters in $f_{\rm DM}$ are $J_{0}$, which sets the NFW scale radius, $h_{0}$, which sets the size of the dark halo’s core, and the normalisation $N$.

[Cole & Binney(2017), Cole & Binney (2017)] found (see Fig. 1) that with this modified dark halo DF their models of the Milky Way have a similarly good fit to the observations but the central regions are now dominated by baryons. The dark matter fraction is reduced and consistent with the results of surveys of microlensing events. In summary they found the local DM density $\rho_{DM}{\gtrsim}0.012$ $M_{\odot}$ pc*-3*, stellar disc scale radius $R_{d}\sim 2.9$ kpc and stellar disc mass M${}_{d}\gtrsim 4\times 10^{10}$ $M_{\odot}$.

4 Using RAVE/TGAS

Our current modelling uses RAVE data in 8 spatial bins at $R_{0}\pm$1 kpc and at 0, 0.3, 0.6, 1.0, 1.5 kpc in $\left|z\right|$ to compute velocity distributions predicted by the DF at the mean positions. In order to take advantage of Gaia DR1 ([Gaia2016, Gaia Collaboration 2016]) we are developing a method using the RAVE/TGAS observations ([Kunder2017, Kunder et al. 2017]). In order to do this we will develop a full selection function $S(s)$ for this sample. Deriving an a priori calculation of $S(s)$ needs a full chemodynamical model of the MW disc $S(s,τ,[Fe/H])$ such as used in [Schönrich(2014), Schönrich & Bergemann 2014] and in addition we need to model the exact distribution of stars in age and metallicity in the solar neighbourhood. [Schönrich2017, Schönrich and Aumer 2017] did this using the RAVE selection function of [Wojno, J., Kordopatis, G., Piffl, T., Binney, J., Steinmetz, M., Matijevič, G., Bland-Hawthorn, J., Sharma, S., McMillan, P., Watson, F., Reid, W., Kunder, A., Enke, H., Grebel, E. K., Seabroke, G., Wyse, R. F. G., Zwitter, T., Bienaymé, O., Freeman, K. C., Gibson, B. K., Gilmore, G., Helmi, A., Munari, U., Navarro, J. F., Parker, Q. A., Siebert, A., Wojno et al. (2017)] and found that at fixed metallicity $S(s,\tau)$ falls off exponentially with scale 0.12 kpc at $s>0.2$ kpc.

The selection function for TGAS is biased towards younger stars which are more likely to be seen so the kinematics will appear cooler than they really are. Our models need to take this into account. They already have age but not metallicity and so we can add metallicity by use of a suitable metallicity DF. Then we can compute the likelihoods of our model based on the resulting selection function.

5 Action based modelling software library

Our modelling is currently being rewritten using AGAMA (Action-based Galaxy Modelling Architecture) which is a library of low-level programs containing interfaces and generic routines required to create the functions described here. The main sets of functions include gravitational potential and density interfaces, action/angle interface, interface for creating gravitationally self-consistent multicomponent galaxy models etc. The code can be downloaded from https://github.com/GalacticDynamics-Oxford/Agama