Local velocity distribution#

pianistpy.pdf_v(), pianistpy.LosTable

pdf_v(source, r, v) returns the local velocity distribution \(p(V\,|\,R)\) at user-specified 3D radii. \(p(V\,|\,R)\) is the probability density of the velocity \(V = |\mathbf{V}|\) of the model’s mass at the distance \(R\) from the center. As in Chiang+(2026d), the density at a radius is written as an integral over the angle \(\varphi\) between the velocity and the radial direction, with the angle density \(p_\varphi\), and over the velocity magnitude, with its density \(p_V\),

\[\varrho(R) = 2\pi \int_0^\pi {\rm d}\varphi\; p_\varphi(\varphi\,|\,R) \int_0^{V_{\rm esc}(R)} {\rm d}V\; p_V(V\,|\,R, \varphi), \qquad V_{\rm esc}(R) = \sqrt{2P(R)},\]

and the velocity distribution is the product \(p_\varphi\,p_V\) integrated over the angle and multiplied by \(2\pi/\varrho(R)\):

\[p(V\,|\,R) = \frac{2\pi}{\varrho(R)} \int_0^\pi p_\varphi(\varphi\,|\,R)\, p_V(V\,|\,R, \varphi)\,{\rm d}\varphi .\]

\(p(V\,|\,R)\) integrates to one over \(V \geq 0\). For an isotropic, constant-anisotropy or Type I model, the support of \(p(V\,|\,R)\) ends at the escape velocity \(V_{\rm esc}(R)\). A Type II model also holds unbound orbits, and the velocities of a Type II model reach \(V_{\rm esc}/\sqrt{1 - R^2/R_{\rm a}^2}\). For a Type II model, a radius at or beyond the anisotropy radius \(R_{\rm a}\) raises InputError. The radius \(R\) of Chiang+(2026d) is r in the package.

Arguments and units#

  • source is a Model or the model’s line-of-sight table (LosTable).

  • r and v are arrays that broadcast against each other. r holds radii in the model’s unit of length: internal units by default, parsecs or an Astropy quantity in the physical system. v holds velocities in the model’s unit of velocity (km/s in the physical system).

  • The result has the broadcast shape of r and v and is a probability density per unit velocity of the model’s unit system. The result is zero at negative velocities and beyond the end of the support.

For a superposition, the components are mixed by their density at the radius, \(p = \sum_k w_k \varrho_k p_k / \sum_k w_k \varrho_k\).

Accuracy#

Each distinct radius is computed once, with three checks:

  • the norm of the distribution against the density;

  • the second moment of the distribution against \(\sigma_{\rm r}^2 (3 - 2\beta_\sigma)\) from Radial velocity dispersion;

  • the change of the distribution when the number of nodes of the distribution’s angle integral doubles.

The number of nodes of the angle integral doubles, at most three times, until the change of the distribution is under velocity_tolerance: \(10^{-6}\) with the default settings, \(10^{-4}\) with Settings(quick=True). A distribution whose relative error is above \(10^{-4}\) is returned with an AccuracyWarning. A distribution whose relative error is above \(10^{-2}\) raises AccuracyError, unless the model’s settings set write_anyway=True.

pianistpy.observables.stages.velocity_distribution() returns the velocity distribution of a one-component model at one radius, in internal units in every unit system. The returned object also holds the record of the distribution’s checks, the escape velocity, the end of the distribution’s support and the share of the distribution above the escape velocity (unbound_fraction, nonzero only for a Type II model).

Example#

import numpy as np
import pianistpy as pp

model = pp.solve(
    pp.Component(pp.DoublePowerLaw.hernquist(), pp.Anisotropy(beta_c=0.25)),
    pp.Truncation(k_t=0.05),
    settings=pp.Settings(quick=True),
)
v = np.linspace(0.0, 1.4, 1401)                           # beyond the escape velocity at each of the three radii
p = pp.pdf_v(model, np.array([[0.1], [1.0], [10.0]]), v)   # (3, 1401): radii by velocities
p.sum(axis=1) * (v[1] - v[0])                              # one at each radius, to the accuracy of the sum on this grid

physical = model.to("physical", concentration=10.0, mass_0=1e9, r_0=500.0)
pp.pdf_v(physical, 500.0, [10.0, 50.0, 100.0])             # per km/s, at 500 pc