Inferring the Gravitational Potential of the Milky Way

Casey Chu1,2, Yoram Lithwick1, Fabio Antonini1
1Northwestern University, 2Harvey Mudd College

Hello! I'm Casey Chu, a student a Harvey Mudd College. I was an REU student at CIERA at Northwestern University during the summer of 2015. Below is a summary of what I did, and the C++ code I wrote is available as a GitHub repository. Feel free to email me at to chat about this research or this research experience!

How do we study something we can't see?

One way to study dark matter is to observe its effect on visible objects: stars. In particular, the distribution of dark matter in our galaxy generates a gravitational potential that affects the motion of the stars we observe.

The spacecraft Gaia is currently recording the position and velocity of one billion stars in the Milky Way (Perryman et al. 2001), so we are in perfect position to extract this information that's statistically encoded in the motion of stars.

We've created a simple, numerical algorithm to extract the gravitational potential from Gaia-like data. Given a set of observed positions and velocities of particles (stars), we are able to infer the potential in which they are traveling.

The fundamental assumption: phase-mixing

The algorithm we developed relies on the system being phase-mixed. Figure 1 depicts the phase-mixing of 200 particles that are initially close together in phase space, evolving in a potential of \Phi(x) = \frac{1}{2}|x|^\alpha for \alpha = 1.5.

The phase-mixing of 200 particles. The arrows
trace the phase-space flow generated by the potential.

Notice that as time progresses, the particles approach a configuration that is macroscopically steady-state, or phase-mixed.

One important property of phase-mixing is that if we evolve the same set of particles in different potentials (e.g., different values of \alpha), then the particles tend towards different steady-state configurations.

The insight behind the inference algorithm

Let's assume that we know that the particles we observe are in a steady-state configuration. Then:

We demonstrate this in Figure 2, where a set of particles, drawn from an unknown potential, is evolved in three different trial potentials (different values of \alpha).

The evolution of observed particles in different potentials.

In this simple example, since the configuration of the particles evolved in the \alpha = 1.5 potential most resembles the original, observed distribution, we may infer that the particles most likely come from an \alpha = 1.5 potential.

The likelihood function

Instead of choosing the best potential by eye, we can construct a function \ell that quantifies how well a trial potential preserves the original configuration.

First, we define a function that encodes the distribution of our observed particles after having been evolved for a time t under a potential \Phi. We use a technique called kernel density estimation to approximate the phase-space density at time t; define a function f to be f(x,v;t \,|\, \Phi) = \frac{1}{n}\sum_i K(x_i(t) - x)\, K(v_i(t) - v), where x_i(t) and v_i(t) denote the position and velocity of the ith particle evolved to time t under the potential \Phi, and K(\cdot) is a kernel, e.g. a Gaussian centered at 0. The result is that this function f(x,v;t\,|\,\Phi) gives an approximation to the phase-space density of particles at x and v and time t.

Next, we define a time-averaged version of f: f(x,v \,|\, \Phi) = \frac{1}{T}\int_0^T f(x,v;t\,|\,\Phi) \,dt. This approximates the steady-state phase-space density of particles evolved under a potential \Phi.

Wwe can interpret this density function as a probability distribution function. Then the probability of observing the positions and velocities we observed assuming a particular potential is \Pr(x_1,v_1,\ldots,x_n,v_n\,|\,\Phi) = \prod_i f(x_i, v_i \,|\, \Phi), assuming that the particles were drawn independently.

Finally, we define the likelihood \ell to be a function of the potential instead of the parameters: \ell(\Phi \,|\, x_1,v_1,\ldots,x_n,v_n) = \Pr(x_1,v_1,\ldots,x_n,v_n\,|\,\Phi) The satisfies our original goal: this function \ell is greatest for potentials that preserve the original configuration.

Summary of the algorithm

Our algorithm to infer the potential from a set of observed positions and velocities is as follows.

  1. Guess a trial potential by choosing a set of values for the parameters of a potential, e.g., via a Monte Carlo method.
  2. Calculate the likelihood of observing the data in the chosen potential, according to the formula in the previous section.
  3. Repeat steps 1–2 until convergence of the likelihood. The most likely potential has the parameters with the greatest likelihood.

A two-dimensional test case

We have tested the algorithm for the two-dimensional logarithmic potential given by \Phi(x,y) = \log(x^2 + \frac{y^2}{q^2} + R_c^2), where q and R_c are parameters that describe the shape of the potential. It is a simple model for real galaxies.

"Stars" in a logarithmic potential.

We simulated n = 10^4 particles in the logarithmic potential with true parameters of q^2=0.8 and R_c^2=1.5. Figure 3 depicts the likelihood we calculated. It peaks near the true values of the parameters, indicating that the algorithm was successful.

A contour plot of the likelihood of the parameters. Red is higher; contours are at 2^{-2^i} for i=0,\ldots,13.

Advantages and future work

Our method has two novel advantages over existing methods (e.g., Bovy et al. 2010, Magorrian 2014, Han et al. 2015):

Nevertheless, there are several challenges that should be addressed before applying the algorithm to Gaia data:

This material is based upon work supported by the National Science Foundation under Grant No. AST-1359462, a Research Experiences for Undergraduates (REU) grant awarded to CIERA at Northwestern University. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.