Inferring the Gravitational Potential of the Milky Way

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 cchu@hmc.edu 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$.

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:

• If we evolve the particles in the correct potential, they will remain in the same, steady-state configuration.
• Conversely, if we evolve them in an incorrect potential, they will evolve towards a different configuration.

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$).

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 $i$th 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.

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.

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):

• General: This algorithm requires only that the particles are in steady-state, whereas some other methods require also that the potential be integrable.
• Intuitive: It has an intuitive interpretation in terms of phase-mixing.

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

• Computation time: An application to $10^9$ particles will be computationally expensive.
• Noise: The likelihood function is currently quite noisy, which may prevent accurate inference of parameters.
• Lack of error bounds: It's currently unclear how to accurately quantify the uncertainty in the inferred parameters.

• Binney & Tremaine, 2008, Galactic Dynamics: Second Edition.
• Bovy et al., 2010, ApJ, 711, 1157.
• Han et al., 2015, arXiv:1507.00769 (pre-print).
• Magorrian, 2014, MNRAS, 437, 2230.
• Perryman et al., 2001, A&A, 369, 339.

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.