Disk potential using SCF basis-function-expansion

class galpy.potential.DiskSCFPotential(amp=1.0, normalize=False, dens=<function DiskSCFPotential.<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, N=10, L=10, a=1.0, radial_order=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 SCFPotential class and methods (scf_compute_coeffs_axi or scf_compute_coeffs depending on whether \(\rho(R,\phi,z)\) is axisymmetric or not). 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 SCF approach of Hernquist & Ostriker (1992) to solve the Poisson equation for \(\Phi_{\mathrm{ME}}(R,\phi,z)\) rather than solving it on a grid using spherical harmonics and interpolating the solution (as done in Dehnen & Binney 1998).

__init__(amp=1.0, normalize=False, dens=<function DiskSCFPotential.<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, N=10, L=10, a=1.0, radial_order=None, costheta_order=None, phi_order=None, ro=None, vo=None)[source]

NAME:

__init__

PURPOSE:

initialize a DiskSCF Potential

INPUT:

amp - amplitude to be applied to the potential (default: 1); cannot have units currently

normalize - 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=, vo= distance and velocity scales for translation into internal units (default from configuration file)

dens= function of R,z[,phi optional] that gives the density [in natural units, cannot return a Quantity currently]

N=, L=, a=, radial_order=, costheta_order=, phi_order= keywords setting parameters for SCF solution for Phi_ME (see scf_compute_coeffs_axi or scf_compute_coeffs depending on whether \(\rho(R,\phi,z)\) is axisymmetric or not)

Either:

  1. Sigma= Dictionary of surface density (example: {‘type’:’exp’,’h’:1./3.,’amp’:1.,’Rhole’:0.} for amp x exp(-Rhole/R-R/h) )

    hz= 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])

  2. Sigma= function of R that gives the surface density

    dSigmadR= function that gives d Sigma / d R

    d2SigmadR2= function that gives d^2 Sigma / d R^2

    Sigma_amp= amplitude to apply to all Sigma functions

    hz= function of z that gives the vertical profile

    Hz= function of z such that d^2 Hz(z) / d z^2 = hz

    dHzdz= function of z that gives d Hz(z) / d z

In both of these cases lists of arguments can be given for multiple disk components; can’t mix (a) and (b) in these lists; if hz is a single item the same vertical profile is assumed for all Sigma

OUTPUT:

DiskSCFPotential object

HISTORY:

2016-12-26 - Written - Bovy (UofT)