Multipole Expansion Potential

class galpy.potential.MultipoleExpansionPotential(amp=1.0, rho_cos_splines=None, rho_sin_splines=None, rgrid=array([1.00000000e-03, 1.01036227e-03, 1.02083192e-03, ..., 2.93877957e+01, 2.96923201e+01, 3.00000000e+01], shape=(1001,)), tgrid=None, normalize=False, ro=None, vo=None)

Class that implements a gravitational potential computed via multipole expansion of an arbitrary density distribution.

This class decomposes a user-supplied density function into real spherical harmonics on a radial grid, then evaluates the gravitational potential using classical multipole integrals (the real-valued form of Bovy 2026, Chapter 12.3.1, equations 12.78–12.79):

\[\rho(r,\theta,\phi) = \sum_{l=0}^{L-1}\sum_{m=0}^{l}\,\left[\rho^{\cos}_{lm}(r)\,\cos(m\phi) + \rho^{\sin}_{lm}(r)\,\sin(m\phi)\right]\,P_l^m(\cos\theta)\,,\]

where the radial density coefficients are obtained by projection:

\[ \begin{align}\begin{aligned}\rho^{\cos}_{lm}(r) = \alpha_{lm} \int_0^{\pi}\!\int_0^{2\pi} \rho(r,\theta,\phi)\,P_l^m(\cos\theta)\,\cos(m\phi)\,\sin\theta\,\mathrm{d}\phi\,\mathrm{d}\theta\,,\\\rho^{\sin}_{lm}(r) = \alpha_{lm} \int_0^{\pi}\!\int_0^{2\pi} \rho(r,\theta,\phi)\,P_l^m(\cos\theta)\,\sin(m\phi)\,\sin\theta\,\mathrm{d}\phi\,\mathrm{d}\theta\,,\end{aligned}\end{align} \]

with \(\alpha_{lm} = \sqrt{\frac{2l+1}{4\pi}\,\frac{(l-m)!}{(l+m)!}}\). The gravitational potential has an analogous expansion:

\[\Phi(r, \theta, \phi) = \sum_{l=0}^{L-1}\sum_{m=0}^{l}\,\left[\Phi^{\cos}_{lm}(r)\,\cos(m\phi) + \Phi^{\sin}_{lm}(r)\,\sin(m\phi)\right]\,P_l^m(\cos\theta)\,,\]

where the radial potential functions are given by:

\[\Phi^{\cos,\sin}_{lm}(r) = -\frac{4\pi}{2l+1} \left[r^{-(l+1)} \int_0^r a^{l+2} \rho^{\cos,\sin}_{lm}(a) \, da + r^l \int_r^{\infty} a^{1-l} \rho^{\cos,\sin}_{lm}(a) \, da\right]\]

The spherical harmonic decomposition is performed via Gauss-Legendre quadrature (theta) and trapezoidal rule (phi). Radial integrals are evaluated to high precision using spline integration and precomputed cumulative integrals (I_inner, I_outer), which are stored as Bernstein polynomials to ensure smooth derivatives that satisfy the radial Poisson equation exactly. Outside of the radial grid, below the minimum radius, a constant-density extrapolation is used (with density equal to the value at the minimum grid radius), while above the maximum radius, the density is assumed to be zero.

Time-dependent potentials are supported in two ways:

• Via direct initialization: pass rho_cos_splines (and optionally rho_sin_splines) as lists of callables f(r, t) instead of InterpolatedUnivariateSpline instances, together with a tgrid array. Each callable should return the normalized radial density coefficient at the given radius and time (see the rho_cos_splines parameter documentation for details on the normalization convention).

• Via from_density: pass a density function that accepts a t keyword argument, i.e., dens(R, z, phi, t=0.0), together with a tgrid array of time grid points. The density function must be in galpy’s internal units (astropy Quantity output is not supported for time-dependent densities).

In both cases, the radial multipole integrals are precomputed at each time in tgrid, and the resulting radial functions \(\Phi^{\cos,\sin}_{lm}(r, t)\) are interpolated in time using piecewise polynomials. This allows efficient evaluation of the potential, forces, and second derivatives at arbitrary times within the tgrid range during orbit integration.

__init__(amp=1.0, rho_cos_splines=None, rho_sin_splines=None, rgrid=array([1.00000000e-03, 1.01036227e-03, 1.02083192e-03, ..., 2.93877957e+01, 2.96923201e+01, 3.00000000e+01], shape=(1001,)), tgrid=None, normalize=False, ro=None, vo=None)

Initialize a MultipoleExpansionPotential from precomputed density splines (use MultipoleExpansionPotential.from_density to initialize from a density function).

Parameters:

• amp (float or Quantity, optional) – Amplitude to be applied to the potential (default: 1).

• rho_cos_splines (list of list of InterpolatedUnivariateSpline, optional) – Cosine density coefficient splines with the spherical harmonic normalization absorbed, shape [L][M]. rho_cos_splines[l][m] is a spline representing \(\hat{\rho}^{\cos}_{lm}(r) = \beta_{lm}\,\rho^{\cos}_{lm}(r)\) as a function of r, where \(\rho^{\cos}_{lm}(r)\) are the coefficients in the real spherical-harmonic expansion of the density \(\rho(r,\theta,\phi) = \sum_{l,m} [\rho^{\cos}_{lm}(r)\,\cos(m\phi) + \rho^{\sin}_{lm}(r)\,\sin(m\phi)]\,P_l^m(\cos\theta)\) (the real-valued form of Eq. 12.78 in Bovy 2026) and \(\beta_{lm} = (2-\delta_{m0})\,\sqrt{\frac{2l+1}{4\pi}\,\frac{(l-m)!}{(l+m)!}}\) is the real spherical harmonic normalization factor. If None, computes a default Hernquist monopole. For time-dependent potentials, these can be callables f(r, t) instead of InterpolatedUnivariateSpline instances; in that case tgrid must also be provided.

• rho_sin_splines (list of list of InterpolatedUnivariateSpline, optional) – Like rho_cos_splines but for the sine coefficients: rho_sin_splines[l][m] represents \(\hat{\rho}^{\sin}_{lm}(r) = \beta_{lm}\,\rho^{\sin}_{lm}(r)\). If None, set to zero splines matching rho_cos_splines. For time-dependent potentials, these can be callables f(r, t) like rho_cos_splines.

• rgrid (numpy.ndarray, optional) – Radial grid points (1D array). Default: numpy.geomspace(1e-3, 30, 1_001).

• tgrid (numpy.ndarray or None, optional) – Time grid for time-dependent potentials. If provided, rho_cos_splines and rho_sin_splines are interpreted as callables f(r, t) and the potential is time-dependent: BPoly radial integrals are precomputed at each time in tgrid and interpolated in time at evaluation. Default: None (static potential).

• normalize (bool or float, optional) – If True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.

• ro (float or Quantity, optional) – Distance scale for translation into internal units (default from configuration file).

• vo (float or Quantity, optional) – Velocity scale for translation into internal units (default from configuration file).

Notes

• 2026-02-13 - Written - Bovy (UofT)

• 2026-03-06 - Refactored to accept splines; density computation moved to from_density - Bovy (UofT)

• 2026-03-23 - Added time-dependent support via tgrid - Bovy (UofT)

classmethod from_density(dens, L=6, rgrid=array([1.00000000e-03, 1.01036227e-03, 1.02083192e-03, ..., 2.93877957e+01, 2.96923201e+01, 3.00000000e+01], shape=(1001,)), tgrid=None, symmetry=None, costheta_order=None, phi_order=None, amp=1.0, normalize=False, ro=None, vo=None)

Initialize a MultipoleExpansionPotential from a density function.

Parameters:

• dens (callable or Potential) – Density function. Can take 1 arg (r), 2 args (R, z), or 3 args (R, z, phi). Can also be a galpy Potential instance. For time-dependent densities, add a t keyword argument (e.g., dens(R, z, phi, t=0.)).

• L (int, optional) – Maximum spherical harmonic degree + 1 (l goes from 0 to L-1). Default: 6.

• rgrid (numpy.ndarray, optional) – Radial grid points (1D array). Default: numpy.geomspace(1e-3, 30, 1_001).

• tgrid (numpy.ndarray or None, optional) – Time grid for time-dependent potentials. Required when the density function accepts a t parameter. Default: None.

• symmetry (str or None, optional) – 'spherical', 'axisymmetric', or None (general). Determines M.

• costheta_order (int, optional) – Gauss-Legendre quadrature order for theta. Default: max(20, L+1).

• phi_order (int, optional) – Number of uniform phi points for trapezoidal rule. Default: max(20, 2*L+1).

• amp (float or Quantity, optional) – Amplitude to be applied to the potential (default: 1).

• normalize (bool or float, optional) – If True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.

• ro (float or Quantity, optional) – Distance scale for translation into internal units (default from configuration file).

• vo (float or Quantity, optional) – Velocity scale for translation into internal units (default from configuration file).

Returns:

A new MultipoleExpansionPotential instance.

Return type:

MultipoleExpansionPotential

Notes

• 2026-03-06 - Written - Bovy (UofT)

• 2026-03-23 - Added time-dependent support - Bovy (UofT)