This page was generated from a Jupyter notebook. You can download it here.

The Potential of N-body Simulations

galpy can set up and work with the frozen potential of an N-body simulation. This allows you to study the properties of such potentials, investigate orbits, and calculate action-angle coordinates using the galpy framework. Currently, this functionality uses axisymmetrized versions of the snapshots. The use of this functionality requires pynbody.

%matplotlib inline
import warnings

warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)
import numpy
from matplotlib import pyplot as plt
import pynbody
from galpy.potential import (
    SnapshotRZPotential,
    InterpSnapshotRZPotential,
    KeplerPotential,
)

Simple example: SnapshotRZPotential

As a first, simple example we look at the potential of a single simulation particle, which should correspond to galpy’s KeplerPotential. We can create such a single-particle snapshot using pynbody and then get its potential in galpy:

s = pynbody.new(star=1)
s["mass"] = 1.0
s["eps"] = 0.0
sp = SnapshotRZPotential(s, num_threads=1)

With these definitions, this snapshot potential should be the same as KeplerPotential with an amplitude of one, which we can verify:

from galpy.potential import KeplerPotential

kp = KeplerPotential(amp=1.0)
print("Potentials: ", sp(1.1, 0.0), kp(1.1, 0.0), sp(1.1, 0.0) - kp(1.1, 0.0))
print(
    "R forces: ",
    sp.Rforce(1.1, 0.0),
    kp.Rforce(1.1, 0.0),
    sp.Rforce(1.1, 0.0) - kp.Rforce(1.1, 0.0),
)

Potentials:  -0.9090909090909091 -0.9090909090909091 0.0
R forces:  -0.8264462809917354 -0.8264462809917353 -1.1102230246251565e-16

SnapshotRZPotential instances can be used wherever other galpy potentials can be used (note that the second derivatives have not been implemented, such that functions depending on those will not work). For example, we can plot the rotation curve:

sp.plotRotcurve();

InterpSnapshotRZPotential

Because evaluating the potential and forces of a snapshot is computationally expensive, most useful applications of frozen N-body potentials employ interpolated versions of the snapshot potential. These can be set up in galpy using InterpSnapshotRZPotential.

To illustrate its use we will make use of one of pynbody’s example snapshots, g15784. This snapshot is used in the pynbody quickstart tutorial and can be downloaded from Zenodo.

Once you have downloaded the testdata, load and prepare the snapshot:

s = pynbody.load("gasoline_ahf/g15784.lr.01024.gz")

# Get the main galaxy, center, and align face-on
h = s.halos()
h0 = h[0]
pynbody.analysis.faceon(h0)

# Convert to physical units with G=1
s.physical_units(mass="kpc km**2 s**-2 G**-1")

pynbody.halo : Unable to load AHF substructure file; continuing without. To expose the underlying problem as an exception, pass ignore_missing_substructure=False to the AHFCatalogue constructor

We can now load an interpolated version of this snapshot’s potential into galpy:

spi = InterpSnapshotRZPotential(
    h0,
    rgrid=(numpy.log(0.01), numpy.log(20.0), 101),
    logR=True,
    zgrid=(0.0, 10.0, 101),
    interpPot=True,
    zsym=True,
)

where we further assume that the potential is symmetric around the mid-plane (\(z=0\)). Depending on the size of the simulation, this instantiation can take a long time. This potential instance has physical units (and thus the rgrid= and zgrid= inputs are given in kpc if the simulation’s distance unit is kpc). For example, if we ask for the rotation curve:

spi.plotRotcurve(
    Rrange=[0.01, 19.9],
    xlabel=r"$R\,(\mathrm{kpc})$",
    ylabel=r"$v_c(R)\,(\mathrm{km\,s}^{-1})$",
);

Converting to natural units

Because galpy works best in a system of natural units, we can convert the interpolated snapshot potential to natural units using normalize(). For example, using the circular velocity at \(R = 10\) kpc:

print("Vcirc at R=10 kpc:", spi.vcirc(10.0))
spi.normalize(R0=10.0)

Vcirc at R=10 kpc: 294.6732646339408

We can then plot the rotation curve again, keeping in mind that the distance unit is now \(R_0\):

spi.plotRotcurve(Rrange=[0.01, 1.99])

[<matplotlib.lines.Line2D at 0x7f2a22319bd0>]

In particular, spi.vcirc(1.) now returns 1.0. We can also plot the potential:

spi.plot(rmin=0.01, rmax=1.9, nrs=51, zmin=-0.99, zmax=0.99, nzs=51)
plt.gca().set_aspect(1.0);

This simulation’s potential is quite spherical, which is confirmed by looking at the flattening:

print(f"Flattening at R=1.0, z=0.1: {spi.flattening(1.0, 0.1)}")
print(f"Flattening at R=1.5, z=0.1: {spi.flattening(1.5, 0.1)}")

Flattening at R=1.0, z=0.1: 0.8631229510097888
Flattening at R=1.5, z=0.1: 0.9430920419468051

The epicycle and vertical frequencies can also be interpolated by setting the interpepifreq=True or interpverticalfreq=True keywords when instantiating the InterpSnapshotRZPotential object.

Estimating action-angle coordinates in the N-body simulation

Once we have set up the interpolated potential of an N-body simulation, we can use galpy’s action-angle machinery to compute actions for simulation particles. This is a powerful way to characterize the orbital structure of the simulated galaxy. Here we use actionAngleStaeckel, the most accurate general method for axisymmetric potentials. We need to re-instantiate the interpolated potential with enable_c=True so the action-angle calculation can use the fast C implementation:

spi_c = InterpSnapshotRZPotential(
    h0,
    rgrid=(numpy.log(0.01), numpy.log(20.0), 101),
    logR=True,
    zgrid=(0.0, 10.0, 101),
    interpPot=True,
    zsym=True,
    enable_c=True,
)
spi_c.normalize(R0=10.0)

Now we load a pristine copy of the simulation (because the normalization above can lead to some inconsistent behavior in pynbody) and select star particles in a radial ring near \(R = 8\,\mathrm{kpc}\):

sc = pynbody.load("gasoline_ahf/g15784.lr.01024.gz")
hc = sc.halos()
hc0 = hc[0]
pynbody.analysis.faceon(hc0)
sc.physical_units()

# Select stars in a radial ring near R = 8 kpc
sn = pynbody.filt.BandPass("rxy", "7 kpc", "9 kpc")
R, vR, vT, z, vz = [
    numpy.ascontiguousarray(hc0.s[sn][x]) for x in ("rxy", "vr", "vt", "z", "vz")
]

# Convert to natural units
ro, vo = 10.0, spi_c.vcirc(1.0, use_physical=False) * 294.67
R /= ro
z /= ro
vR /= vo
vT /= vo
vz /= vo
print(f"Selected {len(R)} star particles")

pynbody.halo : Unable to load AHF substructure file; continuing without. To expose the underlying problem as an exception, pass ignore_missing_substructure=False to the AHFCatalogue constructor

Selected 12926 star particles

We can integrate a random orbit from this selection to see what it looks like:

from galpy.orbit import Orbit

numpy.random.seed(1)
ii = numpy.random.permutation(len(R))[0]
o = Orbit([R[ii], vR[ii], vT[ii], z[ii], vz[ii]])
ts = numpy.linspace(0.0, 100.0, 1001)
o.integrate(ts, spi_c)
o.plot();

Now compute actions for all selected particles using actionAngleStaeckel:

from galpy.actionAngle import actionAngleStaeckel

aAS = actionAngleStaeckel(pot=spi_c, delta=0.45, c=True)
jr, lz, jz = aAS(R, vR, vT, z, vz)

Plot the distribution of angular momentum \(L_z\) vs. radial action \(J_R\):

from galpy.util import plot as galpy_plot

galpy_plot.scatterplot(
    lz,
    jr,
    "k.",
    xlabel=r"$L_z$",
    ylabel=r"$J_R$",
    xrange=[0.0, 1.3],
    yrange=[0.0, 0.6],
);

The curved shape is a selection effect: low angular-momentum stars can only enter the selected radial ring if they have highly eccentric orbits (and therefore large \(J_R\)). The density gradient in angular momentum reflects the falling surface density of the disk. We can also look at the distribution of radial and vertical actions:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.plot(jr, jz, "k,")
ax.set_xlabel(r"$J_R$")
ax.set_ylabel(r"$J_z$")
ax.set_xlim(0.0, 0.4)
ax.set_ylim(0.0, 0.2);

For more on action-angle coordinates and the different methods available for computing them, see the action-angle tutorials.