pybhpt: Python Library for Kerr Perturbations

• pybhpt is an open-source Python library for first-order black-hole perturbation theory that computes waveforms and invariant quantities for generic eccentric and precessing Kerr orbits.
• It employs a frequency-domain framework to solve the Teukolsky equation using spectral source integration, MST series, and Hertz potential–based metric reconstruction.
• The library links self-force calculations with conservative Hamiltonian data, enabling cross-framework comparisons among effective-one-body and post-Newtonian descriptions.

Searching arXiv for 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 (Nasipak, 10 Jul 2025, Chen et al., 10 May 2026).

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 $\psi_0$0, eccentricities up to $\psi_0$1, and inclinations up to $\psi_0$2 (Nasipak, 10 Jul 2025).

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_0$3 modes. On the conservative side, it computes $\psi_0$4, 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$5 and $\psi_0$6, related to the metric perturbation by linear differential operators. The rescaled scalars satisfy the Teukolsky equations with spin $\psi_0$7, sourced by the point-particle stress-energy tensor.

For bound motion, the frequency spectrum is discrete:

$\psi_0$8

with integers $\psi_0$9, $\psi_4$0, and $\psi_4$1. For equatorial motion, $\psi_4$2; for circular or spherical motion, $\psi_4$3. The radial Teukolsky equation is

$\psi_4$4

with

$\psi_4$5

Retarded boundary conditions are imposed at the horizon and at null infinity using homogeneous solutions $\psi_4$6 and $\psi_4$7 normalized by their asymptotics. The amplitudes are computed through Green’s-function integrals,

$\psi_4$8

with

$\psi_4$9

and

$\psi_0$0

Extended homogeneous solutions are then assembled as

$\psi_0$1

These formulas define the core frequency-domain structure on which the remainder of the library is built (Nasipak, 10 Jul 2025).

The orbital sector is parameterized either by $\psi_0$2 or by orbital elements $\psi_0$3. Introducing Mino time $\psi_0$4 via $\psi_0$5, the motion separates into librations in $\psi_0$6 and $\psi_0$7 and linear drifts in $\psi_0$8 and $\psi_0$9, with fundamental coordinate-time frequencies

$s=+2$0

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 $s=+2$1 or $s=+2$2, leading to ingoing radiation gauge and outgoing radiation gauge. Two build on Aksteiner–Andersson–Bäckdahl identities and use both $s=+2$3 and $s=+2$4; the paper labels these symmetric radiation gauge and antisymmetric radiation gauge.

The reconstructed fields are

$s=+2$5

$s=+2$6

$s=+2$7

$s=+2$8

IRG satisfies $s=+2$9 and $\psi_4$0; ORG satisfies $\psi_4$1 and $\psi_4$2. SRG and ARG are traceless combinations rather than standard gauges. The Aksteiner–Andersson–Bäckdahl identity gives

$\psi_4$3

with $\psi_4$4 in Boyer–Lindquist coordinates. ARG requires time derivatives of Hertz potentials, and static modes with $\psi_4$5 are excluded from ARG in this work (Nasipak, 10 Jul 2025).

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,

$\psi_4$6

For bound geodesics,

$\psi_4$7

The completed solution in each domain and in a target gauge $\psi_4$8 is

$\psi_4$9

This organization is central to pybhpt’s redshift extraction, because continuity of $s=-2$0 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 $s=-2$1 by $s=-2$2. Expanding,

$s=-2$3

The generalized redshift is the infinite coordinate-time average at fixed frequencies $s=-2$4, primary mass $s=-2$5, secondary mass $s=-2$6, and spin $s=-2$7:

$s=-2$8

For bound biperiodic motion, the averaging is written as phase-space integrals in $s=-2$9:

$h_{\mu\nu}$0

$h_{\mu\nu}$1

The resulting quantity is described as a quasi-invariant across physically reasonable gauges (Nasipak, 10 Jul 2025).

Regularization uses a mode-sum compatible with locally Lorenz singular structure:

$h_{\mu\nu}$2

In practice, the leading locally Lorenz parameter $h_{\mu\nu}$3 suffices to render the sum convergent, while pybhpt accelerates convergence by fitting the large-$h_{\mu\nu}$4 tail to determine effective higher-order $h_{\mu\nu}$5. For generic Kerr orbits,

$h_{\mu\nu}$6

with $h_{\mu\nu}$7 and $h_{\mu\nu}$8 given in terms of Kerr metric functions and constants of motion, and $h_{\mu\nu}$9 the complete elliptic integral.

The Hamiltonian connection is

$1\mathrm{PA}$0

with actions $1\mathrm{PA}$1, and

$1\mathrm{PA}$2

Thus pybhpt provides $1\mathrm{PA}$3, 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}$4 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 (Nasipak, 10 Jul 2025).

Module	Role
geo	Geodesic functions, turning points, and fundamental frequencies; spectral evaluation of $1\mathrm{PA}$5 and $1\mathrm{PA}$6
swsh	Spin-weighted spheroidal harmonics via spherical expansions; coupling coefficients $1\mathrm{PA}$7 and eigenvalues $1\mathrm{PA}$8
radial	Homogeneous Teukolsky ODEs with Zenginoğlu stabilization; Teukolsky–Starobinsky mappings between $1\mathrm{PA}$9 and $a/M$0; monodromy-based $a/M$1 optional
teuk	Inhomogeneous amplitudes $a/M$2 via SSI; robust error tracking and thresholds
hertz	Mode functions $a/M$3 using identities; angular and radial derivatives for reconstruction
metric	Analytic reconstruction coefficients for applying GHP operators to Hertz modes; assembling $a/M$4
redshift	Coefficients to form $a/M$5, worldline projection, mode-sum assembly, and regularization/fitting to extract $a/M$6
flux	GW fluxes at infinity and horizon built from $a/M$7 modes

The library’s typical workflow is geodesics $a/M$8 $a/M$9 modes $\psi_0$00 Hertz potentials $\psi_0$01 metric reconstruction $\psi_0$02 $\psi_0$03 $\psi_0$04 $\psi_0$05 $\psi_0$06 $\psi_0$07. In the angular sector, spin-weighted spheroidal harmonics are expanded in spin-weighted spherical harmonics,

$\psi_0$08

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 $\psi_0$09 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 $\psi_0$10-derivatives,

$\psi_0$11

where

$\psi_0$12

Hertz potentials are then built in mode-sum form, reconstruction operators are applied in a spherical-harmonic basis, and the worldline quantity $\psi_0$13 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 $\psi_0$14 singularity near the particle, producing $\psi_0$15 as $\psi_0$16. ARG modes are more singular, with $\psi_0$17–$\psi_0$18 behavior in half-string form, reflecting an $\psi_0$19–$\psi_0$20 divergence, and ARG excludes static modes with $\psi_0$21. SRG inherits the “worst” asymptotic behavior of IRG and ORG at boundaries but remains $\psi_0$22 near the particle (Nasipak, 10 Jul 2025).

The main numerical limitation identified in the paper comes from spectral source integration at very high mode numbers. For large $\psi_0$23, $\psi_0$24, and $\psi_0$25, especially at high eccentricity and small $\psi_0$26, oscillatory cancellation limits accuracy; pybhpt drops modes when the relative error exceeds a threshold. This limits accuracy for $\psi_0$27 near the innermost stable orbit. The paper reports that eccentric equatorial redshift values agree with van de Meent–Shah (2015) at the $\psi_0$28–$\psi_0$29 level across gauges. It also reports new precessing results with spins up to $\psi_0$30 and, for near-extremal Kerr, the first observation of negative $\psi_0$31 near the ISCO, implying a surface in $\psi_0$32 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 $\psi_0$33 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 (Chen et al., 10 May 2026).

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 $\psi_0$34 low-order modes, the HeunC framework achieves relative errors of order $\psi_0$35, with HeunC runtime $\psi_0$36–$\psi_0$37 s versus pybhpt $\psi_0$38–$\psi_0$39 s, corresponding to a $\psi_0$40–$\psi_0$41 speedup. For a highly oscillatory single-mode case, HeunC converges in $\psi_0$42 ms, while pybhpt’s uniform trapezoidal rule over $\psi_0$43 took $\psi_0$44 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.