vcdisk.vcdisk

vcdisk.vcdisk(rad, sb, z0=0.3, rhoz='cosh', rhoz_args=None, flaring=False, zsamp='log', rsamp='log')

Circular velocity of a thick disk of arbitrary surface density.

This function uses the method of [Casertano83] to calculate the radial force on the disk plane and then the circular velocity.

Parameters

  • rad (list or numpy.ndarray) – array of radii in \(\rm kpc\).

  • sb (list or numpy.ndarray) – array of surface densities in \(\rm M_\odot / kpc^2\).

  • z0 (float, optional) – disk scaleheight in kpc. Default: 0.3 \(\rm kpc\).

  • rhoz (str or callable, optional) – vertical density function. Can either be one of two hardcoded options, 'cosh' (default) or 'exp', or it can be any user-defined function with input and output numpy.ndarray. The function should define the vertical surface density in \(\rm M_\odot / kpc^2\) and it should be normalized such that rhoz(0)=1. It can have additional arguments handled by rhoz_args. See Vertical density for details.

  • rhoz_args (dict, optional) – dictionary of arguments of the user-defined function rhoz.

  • flaring (bool, optional) – whether the vertical density function rhoz depends also on radius, i.e. if \(\rho_z=\rho_z(z,R)\). If True, then rhoz must be a callable function with signature rhoz(z, R) or rhoz(z, R, **rhoz_args). This option assumes that the vertical density provided is normalized, i.e. \(\int {\rm d}z\, \rho_z(z,R)=1\,\,\,\forall R\).

  • zsamp (str or numpy.ndarray, optional) – sampling in the z-direction. Can be either 'log' (default) for logarithmically spaced values, 'lin' for linearly spaced values, or a user-defined np.array. If 'log' or 'lin' an numpy.ndarray is created in the range [z0/10, z0*10] with size 100.

  • rsamp (str, optional) – sampling in the \(R\)-direction. Can be either 'log' (default) for logarithmically spaced values, 'lin' for linearly spaced values, or 'nat' for the natural spacing of the data in input rad.

Returns

array of \(V_{\rm disk}\) velocities in \(\rm km/s\).

Return type

numpy.ndarray

See also

gipsy.rotmod, galpynamics

Notes

Background

The circular velocity on the mid-plane of a disk galaxy can be written as

\[V_{\rm disk}(r) = \sqrt{-r\,F_{r,\rm disk}(r)},\]

where \(F_{r,\rm disk}(r)\) is the radial force on the plane of the disk, which is calculated as

\[F_{r,\rm disk}(r) = 4\pi G \int_0^\infty \mathrm{d}u \int_0^\infty \mathrm{d}z \,\,2\sqrt{\frac{u}{rp}} \,\frac{\mathcal{K}(p)-\mathcal{E}(p)}{\pi}\, \frac{\partial \rho(u,z)}{\partial u}.\]

Here \(p = x-\sqrt{x^2-a}\), \(x = (r^2+u^2+z^2)/(2ru)\), and \(\mathcal{K}\) and \(\mathcal{E}\) are the complete elliptic integrals respectively of the first and second kinds (see Eqs. 4-5-6 and Appendix A of [Casertano83] for full details).

This program follows the notation of gipsy in defining a signed circular velocity which has the opposite sign of the radial force, i.e.

\[V_{\rm disk}(r) = -{\rm sign}\left[F_{r,\rm disk}(r)\right]\sqrt{r\,|F_{r,\rm disk}(r)|}.\]

This convention is such that \(V_{\rm disk}\) is negative if there is a net force away from the galaxy center. This can happen, for instance, in cases where the surface density has a depression close to the center, as it is often observed in the neutral hydrogen disks of nearby spiral galaxies.

Vertical density

The radial force on the disk plane depends both on the radial surface density as well as on the vertical density distribution of the disk. While the radial surface density of can be typically obtained from observations of the surface brightness, in rotation curve studies the vertical density distribution is often unknown.

Here we provide the user with several choices for the vertical density profile of the disk.

Vertical density: constant scaleheight

The simplest and the most common assumption is that the axisymmetric disk has a constant scaleheight, which implies that the density is separable as \(\rho(R,z)=\rho_R(R)\rho_z(z)\). Thus, defining the surface density on the disk plane as \(\Sigma(R)=\int {\rm d}z\, \rho(R,z)\), we can write \(\rho(R,z)=\Sigma(R)\rho_z(z)/N\), where \(N\) is the normalization of the vertical density \(N=\int {\rm d}z\, \rho_z(z)\).

We provide full implementations of two popular choices for \(\rho_z(z)\) that are often used to describe disk galaxies, and we also provide an interface for the user to supply their own vertical density profile.

  • 'cosh' (default): this corresponds to \(\rho_z(z)=\cosh^{-2}(z/z_0)\), with \(N=2z_0\), i.e. the so-called [vdKruitSearle81] disk,

  • 'exp': this corresponds to \(\rho_z(z)=\exp(-z/z_0)\), with \(N=2z_0\).

  • the user may also choose to specify a custom vertical density profile by passing a callable function as the rhoz argument.

Note

If the user specifies a callable rhoz, vcdisk.vcdisk() will try to normalize it for you using scipy.integrate.quad() to solve the integral \(\int {\rm d}z\, \rho_z(z)\). Do look out for potential issues related to this normalization.

Warning

If the user-defined function rhoz has signature rhoz(z, **rhoz_args), then make sure that the rhoz_args dictionary supplied to vcdisk.vcdisk() is in the correct positional order, since keys get lost when calculating the normalization integral.

Vertical density: flaring disks

Another possibility is that the vertical density \(\rho_z\) depends also on radius, so that \(\rho(R,z)\) is not separable. Often this happens since galactic disks flare in the outer regions, where the gravitational support is lower, thus one may wish to compute the circular velocity in the case where the scaleheight depends on radius, \(z_0=z_0(R)\).

By activating the option flaring=True, vcdisk.vcdisk() computes the potential with a \(\rho_z(z,R)\) specified via the rhoz parameter. As an example, this could be something like \(\rho_z(z,R) = e^{-z/z_0(R)}/2z_0(R)\).

Note

If flaring=True, vcdisk.vcdisk() expects the callable rhoz function to be already normalized, such that \(\int {\rm d}z\, \rho_z(z,R)=1\,\,\,\forall R\).

Implementation details

vcdisk.vcdisk() uses Simpson’s method, implemented in scipy.integrate.simpson(), to compute the double integral quickly and efficiently. The scipy library is also used for the implementation of the elliptic integrals scipy.special.ellipk() and scipy.special.ellipe().

The arbitrary sampling in the \(z\)-direction unfortunately seems to somewhat have a (small) impact on the derived \(V_{\rm disk}\). This is probably because of numerical noise in the double integral when Simpson’s method is used on very steep functions such as 3-D densities of galaxy disks. After extensive tests where I compared the results with analytic formulae for a thin disk and with gipsy’s implementation, I have come to the conclusion that a logarithmic sampling in the \(z\)-direction in the range \(z\in[z_0/10, 10\,z_0]\) with 100 samples is a good compromise between speed and accuracy.

References

Casertano83(1,2)

Casertano, S. 1983, MNRAS, 203, 735. Rotation curve of the edge-on spiral galaxy NGC 5907: disc and halo masses. doi:10.1093/mnras/203.3.735

vdKruitSearle81

van der Kruit, P. C. & Searle, L. 1981, A&A, 95, 105. Surface photometry of edge-on spiral galaxies. I - A model for the three-dimensional distribution of light in galactic disks.

Example

>>> import numpy as np
>>> from vcdisk import vcdisk
>>> md, rd = 1e10, 2.0                        # mass, scalelength of the disk
>>> r = np.linspace(0.1, 30.0, 50)            # radii samples
>>> sb = md / (2*np.pi*rd**2) * np.exp(-r/rd) # exponential disk surface density
>>> vcdisk(r, sb)
array([ 6.36888305, 38.63854589, ..., 38.26020293, 37.82648552])