Progressive Approximations in Mass Modeling

I have said I wasn’t going to attempt to teach an entire graduate course on galaxy dynamics in this forum, and I’m not. But I can give some pointers for those who want to try it for themselves. It also provides some useful context for fans of Deur’s approach.

The go-to textbook for this topic is Galactic Dynamics by Binney & Tremaine. The first edition was published in 1987, conveniently when I switched to grad school in astronomy. It was already a deep and well-developed field at that time; this is a compendium of considerable scientific knowledge.

Fun story: a colleague in a joint physics & astronomy department once complained to me that she wanted to develop a course in galaxy dynamics, which is a staple of graduate programs in astronomy & astrophysics. However, there was a certain senior colleague who objected, saying that since it was astronomy, it couldn’t possibly be a rigorous course worthy of a full semester graduate course. This is a casual bias that astronomers often encounter when talking to physicists, many of whom have attitudes about the subject that were trapped in amber sometime in the Jurassic. I suggested that she walk into his office and drop a copy of Galactic Dynamics on his desk from on high, as (1) it would make a hefty impact, and (2) no one who so much as skims this book could persist in this toxic attitude.

She later reported that she had done this, and it had worked.

Galactic Dynamics is not a starter book. It is the textbook we use when teaching the graduate course that this is not. A useful how-to guide for the specific material I’ll discuss here is provided by Federico Lelli. In brief, to model the gravitational potential of an observed distribution of matter, we can make one of the following series of approximations:

All science is an approximation at some level. The most crude approximation we can employ here is to imagine that all of the mass resides at a central point. In this limit, the potential is simply

V² = GM/R

where V is the orbital speed of a test particle on a circular orbit, G is Newton’s constant, M is the mass, and R is the distance from the point mass. Galaxies are not point masses, so this is a terrible approximation, as can be seen by the divergent V ~ R^-1/2 behavior as R → 0 (the dotted line above).

The next bad approximation one can make is a spherical cow: assume the mass is distributed in a sphere that is projected as the image we see on the sky. This at least incorporates the fact that the mass is not all concentrated at a point, so

V² = GM(R)/R

acknowledges that the mass M is spread out as a function of radius. This is a spherical cow. Since we cannot see dark matter, we almost always assume it to be a spherical cow.

For the luminous disk of a spiral galaxy, a common approximation is the so-called exponential disk:

Σ(R) = Σ₀ e^-R/R_d

where Σ₀ is the central surface density of stars and R_d is the scale length of the disk – the characteristic size over which the surface brightness declines exponentially. This can be integrated by parts to obtain an expression for the enclosed mass M(R) which I leave as an exercise for the eager reader. This provides a handy analytic formula, the rotation curve of which is illustrated above by the dashed line.

Spiral galaxies are fairly thin when seen edge-on, so the spherical cow is not a great approximation. In a classic paper, Freeman (1970) solved the Poisson equation for the case of a razor-thin exponential disk, where one meets modified Bessel functions of the first and second kind (denoted “ikik” above). These must be solved numerically, but one can make a tabulation for use with any choice of disk mass and scale length. Such a thin disk is illustrated by the grey line above for a choice of stellar mass and scale length appropriate to NGC 6946.

Spiral galaxies are not razor thin of course. We only see a projected image on the sky, so for a galaxy like NGC 6946, we may have a good measurement of its azimuthally averaged light (and presumable stellar mass) distribution Σ(R) but we have no idea how thick it is. Here, we have to make an educated guess based on observations of edge-on galaxies. A ballpark average is R:z = 8:1, but some galaxies are thicker and others thinner, so this becomes an approximation with an associated uncertainty. This uncertainty cannot be unambiguously eliminated; it is one of the known unknowns that comprise the inevitable systematic errors in astronomy. Fortunately, allowing for a finite thickness only takes the harsh edge off of the thin disk case, and the assumption one chooses makes little difference to the result (compare the lines labeled thick and thin above).

The exponential disk formula Σ(R) is an azimuthal average over an image like that of NGC 6946. This approximation captures none of the spiral structure: it only tells us about the average rate at which the surface brightness falls off. It also imposes a smooth shape on that fall off that our eyes can see is not necessarily a great approximation. So the next level of approximation is to solve the Poisson equation numerically for the observed surface brightness profile, Σ(R), not just the exponential approximation thereto. This is the blue line in the bottom right graph above.

There are important differences between using the numerical solution for the observed light distribution and the exponential disk approximation. This has been known since the 1980s, but the analytic expression is so convenient that people need an occasional reminder not to trust it too much. Jerry Sellwood felt the need to provide this reminder in 1999:

Small apparent differences in the shape of the mass profile (left) correspond to pronounced differences in the rotation curve (right). I chose the example of NGC 6946 in part because the exponential approximation for it is pretty good. Nevertheless, the details matter, so the best practice is to build numerical mass models, as we did for SPARC.

Building numerical mass models is tractable for external galaxies, where we can see the entire light distribution. It is not possible for our own Milky Way, since we are located within it and cannot see it as a whole. Consequently, the vast majority of Milky Way models rely on the exponential approximation; so far as I’m aware, I’m the only one who has built a model that attempts to get beyond this.

Numerical mass models are still an approximation. We’re assuming that the gravitational potential is static and azimuthally symmetric. Taking the next step would require abandoning these assumptions to model the spiral arms. The Poisson equation can handle that, but it becomes dicey because the arms rotate with some pattern speed (generally unknown) and may grow or dissolve or reform on some unknown timescale. The potential at any given point is time variable even in equilibrium, so we need not just a numerical solution but a live numerical simulation to keep track of it. That can be done, but it has to be done on a case by case basis, and the answer will depend somewhat on additional assumptions that have to be introduced to run the simulation, like specifying a dark matter halo.

One can generalize further to consider the full 3D potential, e.g., to allow for asymmetry in the z-direction as well as in azimuth. One can further imagine non-equilibrium processes, such as an external perturbations. There is good evidence that the Milky Way suffers both of these effects, the passage of the Large Magellanic Cloud being one obvious and apparently large perturbation. So we are in the awkward position that the Gaia data now oblige us to consider the entire run of possible effects through non-equilibrium processes in a mass distribution that is not completely symmetric in any of the three spatial dimensions, but for the main mass component we are stuck with the inadequate approximation of an exponential disk.

Geometry appears to play a crucial role in the approach of Deur to the acceleration discrepancy problem. The essential claim is that the discrepancy correlates with flattening, with highly flattened systems like spirals evincing the classic discrepancy while spherical systems like E0 galaxies showing none. Big if true!

A useful plot appears on slide 44:

This is the one example shown that goes into the plot of many determinations of the slope a on the following slide. It being the only one, it is the only thing I have to evaluate without chasing down every other case. Looking at this, I am not inclined to do so.

At first it looks persuasive: the best fit slope is clear. There is no reason why the discrepancy should depend on the projected ellipticity of a triaxial 3D blob of stars, so this must be telling us something important. I’d be on board with that if it were true, but I’ve seen too many non-correlations masquerading as correlations to believe this one. The fitted slope is strongly influenced by the one point at large ellipticity; absent that, a slope of zero works fine. Mostly what I see here is a lot of scatter, which is normal in extragalactic astronomy. Since there are only a few points at high and low ellipticity, we don’t know what would happen if we went out and got more data. But I bet that what would happen is that the high ellipticity points would wind up looking like those in the middle: a big blob of scatter, with no significant correlation.

I’d kinda like to be wrong about this one, so I won’t even get into the theory side, which I find sorta compelling but ultimately unpersuasive. Why are gravitons confined to a disk? What happens way far out? Surely the flatness of the disk at tens of kpc is not dictating the flatness at 1000 kpc.

Surely.