Disk potential using multipole expansion

class galpy.potential.DiskMultipoleExpansionPotential(amp=1.0, normalize=False, dens=<function DiskMultipoleExpansionPotential.<lambda>>, Sigma={'amp': 1.0, 'h': 0.3333333333333333, 'type': 'exp'}, hz={'h': 0.037037037037037035, 'type': 'exp'}, Sigma_amp=None, dSigmadR=None, d2SigmadR2=None, Hz=None, dHzdz=None, L=10, rgrid=array([1.00000000e-03, 1.01036227e-03, 1.02083192e-03, ...,        2.93877957e+01, 2.96923201e+01, 3.00000000e+01], shape=(1001,)), symmetry=None, costheta_order=None, phi_order=None, ro=None, vo=None)[source]

Class that implements a basis-function-expansion technique for solving the Poisson equation for disk (+halo) systems. We solve the Poisson equation for a given density \(\rho(R,\phi,z)\) by introducing K helper function pairs \([\Sigma_i(R),h_i(z)]\), with \(h_i(z) = \mathrm{d}^2 H(z) / \mathrm{d} z^2\) and search for solutions of the form

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

where \(r\) is the spherical radius \(r^2 = R^2+z^2\). We can solve for \(\Phi_{\mathrm{ME}}(R,\phi,z)\) by solving

\[\frac{\Delta \Phi_{\mathrm{ME}}(R,\phi,z)}{4\pi G} = \rho(R,\phi,z) - \sum_i\left\{ \Sigma_i(r)\,h_i(z) + \frac{\mathrm{d}^2 \Sigma_i(r)}{\mathrm{d} r^2}\,H_i(z)+\frac{2}{r}\,\frac{\mathrm{d} \Sigma_i(r)}{\mathrm{d} r}\left[H_i(z)+z\,\frac{\mathrm{d}H_i(z)}{\mathrm{d} z}\right]\right\}\,.\]

We solve this equation by using the MultipoleExpansionPotential class. This technique works very well if the disk portion of the potential can be exactly written as \(\rho_{\mathrm{disk}} = \sum_i \Sigma_i(R)\,h_i(z)\), because the effective density on the right-hand side of this new Poisson equation is then not ‘disky’ and can be well represented using spherical harmonics. But the technique is general and can be used to compute the potential of any disk+halo potential; the closer the disk is to \(\rho_{\mathrm{disk}} \approx \sum_i \Sigma_i(R)\,h_i(z)\), the better the technique works.

This technique was introduced by Kuijken & Dubinski (1995) and was popularized by Dehnen & Binney (1998). The current implementation is a slight generalization of the technique in those papers and uses the MultipoleExpansionPotential to solve the Poisson equation for \(\Phi_{\mathrm{ME}}(R,\phi,z)\).

__init__(amp=1.0, normalize=False, dens=<function DiskMultipoleExpansionPotential.<lambda>>, Sigma={'amp': 1.0, 'h': 0.3333333333333333, 'type': 'exp'}, hz={'h': 0.037037037037037035, 'type': 'exp'}, Sigma_amp=None, dSigmadR=None, d2SigmadR2=None, Hz=None, dHzdz=None, L=10, rgrid=array([1.00000000e-03, 1.01036227e-03, 1.02083192e-03, ...,        2.93877957e+01, 2.96923201e+01, 3.00000000e+01], shape=(1001,)), symmetry=None, costheta_order=None, phi_order=None, ro=None, vo=None)[source]

Initialize a DiskMultipoleExpansionPotential.

Parameters:

  • amp (float, optional) – Amplitude to be applied to the potential (default: 1); cannot have units currently.

  • 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.

  • dens (callable) – Function of R,z[,phi optional] that gives the density [in natural units, cannot return a Quantity currently].

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

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

  • symmetry (str or None, optional) – 'spherical', 'axisymmetric', or None (general). If None and the density is axisymmetric, 'axisymmetric' is used automatically.

  • costheta_order (int, optional) – Gauss-Legendre quadrature order for theta.

  • phi_order (int, optional) – Number of uniform phi points for trapezoidal rule.

  • Sigma (dict or callable) – Either a dictionary of surface density (example: {‘type’:’exp’,’h’:1./3.,’amp’:1.,’Rhole’:0.} for amp x exp(-Rhole/R-R/h) ) or a function of R that gives the surface density.

  • hz (dict or callable) – Either a dictionary of vertical profile, either ‘exp’ or ‘sech2’ (example {‘type’:’exp’,’h’:1./27.} for exp(-|z|/h)/[2h], sech2 is sech^2(z/[2h])/[4h]) or a function of z that gives the vertical profile.

  • Sigma_amp (float, optional) – Amplitude to apply to all Sigma functions.

  • dSigmadR (callable, optional) – Function that gives d Sigma / d R.

  • d2SigmadR2 (callable, optional) – Function that gives d^2 Sigma / d R^2.

  • Hz (callable, optional) – Function of z such that d^2 Hz(z) / d z^2 = hz.

  • dHzdz (callable, optional) – Function of z that gives d Hz(z) / d z.

  • 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

  • Either specify (Sigma,hz) or (Sigma_amp,Sigma,dSigmadR,d2SigmadR2,hz,Hz,dHzdz)

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