GrayHawk: Gray-Body Factor Solver

• GrayHawk is a public Mathematica-based numerical code that calculates gray-body factors for massless fields in static, spherically symmetric spacetimes.
• It solves a one-dimensional scattering problem using a shooting method and asymptotic fitting to deduce transmission coefficients filtered by curvature-induced potentials.
• The tool supports various black-hole metrics and traversable wormhole geometries, serving as a modular backend for Hawking radiation spectral analysis.

GrayHawk is a public Mathematica-based numerical code for computing gray-body factors, or transmission coefficients, for massless fields propagating on static, spherically symmetric spacetimes. In its original release, it was designed for four-dimensional, asymptotically flat black holes and for massless scalar, Dirac, electromagnetic, and gravitational fields with spins $s=0,\tfrac12,1,2$; it was also intended as a general, easy-to-extend engine for Hawking-radiation calculations and as a backend for BlackHawk v3.0. GrayHawk v2 enlarges that scope in two disjoined directions by enabling a fully numerical computation of the tortoise-coordinate integral and by extending the scattering problem to wormhole solutions (Calzá, 6 Feb 2025, Calzá, 4 Jun 2026, Arbey et al., 4 Jun 2026).

1. Scope, physical role, and observables

GrayHawk addresses the scattering problem that underlies Hawking emission. For a given angular mode, the gray-body factor $\Gamma_{s,l}(\omega)$ encodes the frequency-dependent transmission probability through the curvature-induced potential barrier outside the horizon. Hawking’s derivation yields a thermal spectrum at the horizon, but the spectrum observed at infinity is filtered by that barrier. In the notation used in the original presentation,

$$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$

with $+$ for fermions, $-$ for bosons, $\omega=E_i$, and $n_i$ the internal degrees of freedom (Calzá, 6 Feb 2025).

The first public version is restricted to static, spherically symmetric, asymptotically flat black-hole metrics in $4$D and to massless fields. Within that domain, the intended use cases include precision modeling of Hawking evaporation, especially for primordial black holes, comparison of general-relativistic black holes with modified or quantum-corrected metrics, and the derivation of observational constraints from emission spectra. The v2 paper generalizes the description from black holes alone to curved, static, spherically symmetric spacetimes, and explicitly includes traversable wormholes and black-bounce metrics in scattering and echo studies (Calzá, 4 Jun 2026).

This positioning gives GrayHawk a specific role in the computational ecosystem of black-hole phenomenology. It is not primarily a full evaporation code; rather, it is a solver for the transmission problem whose outputs can be inserted into Hawking-spectrum calculations. That division of labor becomes explicit in its later use inside BlackHawk v3.0, where GrayHawk-generated tables supply the greybody factors for newly implemented spherically symmetric metrics (Arbey et al., 4 Jun 2026).

2. Formalism and wave equation

GrayHawk starts from the general static, spherically symmetric line element

$$ds^2 = -G(r)\,dt^2 + \frac{dr^2}{F(r)} + H(r)\,d\Omega^2,$$

with asymptotic-flatness conditions

$$F(r)\xrightarrow[r\to\infty]{}1,\qquad G(r)\xrightarrow[r\to\infty]{}1,\qquad H(r)\sim r^2.$$

A special subclass is the time-radius symmetric, or $\Gamma_{s,l}(\omega)$0-symmetric, case with $\Gamma_{s,l}(\omega)$1 and $\Gamma_{s,l}(\omega)$2. Schwarzschild, Reissner–Nordström, Bardeen, and Hayward belong to that subclass, whereas some black-bounce metrics do not (Calzá, 6 Feb 2025).

The analytical framework is based on the Newman–Penrose formalism and a Teukolsky-type treatment of massless fields. In the original formulation, a single master equation is derived for $\Gamma_{s,l}(\omega)$3. Separation of variables is performed with

$\Gamma_{s,l}(\omega)$4

where $\Gamma_{s,l}(\omega)$5 are spin-weighted spherical harmonics and the separation constant is $\Gamma_{s,l}(\omega)$6. The v2 paper presents the same general structure and states that the formalism focuses on massless fields with spins $\Gamma_{s,l}(\omega)$7, while also noting in its limitations that spin-$\Gamma_{s,l}(\omega)$8 is not yet implemented in the public code (Calzá, 4 Jun 2026).

The radial problem is recast through the generalized tortoise coordinate

$\Gamma_{s,l}(\omega)$9

into a Schrödinger-like equation

$$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$0

The effective potential $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$1 depends on the spin, the angular quantum number, and the metric functions $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$2, $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$3, and $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$4. For asymptotically flat black holes with a regular horizon, $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$5 as $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$6, so the asymptotic solutions are plane waves. GrayHawk imposes a purely ingoing, unit-normalized wave at the horizon and extracts the asymptotic coefficients at spatial infinity, leading to the identification

$$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$7

This formulation makes the code a direct solver of a one-dimensional scattering problem rather than a purely thermodynamic Hawking-spectrum package (Calzá, 6 Feb 2025).

3. Numerical strategy and code architecture

The central numerical task is to solve the Schrödinger-like radial equation for each frequency and angular mode, under the purely ingoing horizon boundary condition, and to infer the transmission coefficient from the far-field asymptotics. The algorithm is described as a shooting method with asymptotic fitting: impose the near-horizon ingoing wave, integrate outward, fit the large-$$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$8 solution to a plane-wave superposition, and compute $$\frac{d^2 N_i}{dt\,dE_i} = \frac{1}{2\pi} \sum_{l,m} \frac{n_i\,\Gamma^{\,s}_{l,m}(\omega)}{e^{\omega/T} \pm 1},$$9 (Calzá, 6 Feb 2025).

The original public release is organized into two Mathematica notebooks. GrayHawk.nb is the main engine that computes gray-body factors over a table of energies, and CalibratorGH.nb is a helper notebook used to tune the sampling and check numerical stability. The main notebook is structured around parameter assignment, metric definition, tortoise-coordinate construction and inversion, geometric-potential construction, and repeated ODE solves over a user-defined energy list. Default settings correspond to Schwarzschild, $+$0, and $+$1 linearly spaced energies between $+$2 and $+$3 in units $+$4 (Calzá, 6 Feb 2025).

Accurate construction of the tortoise map is a critical step. In the first version, GrayHawk attempts to evaluate

$+$5

with Mathematica’s Integrate when possible, then samples $+$6, numerically inverts the relation to obtain $+$7, and interpolates that inverse map for use in the potential. The default sampling splits the radial domain into a near-horizon logarithmic region from $+$8 to $+$9 with $-$0 points and a far linear region from $-$1 to $-$2 with $-$3 points. The far-zone fit is performed on a subinterval of length x_far, and the calibrator is explicitly intended to detect two common pathologies: insufficiently oscillatory low-energy solutions in the far region and high-energy fits with $-$4 caused by a starting point not close enough to the horizon. The documentation also notes the need to respect the Nyquist–Shannon sampling theorem in the far zone (Calzá, 6 Feb 2025).

GrayHawk v2 preserves the two-notebook architecture as GrayHawkv2.nb and CalibratorGHv2.nb, but adds a flag-controlled fully numerical tortoise construction. The new option is introduced precisely because the analytic integral can fail or be impractical for metrics such as ZLMY and CGSV. Near the horizon, the integrand divergence is handled with a first-order Taylor expansion that reproduces the correct logarithmic behavior, after which the code builds a numerical table $-$5, inverts it, and interpolates $-$6 on a finite interval chosen so that the potential is practically zero at the boundaries (Calzá, 4 Jun 2026).

4. Implemented geometries and extensibility

The original code is preloaded with seven asymptotically flat, static, spherically symmetric black-hole metrics. All are written in the form $-$7, with mass parameter $-$8 set to $-$9 in the code and, in several cases, an additional deformation or regularization parameter. The precompiled set is as follows (Calzá, 6 Feb 2025).

Metric	Parameters	Notes
Schwarzschild	$\omega=E_i$0	Default metric; $\omega=E_i$1
Reissner–Nordström	$\omega=E_i$2	Black-hole configurations assumed
Hayward	$\omega=E_i$3	Regular; reduces to Schwarzschild as $\omega=E_i$4
Bardeen	$\omega=E_i$5	Regular; horizon requires $\omega=E_i$6
Simpson–Visser	$\omega=E_i$7	Black-hole regime for $\omega=E_i$8
Peltola–Kunstatter	$\omega=E_i$9	LQG-inspired regular metric
D’Ambrosio–Rovelli	$n_i$0	Regular black-hole to white-hole geometry

These choices show the intended balance between standard general-relativistic metrics, regular black holes, black-bounce geometries, and loop-quantum-gravity-inspired constructions. In the original description, several of the regular metrics are characterized by de Sitter cores or bounce interiors, and the Simpson–Visser family is explicitly noted to interpolate between black holes and traversable wormholes, although GrayHawk v1 is meant to be used in the black-hole regime (Calzá, 6 Feb 2025).

Extensibility is deliberately minimal. To add a new spherically symmetric black hole, the user provides $n_i$1, $n_i$2, and $n_i$3 for a metric of the required asymptotic form, after which the code finds the horizon as the largest real root of $n_i$4, constructs and inverts the tortoise coordinate, and builds the potential from the general formulas. The original paper gives Frolov’s regular black hole as a worked example of this modification path (Calzá, 6 Feb 2025).

GrayHawk v2 broadens the preloaded pool in two ways. First, it adds black-hole metrics for which analytic tortoise coordinates are unavailable or inconvenient, including Culetu–Ghosh–Simpson–Visser, Dymnikova, and ZLMY. Second, it introduces wormhole geometries such as Ellis–Bronnikov, Dadhich–Kar–Mukherji–Visser, and thin-shell wormholes built by a general cut-and-paste prescription. For smooth metrics the code can determine whether the chosen parameters describe a black hole or a wormhole; for thin-shell mode, enabled through SmoothOrNot=0, it asks the user to verify that the chosen cut radius indeed produces a thin-shell wormhole (Calzá, 4 Jun 2026).

5. Validation, performance, and relation to BlackHawk

The original validation program compares GrayHawk against highly accurate photon spectra previously computed with a Frobenius-based method for Hayward, Bardeen, Simpson–Visser, Peltola–Kunstatter, and D’Ambrosio–Rovelli black holes. In those tests, with photon spectra summed up to $n_i$5, the curves are reported to be visually indistinguishable. Quantitatively, the peak-energy discrepancy is $n_i$6, the peak-intensity discrepancy is $n_i$7, and over more than two decades around the peak the residuals remain below $n_i$8. The same source also gives indicative runtimes: default Schwarzschild settings produce output in $n_i$9 seconds on a typical laptop, while more demanding cases such as spin-$4$0, $4$1 around a D’Ambrosio–Rovelli black hole can take $4$2 seconds (Calzá, 6 Feb 2025).

GrayHawk’s role in BlackHawk v3.0 is more specific. There, GrayHawk is described as the companion code that computes the greybody factors required for the new spherically symmetric metrics implemented in BlackHawk, while BlackHawk itself remains the user-facing code that tabulates those outputs and builds primary Hawking spectra. The resulting tables cover spins $4$3, angular momentum up to $4$4, $4$5 values of the normalized regularization parameter $4$6, and $4$7 energy points in the range $4$8. Additional GrayHawk-derived runs at $4$9, $$ds^2 = -G(r)\,dt^2 + \frac{dr^2}{F(r)} + H(r)\,d\Omega^2,$$0, $$ds^2 = -G(r)\,dt^2 + \frac{dr^2}{F(r)} + H(r)\,d\Omega^2,$$1, and $$ds^2 = -G(r)\,dt^2 + \frac{dr^2}{F(r)} + H(r)\,d\Omega^2,$$2 are used to construct asymptotic fits outside the main tabulated domain (Arbey et al., 4 Jun 2026).

In that integration pipeline, BlackHawk does not solve the scattering problem at runtime. Instead, GrayHawk-generated tables are stored under src/tables/gamma_tables/..., and Appendix material in the BlackHawk v3.0 paper documents the GrayHawk-based scripts used to generate or extend those tables for Bardeen, Hayward, Simpson–Visser, Peltola–Kunstatter, D’Ambrosio–Rovelli, and BCL metrics. This suggests a stable division of responsibilities: GrayHawk acts as the transmissivity engine, while BlackHawk handles interpolation, extrapolation, and assembly of Hawking emission observables (Arbey et al., 4 Jun 2026).

GrayHawk v2 is validated in a different regime. The paper reports agreement with published transmission coefficients for thin-shell Schwarzschild wormholes and for ZLMY quantum-gravity black-hole and wormhole solutions, with residuals of order $$ds^2 = -G(r)\,dt^2 + \frac{dr^2}{F(r)} + H(r)\,d\Omega^2,$$3 except very close to sharp narrow peaks, where small energy shifts inflate relative discrepancies. Those tests are presented as validation of both the numerical tortoise implementation and the wormhole scattering module (Calzá, 4 Jun 2026).

6. Limitations, versioning, and prospective extensions

The core limitations are explicit. The original GrayHawk is restricted to static, spherically symmetric, asymptotically flat black holes in four dimensions and to massless fields. It does not include rotating geometries such as Kerr or Kerr–Newman, does not handle spacetimes with a cosmological constant, and does not implement massive fields or nonminimal couplings. It also notes practical difficulties for metrics whose tortoise coordinates cannot be integrated analytically in the default implementation (Calzá, 6 Feb 2025).

GrayHawk v2 removes one of those bottlenecks by enabling a fully numerical tortoise-coordinate computation and extends the scattering problem to traversable wormholes, including smooth and thin-shell cases. Even after that extension, however, the scope remains four-dimensional, static, spherically symmetric, and asymptotically flat. The v2 discussion identifies several future directions: massive fields, rotating backgrounds inside GrayHawk itself, dedicated modules for quasinormal modes and resonance searches, improved algorithms for resonances and automated parameter scans, interfaces with waveform-analysis pipelines, and additional exotic compact-object geometries. It also states that adding spin-$$ds^2 = -G(r)\,dt^2 + \frac{dr^2}{F(r)} + H(r)\,d\Omega^2,$$4 fields remains a realistic and natural possibility for later releases (Calzá, 4 Jun 2026).

Within the BlackHawk v3.0 context, an additional limitation concerns dynamics rather than scattering. The code computes instantaneous spectra for a chosen mass and regularization parameter, but does not attempt to model the time evolution of regular black holes unless a consistent dynamical framework is available. A plausible implication is that GrayHawk’s present role is deliberately local in parameter space: it provides high-quality greybody factors for fixed backgrounds, while questions of self-consistent metric evolution remain external to the solver (Arbey et al., 4 Jun 2026).

Across these versions, the defining characteristic is the preservation of a modular structure. In the original release that modularity supports extension to new black-hole metrics; in v2 it is retained while adding numerical tortoise routines and wormhole support; and in the BlackHawk v3.0 pipeline it allows GrayHawk outputs to be reused as precomputed greybody-factor tables. This suggests that GrayHawk is best understood as a reusable scattering backend for spherically symmetric geometries rather than as a single-purpose notebook tied to a fixed menu of metrics (Calzá, 6 Feb 2025, Calzá, 4 Jun 2026).