Defining a Milky Way potential model

In [1]:

# Third-party dependencies
from astropy.io import ascii
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

# Gala
from gala.mpl_style import mpl_style
plt.style.use(mpl_style)
import gala.dynamics as gd
import gala.integrate as gi
import gala.potential as gp
from gala.units import galactic
%matplotlib inline

Introduction

gala provides a simple and easy way to access and integrate orbits in an approximate mass model for the Milky Way. The parameters of the mass model are determined by least-squares fitting the enclosed mass profile of a pre-defined potential form to recent measurements compiled from the literature. These measurements are provided with the documentation of gala and are shown below. The radius units are kpc, and mass units are solar masses:

In [2]:

tbl = ascii.read('data/MW_mass_enclosed.csv')

In [3]:

tbl

Out[3]:
Table length=16
rMencMenc_err_negMenc_err_posref
float64float64float64float64str27
0.0130000000.010000000.010000000.0Feldmeier et al. (2014)
0.12800000000.0200000000.0200000000.0Launhardt et al. (2002)
8.189502860861.54994562473.84858963492.61Bovy et al. (2012)
8.3110417867208.04475949382.74387023236.02McMillan (2011)
8.4102421035407.016733918715.615468328224.5Koposov et al. (2010)
19.0208023299175.044317988008.434833267089.9Kuepper et al. (2015)
50.0539884832748.019995734543.3268490735257.0Wilkinson & Evans (1999)
50.0529886965173.09997867269.6638536752776.2Sakamoto et al. (2003)
50.0399914690707.0109976539941.072696676468.3Smith et al. (2007)
50.0419910425326.039991469077.038172735113.6Deason et al. (2012)
60.0399914690958.069985070910.664344945146.9Xue et al. (2008)
80.0689852841359.0299936018002.0110361048549.0Gnedin et al. (2010)
100.01.39970141731e+12899808054059.0831336271726.0Watkins et al. (2010)
120.0539884832260.0199957345315.0123854764645.0Battaglia et al. (2005)
150.0750000000000.0250000000000.0250000000000.0Deason et al. (2012)
200.0679854974257.0409912558030.0313652012195.0Bhattacherjee et al. (2014)

Let’s now plot the above data and uncertainties:

In [4]:

fig, ax = plt.subplots(1, 1, figsize=(4,4))

ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']),
            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa',
            capthick=0, linestyle='none', elinewidth=1.)

ax.set_xlim(1E-3, 10**2.6)
ax.set_ylim(7E6, 10**12.25)

ax.set_xlabel('$r$ [kpc]')
ax.set_ylabel('$M(<r)$ [M$_\odot$]')

ax.set_xscale('log')
ax.set_yscale('log')

fig.tight_layout()
../_images/potential_define-milky-way-model_6_0.png

We now need to assume some form for the potential. For simplicity and within reason, we’ll use a four component potential model consisting of a Hernquist (1990) bulge and nucleus, a Miyamoto-Nagai (1975) disk, and an NFW (1997) halo. We’ll fix the parameters of the disk and bulge to be consistent with previous work (Bovy 2015 - please cite that paper if you use this potential model) and vary the scale mass and scale radius of the nucleus and halo, respectively. We’ll fit for these parameters in log-space, so we’ll first define a function that returns a gala.potential.CCompositePotential object given these four parameters:

In [5]:

def get_potential(log_M_h, log_r_s, log_M_n, log_a):
    mw_potential = gp.CCompositePotential()
    mw_potential['bulge'] = gp.HernquistPotential(m=5E9, c=1., units=galactic)
    mw_potential['disk'] = gp.MiyamotoNagaiPotential(m=6.8E10*u.Msun, a=3*u.kpc, b=280*u.pc,
                                                     units=galactic)
    mw_potential['nucl'] = gp.HernquistPotential(m=np.exp(log_M_n), c=np.exp(log_a)*u.pc,
                                                 units=galactic)
    mw_potential['halo'] = gp.NFWPotential(m=np.exp(log_M_h), r_s=np.exp(log_r_s), units=galactic)

    return mw_potential

We now need to specify an initial guess for the parameters - let’s do that (by making them up), and then plot the initial guess potential over the data:

In [6]:

# Initial guess for the parameters- units are:
#     [Msun, kpc, Msun, pc]
x0 = [np.log(6E11), np.log(20.), np.log(2E9), np.log(100.)]
init_potential = get_potential(*x0)

In [7]:

xyz = np.zeros((3, 256))
xyz[0] = np.logspace(-3, 3, 256)

fig, ax = plt.subplots(1, 1, figsize=(4,4))

ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']),
            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa',
            capthick=0, linestyle='none', elinewidth=1.)

fit_menc = init_potential.mass_enclosed(xyz*u.kpc)
ax.loglog(xyz[0], fit_menc.value, marker='', color="#3182bd",
          linewidth=2, alpha=0.7)

ax.set_xlim(1E-3, 10**2.6)
ax.set_ylim(7E6, 10**12.25)

ax.set_xlabel('$r$ [kpc]')
ax.set_ylabel('$M(<r)$ [M$_\odot$]')

ax.set_xscale('log')
ax.set_yscale('log')

fig.tight_layout()
../_images/potential_define-milky-way-model_11_0.png

It looks pretty good already! But let’s now use least-squares fitting to optimize our nucleus and halo parameters. We first need to define an error function:

In [8]:

def err_func(p, r, Menc, Menc_err):
    pot = get_potential(*p)
    xyz = np.zeros((3,len(r)))
    xyz[0] = r
    model_menc = pot.mass_enclosed(xyz).to(u.Msun).value
    return (model_menc - Menc) / Menc_err

Because the uncertainties are all approximately but not exactly symmetric, we’ll take the maximum of the upper and lower uncertainty values and assume that the uncertainties in the mass measurements are Gaussian (a bad but simple assumption):

In [9]:

err = np.max([tbl['Menc_err_pos'], tbl['Menc_err_neg']], axis=0)
p_opt, ier = leastsq(err_func, x0=x0, args=(tbl['r'], tbl['Menc'], err))
assert ier in range(1,4+1), "least-squares fit failed!"
fit_potential = get_potential(*p_opt)

Now we have a best-fit potential! Let’s plot the enclosed mass of the fit potential over the data:

In [10]:

xyz = np.zeros((3, 256))
xyz[0] = np.logspace(-3, 3, 256)

fig, ax = plt.subplots(1, 1, figsize=(4,4))

ax.errorbar(tbl['r'], tbl['Menc'], yerr=(tbl['Menc_err_neg'], tbl['Menc_err_pos']),
            marker='o', markersize=2, color='k', alpha=1., ecolor='#aaaaaa',
            capthick=0, linestyle='none', elinewidth=1.)

fit_menc = fit_potential.mass_enclosed(xyz*u.kpc)
ax.loglog(xyz[0], fit_menc.value, marker='', color="#3182bd",
          linewidth=2, alpha=0.7)

ax.set_xlim(1E-3, 10**2.6)
ax.set_ylim(7E6, 10**12.25)

ax.set_xlabel('$r$ [kpc]')
ax.set_ylabel('$M(<r)$ [M$_\odot$]')

ax.set_xscale('log')
ax.set_yscale('log')

fig.tight_layout()
../_images/potential_define-milky-way-model_17_0.png

This potential is already implemented in gala in gala.potential.special, and we can import it with:

In [11]:

from gala.potential import MilkyWayPotential

In [12]:

potential = MilkyWayPotential()
potential

Out[12]:

<CompositePotential disk,bulge,nucleus,halo>