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

SCF and Multipole Expansions

galpy supports two complementary basis-function-expansion techniques for representing the gravitational potential of arbitrary density distributions:

• SCFPotential uses the self-consistent field (SCF) method of Hernquist & Ostriker (1992), which expands a density in a biorthogonal basis set built from the Hernquist potential. This works well for smooth, roughly spheroidal systems.

• MultipoleExpansionPotential uses a spherical-harmonic (multipole) expansion where the radial functions are computed by direct integration. This is often more accurate than SCF for non-spherical systems, and it naturally supports time dependence.

For background on multipole expansions, see Chapter 12.3.1 of Bovy (2026). The SCF method is based on Hernquist & Ostriker (1992).

Both approaches let you turn any density function into a potential that can be used for orbit integration, action-angle calculations, and all other galpy operations. This tutorial demonstrates how to use both, and also covers the disk-specific variants DiskSCFPotential and DiskMultipoleExpansionPotential.

%matplotlib inline
import warnings

warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)
import numpy
from matplotlib import pyplot as plt
from galpy import potential
from galpy.potential import (
    SCFPotential,
    MultipoleExpansionPotential,
    DiskSCFPotential,
    TriaxialNFWPotential,
    SolidBodyRotationWrapperPotential,
    HernquistPotential,
    DoubleExponentialDiskPotential,
    scf_compute_coeffs_spherical,
)
from galpy.orbit import Orbit

SCFPotential from a density

The SCFPotential.from_density class method builds an SCF expansion from any density function. As an example, we use a TriaxialNFWPotential with c=1.4, which stretches the density along the \(z\)-axis and makes the isodensity surfaces oblate (flattened in \(R\)). We use a large SCF scale length a=50 and expansion orders N=80, L=40.

tri_nfw = TriaxialNFWPotential(normalize=1.0, c=1.4, a=1.0)

scf_tri = SCFPotential.from_density(
    tri_nfw.dens, 80, L=40, a=50.0, symmetry="axisymmetry"
)

Here symmetry='axisymmetry' indicates that we are assuming axisymmetry in the basis-function expansion; other valid values are symmetry='spherical' when assuming spherical symmetry or symmetry=None for the general, non-axisymmetric computation.

We can check how well the SCF expansion reproduces the true density by comparing along the \(R = z\) line:

xs = numpy.linspace(0.01, 3.0, 1001)
plt.loglog(xs, [tri_nfw.dens(x, x) for x in xs], label="True density")
plt.loglog(xs, scf_tri.dens(xs, xs), "--", label="SCF expansion")
plt.xlabel(r"$R = z$")
plt.ylabel(r"Density")
plt.legend(frameon=False)
plt.title("TriaxialNFW: SCF density comparison along R = z");

The SCF expansion matches the true density very well over several orders of magnitude.

MultipoleExpansionPotential from a density

MultipoleExpansionPotential.from_density provides an alternative expansion using spherical harmonics. We build one for the same oblate TriaxialNFW density.

Tip

Because MultipoleExpansionPotential supports time-dependent densities, it checks whether the density function accepts a t parameter. If you pass a method like tri_nfw.dens (which has a t keyword), it will be interpreted as time-dependent and require a tgrid. To avoid this, pass the potential object directly (as below) or wrap the density in a lambda that strips the t parameter: lambda R, z, phi: pot.dens(R, z, phi).

mep_tri = MultipoleExpansionPotential.from_density(
    tri_nfw, L=40, symmetry="axisymmetry"
)

Now we can compare all three densities (true, SCF, and multipole) in a single plot:

xs = numpy.linspace(0.01, 3.0, 1001)
plt.loglog(xs, [tri_nfw.dens(x, x) for x in xs], label="True density")
plt.loglog(xs, scf_tri.dens(xs, xs), "--", label="SCF expansion")
plt.loglog(xs, mep_tri.dens(xs, xs), ":", label="Multipole expansion")
plt.xlabel(r"$R = z$")
plt.ylabel(r"Density")
plt.legend(frameon=False)
plt.title("True vs. SCF vs. Multipole density along R = z");

Both expansion methods faithfully reproduce the gravitational field and can be used for orbit integration, action-angle calculations, and all other galpy operations.

For example, we can integrate an orbit in the SCFPotential approximation to the triaxial NFW potential and compare it to the orbit in the true potential:

o = Orbit([1.0, 0.1, 1.1, 0.1, 0.3, 0.0])
ts = numpy.linspace(0.0, 100.0, 10001)
o.integrate(ts, tri_nfw)
o.plot()
o.integrate(ts, scf_tri)
o.plot(overplot=True)
o.integrate(ts, mep_tri)
o.plot(overplot=True)
plt.title("Orbit in TriaxialNFW vs. SCF vs. Multipole");

Computing SCF coefficients manually

Instead of using from_density, you can compute the SCF expansion coefficients yourself and pass them directly to SCFPotential. This is useful when you want to inspect the coefficients or store them for later use.

The Hernquist profile is the lowest-order SCF basis function, so only the zeroth-order coefficient should be non-zero in the following example:

hp = HernquistPotential(amp=1.0, a=2.0)
Acos, Asin = scf_compute_coeffs_spherical(hp.dens, 10, a=2.0)
print("Acos coefficients:")
print(Acos)
print(
    "\nAs expected, only Acos[0,0,0] = 1.0 is non-zero; all others are at machine precision."
)

Acos coefficients:
[[[ 1.00000000e+00]]

 [[ 8.65722568e-18]]

 [[-1.13278379e-16]]

 [[-2.93522527e-18]]

 [[-3.04900866e-17]]

 [[-8.03257697e-20]]

 [[-1.36005913e-17]]

 [[ 1.30909657e-18]]

 [[-5.31153122e-18]]

 [[-6.69176217e-19]]]

As expected, only Acos[0,0,0] = 1.0 is non-zero; all others are at machine precision.

To build an SCFPotential from these coefficients:

sp_hernquist = SCFPotential(Acos=Acos, Asin=Asin, a=2.0)
print(sp_hernquist)

SCFPotential with internal parameters: amp=0.5, Acos=[[[ 2.82094792e-01]]

 [[ 2.44215828e-18]]

 [[-3.19552407e-17]]

 [[-8.28011762e-19]]

 [[-8.60109462e-18]]

 [[-2.26594813e-20]]

 [[-3.83665598e-18]]

 [[ 3.69289325e-19]]

 [[-1.49835529e-18]]

 [[-1.88771125e-19]]], Asin=[[[0.]]

 [[0.]]

 [[0.]]

 [[0.]]

 [[0.]]

 [[0.]]

 [[0.]]

 [[0.]]

 [[0.]]

 [[0.]]], a=2.0 and physical outputs off

Note that the internally-stored coefficients are modified by the normalization factors of the basis functions, so they are not exactly the same as the input coefficients.

DiskSCFPotential

For disk-like density distributions, the standard SCF expansion converges slowly because spherical basis functions are a poor match for thin, flat structures. DiskSCFPotential uses a trick to greatly improve convergence for disky systems.

The disk density is approximated as \(\rho_{\mathrm{disk}}(R,\phi,z) \approx \sum_i \Sigma_i(R)\,h_i(z)\), with \(h_i(z) = \mathrm{d}^2 H(z) / \mathrm{d} z^2\), and the potential is written as

\[\Phi(R,\phi,z) = \Phi_{\mathrm{ME}}(R,\phi,z) + 4\pi G\sum_i \Sigma_i(r)\,H_i(z)\]

where \(r^2 = R^2+z^2\) is the spherical radius. The density giving rise to \(\Phi_{\mathrm{ME}}\) is not strongly confined to a plane and can be obtained using the multipole or SCF methods. This trick is due to Kuijken & Dubinski (1995).

As an example, we represent a DoubleExponentialDiskPotential with \(h_R = 1/3\) and \(h_z = 1/27\):

dp = DoubleExponentialDiskPotential(amp=13.5, hr=1.0 / 3.0, hz=1.0 / 27.0)

dscfp = DiskSCFPotential(
    dens=lambda R, z: dp.dens(R, z),
    Sigma={"type": "exp", "h": 1.0 / 3.0, "amp": 1.0},
    hz={"type": "exp", "h": 1.0 / 27.0},
    a=1.0,
    N=10,
    L=10,
)

Compare the density along the \(R = 10z\) line:

xs = numpy.linspace(0.3, 2.0, 1001)
plt.semilogy(xs, dp.dens(xs, xs / 10.0), label="DoubleExponentialDisk")
plt.semilogy(xs, dscfp.dens(xs, xs / 10.0), "--", label="DiskSCFPotential")
plt.xlabel(r"$R$")
plt.ylabel(r"Density")
plt.legend(frameon=False)
plt.title("Density along R = 10z");

The DiskSCFPotential reproduces the disk density very well. Orbit integration in DiskSCFPotential is also much faster than in DoubleExponentialDiskPotential, because the latter requires expensive numerical integrals at every force evaluation whereas DiskSCFPotential evaluates analytic basis functions:

ts = numpy.linspace(0.0, 100.0, 10001)
ic_disk = [1.0, 0.1, 0.9, 0.0, 0.1, 0.0]

o_dscf = Orbit(ic_disk)
o_dscf.integrate(ts, dscfp)
o_dscf.plot()
o_dp = Orbit(ic_disk)
o_dp.integrate(ts, dp)
o_dp.plot(overplot=True)
plt.title("Orbit in DiskSCFPotential vs. DoubleExponentialDiskPotential");

The orbits agree closely, but DiskSCFPotential is dramatically faster for orbit integration. Compare

%%timeit
o.integrate(ts, dp)
# 4.53 s ± 25.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

2.19 s ± 1.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

and

%%timeit
o.integrate(ts, dscfp)
# 57.2 ms ± 99.6 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

72 ms ± 711 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

DiskMultipoleExpansionPotential

galpy also provides DiskMultipoleExpansionPotential, which applies the same Kuijken & Dubinski disk-subtraction technique but uses a MultipoleExpansionPotential (instead of SCF) to solve the residual Poisson equation. Usage is analogous to DiskSCFPotential and can be more accurate for some density profiles.

Time-dependent MultipoleExpansionPotential

The MultipoleExpansionPotential supports time-dependent densities. If you pass a density function that accepts a t keyword argument together with a tgrid array, galpy precomputes the expansion at each time step and interpolates. This is most useful for modeling potentials with complex time dependence, e.g., a time-dependent shape. In the following, we use a rotating bar as a simple example, but note that adding simple rotation is most easily accomplished using existing (or new) wrappers.

As an example, we build a Hernquist-like density with a \(\cos(2\phi)\) bar perturbation rotating at pattern speed \(\Omega_p = 1.5\):

OmegaP = 1.5


def rotating_bar_dens(R, z, phi=0.0, t=0.0):
    phi_bar = phi - OmegaP * t
    r2 = R**2 + z**2
    rho0 = 1.0 / (r2**0.5 * (1.0 + r2**0.5) ** 3)
    return rho0 * (1.0 + 0.3 * numpy.cos(2.0 * phi_bar))


mep_td = MultipoleExpansionPotential.from_density(
    rotating_bar_dens,
    L=4,
    rgrid=numpy.geomspace(0.01, 20.0, 101),
    symmetry=None,
    tgrid=numpy.linspace(0.0, 4.0 * numpy.pi / OmegaP, 51),
)

We can visualize how the potential varies with azimuth at different times, confirming that the bar pattern rotates:

phis = numpy.linspace(0.0, 2.0 * numpy.pi, 100)
for t in [0.0, 1.0]:
    vals = [mep_td(1.0, 0.0, phi=phi, t=t) for phi in phis]
    plt.plot(phis, vals, label=f"t = {t:.1f}")
plt.xlabel(r"$\phi$")
plt.ylabel(r"$\Phi$")
plt.legend()
plt.title("Time-dependent multipole expansion: rotating bar");

The phase shift between the two curves shows the bar rotating at the specified pattern speed. This time-dependent multipole expansion can be used for orbit integration in the same way as any other galpy potential. For example, we can integrate an orbit in a rotating bar-like perturbation on top of a Hernquist halo. First we build the combined potential:

hp = HernquistPotential(normalize=1.0, a=1.0)
omega = 1.3
epsilon = 0.125
tgrid = numpy.linspace(0, 60, 251)
mp_tdep = MultipoleExpansionPotential.from_density(
    lambda R, z, phi, t=0.0: (
        hp.dens(R, z, use_physical=False)
        * (1 + epsilon * numpy.cos(2 * (phi - omega * t)))
    ),
    L=4,
    rgrid=numpy.geomspace(1e-3, 30, 201),
    tgrid=tgrid,
)

Then, we integrate an orbit in this potential and compare it to the equivalent SolidBodyRotationWrapperPotential approach:

mp_static = MultipoleExpansionPotential.from_density(
    lambda R, z, phi: (
        hp.dens(R, z, use_physical=False) * (1 + epsilon * numpy.cos(2 * phi))
    ),
    L=4,
    rgrid=numpy.geomspace(1e-3, 30, 201),
)
mp_wrapped = SolidBodyRotationWrapperPotential(pot=mp_static, omega=omega)
from galpy.orbit import Orbit

ts = numpy.linspace(0, 50, 5001)
o_wrap = Orbit([1.0, 0.1, 1.1, 0.0, 0.05, 0.3])
o_tdep = Orbit([1.0, 0.1, 1.1, 0.0, 0.05, 0.3])
o_wrap.integrate(ts, mp_wrapped)
o_tdep.integrate(ts, mp_tdep)

fig, axes = plt.subplots(1, 2, figsize=(9, 4.5))
axes[0].plot(o_wrap.R(ts), o_wrap.z(ts), lw=1.0, label="Wrapper")
axes[0].plot(o_tdep.R(ts), o_tdep.z(ts), lw=1.0, label="Time-dep.")
axes[0].set_xlabel(r"$R$")
axes[0].set_ylabel(r"$z$")
axes[0].legend(fontsize=9)
axes[1].plot(
    o_wrap.R(ts) * numpy.cos(o_wrap.phi(ts)),
    o_wrap.R(ts) * numpy.sin(o_wrap.phi(ts)),
    lw=1.0,
    label="Wrapper",
)
axes[1].plot(
    o_tdep.R(ts) * numpy.cos(o_tdep.phi(ts)),
    o_tdep.R(ts) * numpy.sin(o_tdep.phi(ts)),
    lw=1.0,
    label="Time-dep.",
)
axes[1].set_xlabel(r"$x$")
axes[1].set_ylabel(r"$y$")
axes[1].legend(fontsize=9)
axes[1].set_aspect("equal")