---
title: 'pybhpt: Python Library for Kerr Perturbations'
url: https://www.emergentmind.com/topics/pybhpt
type: topic
---

# pybhpt: Python Library for Kerr Perturbations

Searching arXiv for recent papers on pybhpt and related black-hole perturbation theory.
pybhpt is an open, frequency-domain Python library that implements first-order, in the small mass ratio $\epsilon = m_2/m_1$, black-hole perturbation theory for a point particle on generic eccentric, precessing bound orbits in Kerr spacetime. It solves the Teukolsky equation for the maximal-spin Newman–Penrose curvature scalars $\psi_0$ and $\psi_4$, reconstructs the metric perturbation in multiple radiation-gauge variants, computes the generalized redshift invariant, and, through a recent Hamiltonian formulation, supplies the conservative first-order Hamiltonian needed for first post-adiabatic waveform generation. In later work, pybhpt also serves as a benchmark and integration target for alternative frequency-domain Kerr solvers based on a unified confluent-Heun framework [2507.07746; 2605.09250].

## 1. Problem domain and physical content

pybhpt is designed for first-order self-force and waveform calculations in the small-mass-ratio limit for a point particle on generic eccentric, precessing, bound Kerr geodesics. The library computes frequency-domain solutions for the Weyl scalars $\psi_0$ with spin weight $s=+2$ and $\psi_4$ with spin weight $s=-2$, assembles both horizon-side and infinity-side extended homogeneous solutions, reconstructs radiation-gauge metric perturbations $h_{\mu\nu}$, and evaluates the generalized redshift invariant along eccentric, precessing Kerr orbits. Through the identification of the interaction Hamiltonian with the time-averaged first-order redshift correction, it provides conservative Hamiltonian data for $1\mathrm{PA}$ waveforms with conservative self-force effects. The code reaches near-extremal spins, with $a/M$ up to $0.999$, eccentricities up to $e \approx 0.6$, and inclinations up to $49\pi/100$ [2507.07746].

The implemented quantities sit at the interface of several relativistic-perturbation problems. On the dissipative side, pybhpt computes fluxes at infinity and at the horizon from $\psi_4$ modes. On the conservative side, it computes $\langle \tilde z_1\rangle_t$, which directly determines the first-order conservative Hamiltonian. This places the library in a position to connect self-force calculations with effective-one-body and post-Newtonian descriptions. A plausible implication is that pybhpt is useful not only for flux production but also for invariant-based calibration problems across approximation schemes.

## 2. Frequency-domain Kerr perturbation framework

pybhpt adopts a null tetrad aligned with Kerr’s principal null directions and, in the Geroch–Held–Penrose formalism, uses the Kinnersley tetrad. At linear order, the gauge-invariant curvature information is encoded by the maximal spin-weight scalars $\psi_0$ and $\psi_4$, related to the metric perturbation by linear differential operators. The rescaled scalars satisfy the Teukolsky equations with spin $s=\pm 2$, sourced by the point-particle stress-energy tensor.

For bound motion, the frequency spectrum is discrete:
$$
\omega \equiv \omega_{mkn} = m \Omega_\phi + k \Omega_\theta + n \Omega_r,
$$
with integers $m$, $k$, and $n$. For equatorial motion, $k=0$; for circular or spherical motion, $n=0$. The radial Teukolsky equation is
$$
\left[\Delta^{-s}\frac{d}{dr}\left(\Delta^{s+1}\frac{d}{dr}\right) + V_{sjm\omega}(r)\right]\psi^{\mathrm{ret}}_{\pm 2\,jm\omega}(r) = T_{\pm 2\,jm\omega}(r),
$$
with
$$
V_{sjm\omega}(r) = \frac{K^2 - 2is(r-M)K}{\Delta} + 4is\omega r - \lambda_{sjm\omega}, \qquad
K=(r^2+a^2)\omega-ma.
$$

Retarded boundary conditions are imposed at the horizon and at null infinity using homogeneous solutions $R^{\mathcal H}_{sjm\omega}$ and $R^{\mathcal I}_{sjm\omega}$ normalized by their asymptotics. The amplitudes are computed through Green’s-function integrals,
$$
\psi^{\mathrm{ret}}_{\pm 2\,jm\omega}(r \le r_{\min}) = Z^{\mathcal H}_{\pm 2\,jm\omega} R^{\mathcal H}_{\pm 2\,jm\omega}, \qquad
\psi^{\mathrm{ret}}_{\pm 2\,jm\omega}(r \ge r_{\max}) = Z^{\mathcal I}_{\pm 2\,jm\omega} R^{\mathcal I}_{\pm 2\,jm\omega},
$$
with
$$
Z^{\mathcal H/\mathcal I}_{\pm 2\,jm\omega}
=
\int_{r_{\min}}^{r_{\max}}
\frac{\Delta^{\pm 2} R^{\mathcal I/\mathcal H}_{\pm 2\,jm\omega} T_{\pm 2\,jm\omega}}
{\mathcal W_{\pm 2\,jm\omega}}\,dr,
$$
and
$$
\mathcal W_{sjm\omega}
=
\Delta^{s+1}
\left(
R^{\mathcal H}_{sjm\omega}\frac{dR^{\mathcal I}_{sjm\omega}}{dr}
-
R^{\mathcal I}_{sjm\omega}\frac{dR^{\mathcal H}_{sjm\omega}}{dr}
\right).
$$
Extended homogeneous solutions are then assembled as
$$
\psi^{\mathcal J}_{\pm 2\,jm\omega}(r) = Z^{\mathcal J}_{\pm 2\,jm\omega} R^{\mathcal J}_{\pm 2\,jm\omega}(r).
$$
These formulas define the core frequency-domain structure on which the remainder of the library is built [2507.07746].

The orbital sector is parameterized either by $(\mathcal E,\mathcal L_z,\mathcal Q)$ or by orbital elements $(p,e,x\equiv \cos\iota)$. Introducing Mino time $\lambda$ via $\Sigma\,d\lambda=d\tau$, the motion separates into librations in $r$ and $\theta$ and linear drifts in $t$ and $\phi$, with fundamental coordinate-time frequencies
$$
\Omega_r = \Upsilon_r/\Upsilon_t,\qquad
\Omega_\theta = \Upsilon_z/\Upsilon_t,\qquad
\Omega_\phi = \Upsilon_\phi/\Upsilon_t.
$$
This structure is the basis for the spectral source integration used later in the code.

## 3. Metric reconstruction and gauge structure

In vacuum regions, pybhpt reconstructs the metric perturbation from Hertz potentials using Wald’s adjoint-operator identity. The library implements four reconstruction formalisms. Two are the Chrzanowski–Cohen–Kegeles–Wald reconstructions from either $\psi_0$ or $\psi_4$, leading to ingoing radiation gauge and outgoing radiation gauge. Two build on Aksteiner–Andersson–Bäckdahl identities and use both $\psi_0$ and $\psi_4$; the paper labels these symmetric radiation gauge and antisymmetric radiation gauge.

The reconstructed fields are
$$
h^{\mathrm{ORG,rec}}_{\alpha\beta}
=
4\,\mathrm{Re}\,(\mathcal S_4^\dagger \zeta^4 \Phi_{O0})_{\alpha\beta},
$$
$$
h^{\mathrm{IRG,rec}}_{\alpha\beta}
=
4\,\mathrm{Re}\,(\mathcal S_0^\dagger \zeta^4 \Phi_{I4})_{\alpha\beta},
$$
$$
h^{\mathrm{SRG,rec}}_{\alpha\beta}
=
4\,\mathrm{Re}\,\big[\mathcal S_4^\dagger \zeta^4 \Phi_{S0} + \mathcal S_0^\dagger \zeta^4 \Phi_{S4}\big]_{\alpha\beta},
$$
$$
h^{\mathrm{ARG,rec}}_{\alpha\beta}
=
4\,\mathrm{Re}\,\big[\mathcal S_4^\dagger \zeta^4 \Phi_{A0} - \mathcal S_0^\dagger \zeta^4 \Phi_{A4}\big]_{\alpha\beta}.
$$
IRG satisfies $l^\alpha h^{\mathrm{IRG}}_{\alpha\beta}=0$ and $g^{\alpha\beta}h^{\mathrm{IRG}}_{\alpha\beta}=0$; ORG satisfies $n^\alpha h^{\mathrm{ORG}}_{\alpha\beta}=0$ and $g^{\alpha\beta}h^{\mathrm{ORG}}_{\alpha\beta}=0$. SRG and ARG are traceless combinations rather than standard gauges. The Aksteiner–Andersson–Bäckdahl identity gives
$$
M \mathcal L_\xi h^{\mathrm{SRG,rec}}_{\alpha\beta}
=
\frac{4}{3}\,\mathrm{Re}\,
\big[\mathcal S_4^\dagger \zeta^4 \psi_0 - \mathcal S_0^\dagger \zeta^4 \psi_4\big]_{\alpha\beta},
$$
with $\mathcal L_\xi = \partial_t$ in Boyer–Lindquist coordinates. ARG requires time derivatives of Hertz potentials, and static modes with $\omega=0$ are excluded from ARG in this work [2507.07746].

For point particles, the spacetime is divided into two vacuum domains separated by the worldline. Reconstruction is performed on each side, and radiation-gauge reconstructions are glued into a no-string solution. Physical solutions also require stationary-axisymmetric completion terms,
$$
h^{\mathrm{comp}\,\pm}_{\alpha\beta}
=
c_M^\pm \frac{\partial g_{\alpha\beta}}{\partial M}
+
c_J^\pm \frac{\partial g_{\alpha\beta}}{\partial J}.
$$
For bound geodesics,
$$
c_M^- = c_J^- = 0,\qquad
c_M^+ = M\mathcal E,\qquad
c_J^+ = M^2\mathcal L_z.
$$
The completed solution in each domain and in a target gauge $G'$ is
$$
h^{G'\pm}_{\alpha\beta}
=
h^{G,\mathrm{rec}\,\pm}_{\alpha\beta}
+
x_{\alpha\beta}
+
\mathcal L_{\xi^\pm} g_{\alpha\beta}.
$$
This organization is central to pybhpt’s redshift extraction, because continuity of $h_{uu}=h_{\alpha\beta}u^\alpha u^\beta$ across the worldline is the key practical condition for obtaining the regularized invariant.

## 4. Generalized redshift invariant and Hamiltonian interpretation

The generalized redshift invariant is defined in the effective metric $\tilde g_{\mu\nu}=g_{\mu\nu}+\epsilon h^{R(1)}_{\mu\nu}+O(\epsilon^2)$ by $\tilde z = 1/\tilde u^t$. Expanding,
$$
\tilde z = \tilde z_0 + \epsilon \tilde z_1 + O(\epsilon^2),\qquad
\tilde z_1 = -\frac{1}{2}\tilde z_0 h^{R}_{uu},\qquad
\tilde z_0 = (u^t)^{-1}.
$$
The generalized redshift is the infinite coordinate-time average at fixed frequencies $\Omega_i=(\Omega_r,\Omega_\theta,\Omega_\phi)$, primary mass $m_1$, secondary mass $m_2$, and spin $a$:
$$
\langle \tilde z\rangle_t(\Omega_i;m_1,m_2,a)=\left\langle \frac{d\tilde\tau}{dt}\right\rangle_t.
$$
For bound biperiodic motion, the averaging is written as phase-space integrals in $(q_r,q_z)$:
$$
\langle \tilde z^{(0)}\rangle_t
=
\Upsilon_t^{-1}
\int_0^{2\pi}\int_0^{2\pi}
\Sigma_p(q_r,q_z)\,
\frac{dq_z}{2\pi}\frac{dq_r}{2\pi},
$$
$$
\langle \tilde z^{(0)} h^R_{uu}\rangle_t
=
\Upsilon_t^{-1}
\int_0^{2\pi}\int_0^{2\pi}
\Sigma_p h^R_{uu}\,
\frac{dq_z}{2\pi}\frac{dq_r}{2\pi}.
$$
The resulting quantity is described as a quasi-invariant across physically reasonable gauges [2507.07746].

Regularization uses a mode-sum compatible with locally Lorenz singular structure:
$$
h^R_{uu} = \sum_{\ell=0}^\infty \big(h^\ell_{uu} - h^{S,\ell}_{uu}\big).
$$
In practice, the leading locally Lorenz parameter $H_{[0]}$ suffices to render the sum convergent, while pybhpt accelerates convergence by fitting the large-$\ell$ tail to determine effective higher-order $H_{[2n]}$. For generic Kerr orbits,
$$
H_{[0]} = \frac{4}{\pi}\sqrt{\eta/k}\,K(k),\qquad
k=\frac{2\eta}{\eta+\zeta},
$$
with $\eta^2$ and $\zeta$ given in terms of Kerr metric functions and constants of motion, and $K(k)$ the complete elliptic integral.

The Hamiltonian connection is
$$
H(J_i,\phi_i)
=
E^{(0)}(J_i)
+
\frac{\epsilon}{2}\langle \mathcal H^{(1)}\rangle_t(J_i)
+
O(\epsilon^2),
$$
with actions $J_i=\{J_r,J_\theta,J_\phi\}$, and
$$
\langle \mathcal H^{(1)}\rangle_t = m_2\big[\langle \tilde z_1\rangle_t + O(\epsilon)\big].
$$
Thus pybhpt provides $\langle \tilde z_1\rangle_t(J_i)$, from which conservative equations of motion follow through Hamilton’s equations. This is the library’s main link to waveform generation with conservative $1\mathrm{PA}$ effects and to cross-framework comparisons involving self-force, EOB, and PN descriptions.

## 5. Software organization and computational workflow

pybhpt is organized into domain-specific modules that follow the perturbative pipeline from geodesic motion to radiative and conservative observables [2507.07746].

| Module | Role |
|---|---|
| `geo` | Geodesic functions, turning points, and fundamental frequencies; spectral evaluation of $\Delta \hat x$ and $\Upsilon_i$ |
| `swsh` | Spin-weighted spheroidal harmonics via spherical expansions; coupling coefficients $b_{\ell sjm\omega}$ and eigenvalues $\lambda_{sjm\omega}$ |
| `radial` | Homogeneous Teukolsky ODEs with Zenginoğlu stabilization; Teukolsky–Starobinsky mappings between $s$ and $-s$; monodromy-based $\nu$ optional |
| `teuk` | Inhomogeneous amplitudes $Z^{\mathcal J}_{\pm 2\,jm\omega}$ via SSI; robust error tracking and thresholds |
| `hertz` | Mode functions $\Phi^{\mathcal J}_{\pm 2}$ using identities; angular and radial derivatives for reconstruction |
| `metric` | Analytic reconstruction coefficients for applying GHP operators to Hertz modes; assembling $h^{\mathrm{rec}}_{\alpha\beta}$ |
| `redshift` | Coefficients to form $h_{uu}$, worldline projection, mode-sum assembly, and regularization/fitting to extract $\langle \tilde z_1\rangle_t$ |
| `flux` | GW fluxes at infinity and horizon built from $\psi_4$ modes |

The library’s typical workflow is geodesics $\rightarrow$ $\psi_0/\psi_4$ modes $\rightarrow$ Hertz potentials $\rightarrow$ metric reconstruction $\rightarrow$ $h_{uu}$ $\rightarrow$ $\langle \tilde z_1\rangle_t$ $\rightarrow$ $\langle \mathcal H^{(1)}\rangle_t$. In the angular sector, spin-weighted spheroidal harmonics are expanded in spin-weighted spherical harmonics,
$$
{}_{s}S_{jm\omega}(z)
=
\sum_{\ell=\mathrm{basis}} b_{\ell sjm\omega}\,{}_{s}Y_{\ell m}(z),
$$
with the coupling coefficients obtained from a five-term recurrence recast as an eigenvalue problem. In the radial sector, homogeneous solutions are integrated using Zenginoğlu’s transformation and a Prince–Dormand $(8,9)$ method with adaptive stepping from the GNU Scientific Library. Coulomb or Heun boundary series supply initial data, and the complementary spin sector is rebuilt through reduced Teukolsky–Starobinsky identities.

The source amplitudes are evaluated by spectral source integration. The source takes a distributional form with up to two $r$-derivatives,
$$
T_{\pm 2\,jm\omega}(r)
=
\epsilon \int dt \,\frac{dt}{u^t}\,\Delta^{\mp 2}
J_{\pm 2\,jm\omega}(t,r)\,e^{i(\omega t-m\phi_p(t))},
$$
where
$$
J_{\pm 2\,jm\omega}
=
J^{(0)}_{\pm 2\,jm\omega}\delta(r-r_p)
+
\partial_r\!\big[J^{(1)}_{\pm 2\,jm\omega}\delta(r-r_p)\big]
+
\partial_r^2\!\big[J^{(2)}_{\pm 2\,jm\omega}\delta(r-r_p)\big].
$$
Hertz potentials are then built in mode-sum form, reconstruction operators are applied in a spherical-harmonic basis, and the worldline quantity $h_{uu}$ is regularized and averaged.

## 6. Convergence, validation, and limitations

Away from the worldline, convergence is exponential. Near the worldline, the radiation gauges have the expected singular structure. IRG, ORG, and SRG have the expected $O(s^{-1})$ singularity near the particle, producing $h_\ell \sim \ell^0$ as $\ell\to\infty$. ARG modes are more singular, with $\sim \ell^2$–$\ell^3$ behavior in half-string form, reflecting an $O(s^{-3}$–$s^{-4})$ divergence, and ARG excludes static modes with $\omega=0$. SRG inherits the “worst” asymptotic behavior of IRG and ORG at boundaries but remains $O(s^{-1})$ near the particle [2507.07746].

The main numerical limitation identified in the paper comes from spectral source integration at very high mode numbers. For large $m$, $k$, and $n$, especially at high eccentricity and small $p$, oscillatory cancellation limits accuracy; pybhpt drops modes when the relative error exceeds a threshold. This limits accuracy for $e \gtrsim 0.5$ near the innermost stable orbit. The paper reports that eccentric equatorial redshift values agree with van de Meent–Shah (2015) at the $10^{-6}$–$10^{-5}$ level across gauges. It also reports new precessing results with spins up to $a/M=0.999$ and, for near-extremal Kerr, the first observation of negative $\langle \tilde z_1\rangle_t$ near the ISCO, implying a surface in $(p,e,x)$ where the interaction Hamiltonian vanishes.

These features are important for interpretation. The gauge-dependent local fields remain singular in the expected ways, but the regularized redshift invariant and Hamiltonian data retain their utility. A plausible implication is that pybhpt’s most robust outputs are invariant or quasi-invariant quantities assembled after reconstruction and regularization, rather than raw local metric components near the particle.

## 7. Relation to later HeunC-based solvers

A later study on generic Kerr-orbit fluxes describes pybhpt as implementing a hybrid semi-analytical/semi-numerical pipeline: MST series for the radial Teukolsky equation, spectral methods for spin-weighted spheroidal harmonics, and Frobenius expansions of HeunC for boundary data, together with adaptive Runge–Kutta integration and spectral source techniques. That work identifies three typical challenges for pybhpt-style frequency-domain calculations: computational overhead in determining auxiliary parameters such as MST’s renormalized angular momentum $\nu$ and high-accuracy angular eigenvalues, stiffness and loss of accuracy for strong-field or high-frequency modes, and sensitivity of oscillatory source integrals to grid resolution [2605.09250].

The unified confluent-Heun framework reformulates both the angular and radial Teukolsky equations directly as confluent Heun equations and computes global solutions via Motygin’s hybrid analytic-continuation algorithm. In the benchmarks reported there, for the total radiative flux summed over $168$ low-order modes, the HeunC framework achieves relative errors of order $10^{-11}$, with HeunC runtime $9$–$11$ s versus pybhpt $39$–$126$ s, corresponding to a $3$–$13\times$ speedup. For a highly oscillatory single-mode case, HeunC converges in $247$ ms, while pybhpt’s uniform trapezoidal rule over $[0,2\pi]$ took $\sim 790$ ms with larger error. The same study describes integration into existing pipelines, including pybhpt, as straightforward: reuse orbit and frequency infrastructure, replace MST or Runge–Kutta radial propagation with HeunC plus connection coefficients, and adopt adaptive bi-power mapping quadrature.

This comparison places pybhpt in a broader methodological lineage. It remains a complete frequency-domain Python pipeline for Kerr perturbation theory, metric reconstruction, redshift regularization, and conservative Hamiltonian data, while later HeunC-based work suggests a path toward faster and more stable flux backends for strong-field and highly oscillatory regimes.

Source: https://www.emergentmind.com/topics/pybhpt