GammaPBHPlotter: PBH Gamma-Ray Spectra

• GammaPBHPlotter is a Python tool that computes instantaneous per-PBH gamma-ray spectra from Hawking evaporation, including primary, secondary, FSR, and IA emissions.
• It supports both monochromatic and mass-averaged PBH spectra over a wide mass range, integrating models like Gaussian, non-Gaussian, and log-normal mass distributions.
• The code leverages BlackHawk v2.0 for primary and secondary spectral calculations while incorporating analytical methods for final-state radiation and in-flight annihilation to aid observational forecasts.

Searching arXiv for the specified paper and directly related references. GammaPBHPlotter is a public Python code for calculating and plotting the Hawking radiation gamma-ray spectra of primordial black holes (PBHs), with emphasis on the complete gamma-ray output in the mass range relevant for observational searches. It is designed to compute both monochromatic spectra, for a single PBH mass, and mass-averaged spectra, for extended PBH mass distributions, while including four distinct emission components: primary or direct Hawking emission, secondary emission from decay and hadronization of unstable particles, final-state radiation (FSR), and in-flight annihilation (IA) gamma-ray emission. The code is presented as a specialized, open-source tool directly oriented toward phenomenological and observational work on PBH evaporation (Carlini et al., 26 Aug 2025).

1. Scope, regime, and intended use

GammaPBHPlotter targets the calculation of the complete gamma-ray output of PBHs in a mass regime where Hawking radiation may be observable. The paper gives a nominal mass range of $10^{14}$ to $10^{18}\,\mathrm{g}$, while the implementation allows monochromatic masses from $5\times 10^{13}$ to $10^{19}\,\mathrm{g}$. PBHs lighter than $\sim 10^{14}\,\mathrm{g}$ would already have evaporated, whereas more massive PBHs emit at lower temperatures and thus at lower characteristic photon energies.

The code addresses several distinct but closely related tasks. It computes monochromatic spectra for a fixed PBH mass, computes mass-averaged spectra by integrating over PBH mass functions, and resolves the gamma-ray signal into physically distinct components. Its stated use cases include forecasting PBH signals for gamma-ray telescopes such as e-ASTROGAM and AMEGO-X, building PBH evaporation models for dark matter or cosmology studies, comparing different PBH mass functions and their spectral signatures, and producing spectral inputs for abundance constraints based on diffuse gamma-ray backgrounds or point-source limits.

This focus places the code within PBH phenomenology rather than within detector simulation or cosmological transport modeling. The object being computed is the instantaneous per-PBH gamma-ray spectrum; conversion into observational fluxes, line-of-sight integrals, or cosmological backgrounds is left to downstream analysis.

2. Physical framework and relation to Hawking evaporation

GammaPBHPlotter does not rederive the Hawking evaporation formalism internally. Instead, it delegates the primary and secondary emission calculations to BlackHawk v2.0, described as a well-established PBH evaporation code implementing Hawking radiation with Standard Model particle content and numerical greybody factors (Carlini et al., 26 Aug 2025).

For a Schwarzschild black hole of mass $M$, the black-hole temperature is given in the paper’s description as

$$T_\mathrm{BH} = \frac{\hbar c^3}{8\pi G M k_\mathrm{B}},$$

so lower-mass PBHs are hotter and radiate more energetic particles. In the mass range near $\sim 10^{14}\,\mathrm{g}$, the characteristic temperatures reach the MeV-GeV regime.

The underlying differential Hawking emission rate is summarized as

$$\frac{\mathrm{d}^2 N_i}{\mathrm{d}E\,\mathrm{d}t}(E,M) = \frac{\Gamma_i(E,M)}{2\pi\hbar}\, \frac{1}{\exp\left(E/k_\mathrm{B} T_\mathrm{BH}\right)-(-1)^{2s_i}},$$

where $\Gamma_i(E,M)$ denotes the greybody factor for species $10^{18}\,\mathrm{g}$0 with spin $10^{18}\,\mathrm{g}$1. In this formulation, the greybody factors encode the frequency-dependent transmission probabilities through the black-hole potential barrier. BlackHawk provides numerical tabulations of these factors for the relevant Standard Model species.

The description also recalls the mass-loss relation

$10^{18}\,\mathrm{g}$2

BlackHawk can treat temporal evolution, but GammaPBHPlotter is explicitly focused on instantaneous gamma-ray spectra for specified PBH masses rather than on dynamical mass evolution over cosmic time.

The particle content adopted is Standard Model only. Direct emission includes photons, electrons and positrons, neutrinos, quarks, gluons, and, once temperature thresholds are crossed, heavier particles such as muons, hadrons, and $10^{18}\,\mathrm{g}$3 bosons. Any beyond-Standard-Model degrees of freedom are not part of the paper’s description.

3. Gamma-ray components and spectral decomposition

The total gamma-ray spectrum is decomposed into four components. Primary and secondary spectra are imported from BlackHawk; FSR and IA are added explicitly within GammaPBHPlotter.

Component	Origin	Modeling source
Primary/direct	Photons emitted at the horizon	BlackHawk v2.0
Secondary	Decay and hadronization products of unstable Hawking-emitted particles	BlackHawk v2.0 with PYTHIA 8.2
Final-state radiation	Internal bremsstrahlung from relativistic $10^{18}\,\mathrm{g}$4	Analytical QED expression
In-flight annihilation	Relativistic positron annihilation with ambient electrons before thermalization	Formula from Keith et al. (2022)

The primary or direct component consists of photons emitted directly through Hawking radiation. In practice, GammaPBHPlotter ingests the photon spectrum computed by BlackHawk, where the photon greybody factors appropriate to spin-1 fields in the Schwarzschild geometry are already included. The direct component tends to dominate at higher photon energies, near the peak of the Hawking spectrum, for sufficiently hot PBHs.

The secondary component originates from Hawking emission of quarks and gluons, subsequent hadronization into hadronic jets, and decays of unstable hadrons, especially $10^{18}\,\mathrm{g}$5, together with more extended cascade channels. BlackHawk performs this modeling using PYTHIA 8.2 for hadronization and decay. Conceptually, the resulting spectrum is a convolution of primary Hawking parton spectra with fragmentation and decay kernels. In the description, secondary emission is characterized as typically dominant at lower photon energies, producing a broad spectrum extending below the Hawking peak.

FSR is treated as internal bremsstrahlung from relativistic charged leptons, especially positrons, emitted in Hawking evaporation. The code evaluates the FSR gamma-ray spectrum through an analytical QED formula based on the Hawking-emitted $10^{18}\,\mathrm{g}$6 spectrum. The expression is written as

$10^{18}\,\mathrm{g}$7

with $10^{18}\,\mathrm{g}$8, $10^{18}\,\mathrm{g}$9 the kinetic energy of the positron, $5\times 10^{13}$0 the photon energy, and $5\times 10^{13}$1. In this construction, the FSR contribution fills in additional gamma-ray flux at energies below the characteristic lepton energies.

IA is distinct from the traditional $5\times 10^{13}$2 keV line from thermalized positrons. It denotes gamma rays produced when relativistic positrons annihilate with ambient electrons before they have cooled to nonrelativistic energies. GammaPBHPlotter computes this component using a formula from Keith et al. (2022): $5\times 10^{13}$3 where the kernel is

$5\times 10^{13}$4

The description specifies $5\times 10^{13}$5 as the assumed interstellar hydrogen and electron density, Bethe-Bloch energy losses, and an annihilation probability

$5\times 10^{13}$6

IA is therefore environment-dependent by construction, and the default result corresponds to a typical interstellar medium assumption.

4. Monochromatic and mass-averaged PBH spectra

GammaPBHPlotter supports both single-mass PBH populations and extended PBH mass functions. The monochromatic case assumes all PBHs have the same mass $5\times 10^{13}$7. In practice, the user chooses a mass between $5\times 10^{13}$8 and $5\times 10^{13}$9, BlackHawk computes the primary and secondary spectra at that mass, and GammaPBHPlotter adds FSR and IA to produce $10^{19}\,\mathrm{g}$0 per PBH.

The paper’s example for a monochromatic spectrum is a PBH of mass $10^{19}\,\mathrm{g}$1, for which the total spectrum and the individual components are displayed. This example is used to show the relative contributions of direct emission, secondary emission, FSR, and IA.

Mass-averaged spectra are obtained by integrating the instantaneous emission over a PBH mass distribution: $10^{19}\,\mathrm{g}$2 where $10^{19}\,\mathrm{g}$3 is the PBH mass function with the appropriate normalization. The code supports four classes of mass distributions:

1. Monochromatic distribution: a delta function at mass $10^{19}\,\mathrm{g}$4.
2. Gaussian mass distribution from Gaussian density perturbations: derived from the model of Biagetti et al. (2021), Gaussian around a mean mass, with $10^{19}\,\mathrm{g}$5 identified as the standard deviation of the initial density perturbations.
3. Non-Gaussian mass distribution (Biagetti et al.): a PBH mass function derived from non-Gaussian perturbations, giving richer tails and skewness than the Gaussian case.
4. Log-normal mass distribution: a log-normal in PBH mass, controlled by a central mass and a width.

The comparison shown in the second figure contrasts a single $10^{19}\,\mathrm{g}$6 PBH with Gaussian mass-averaged spectra having the same mean mass but different $10^{19}\,\mathrm{g}$7. The stated effect is spectral broadening and reshaping induced by the extended mass function.

5. Software architecture, configuration, and outputs

GammaPBHPlotter is written in Python 3.9 and uses numpy, scipy, matplotlib, tqdm, and colorama (Carlini et al., 26 Aug 2025). It is described as running on Windows, Linux, and macOS. The software also “automatically checks and downloads all missing modules,” and a user manual is provided in the Zenodo package.

The main configurable inputs are the PBH mass or mass-distribution parameters, the selected emission components, the spectral settings, and the environmental parameters entering IA. The description explicitly lists monochromatic mass values and the parameters of Gaussian, non-Gaussian, or log-normal PBH mass functions. It also states that the code includes the four emission components—primary/direct Hawking photons, secondary emission from decay and hadronization, FSR from $10^{19}\,\mathrm{g}$8, and IA—implying that users can choose among them in the code configuration. For IA, the hydrogen and electron density defaults to $10^{19}\,\mathrm{g}$9.

There is no explicit mention of cosmological distance or redshift inputs. The code focuses on per-PBH instantaneous spectra in physical units, and the conversion to observable fluxes is left to the user. This separation makes the code a spectral backend rather than a full end-to-end observational analysis package.

Its principal outputs are differential gamma-ray spectra,

$\sim 10^{14}\,\mathrm{g}$0

resolved into primary, secondary, FSR, IA, and their sum. For extended PBH populations, the same quantities are produced after averaging over the chosen mass distribution. Visualization is carried out with matplotlib in log-log format. The paper identifies two representative outputs: Monochromatic.png, showing the total spectrum and its components for a $\sim 10^{14}\,\mathrm{g}$1 PBH, and Spectrum_comparison.png, comparing a monochromatic spectrum with mass-averaged spectra for Gaussian perturbation models of identical mean mass and varying $\sim 10^{14}\,\mathrm{g}$2.

The Zenodo record specified for the code is

$\sim 10^{14}\,\mathrm{g}$3

listed as Carlini & Cholis, “GammaPBHPlotter: code and user manual”, 2025.

6. Validation status, assumptions, limitations, and observational role

The paper is explicitly short and tool-oriented, and it does not present an extensive standalone validation campaign. Its validation strategy is largely inherited from its components. Primary and secondary spectra come from BlackHawk v2.0, which is described as having been thoroughly compared with previous PBH evaporation calculations. The hadronization and decay treatment relies on PYTHIA 8.2, identified as a standard Monte Carlo framework for such processes. The added FSR contribution uses an analytical QED bremsstrahlung expression, and IA follows the methodology of Keith et al. (2022), developed in the context of PBH sensitivity forecasts (Carlini et al., 26 Aug 2025).

Within that framework, the main improvement over a workflow built from “raw” BlackHawk plus separate scripts is the packaging of the complete gamma-ray spectrum, including FSR and IA, into a single accessible Python tool. A second stated improvement is the ready-made support for several PBH mass distributions motivated by formation models, including those associated with Biagetti et al. A third is the integrated plotting and module-management functionality.

The assumptions and caveats are correspondingly clear. The particle content is Standard Model only. The black holes are treated as non-rotating Schwarzschild PBHs; rotation and charge are not discussed. IA is computed in a homogeneous interstellar medium with $\sim 10^{14}\,\mathrm{g}$4 and Bethe-Bloch energy losses, so the IA spectrum is environment-specific and would change in denser or more rarefied media. The code is restricted to PBH masses between $\sim 10^{14}\,\mathrm{g}$5 and $\sim 10^{14}\,\mathrm{g}$6. The final explosive phase below $\sim 10^{14}\,\mathrm{g}$7 is not emphasized, and at the high-mass end the gamma-ray emission becomes weak because of low black-hole temperatures. The code does not model PBH clustering, halo profiles, or cosmological evolution, and it does not fold the emission with instrument response functions or background models.

The observational connection is nevertheless direct. For a single PBH at distance $\sim 10^{14}\,\mathrm{g}$8, the spectrum may be converted to an observable flux using

$\sim 10^{14}\,\mathrm{g}$9

For a PBH population with number density $M$0 and mass function $M$1, the description gives

$M$2

Such flux constructions allow the code’s outputs to be compared with gamma-ray data or projected sensitivities from instruments including Fermi-LAT, INTEGRAL, COMPTEL, e-ASTROGAM, AMEGO-X, and H.E.S.S. In this sense, GammaPBHPlotter functions as a spectral engine for PBH phenomenology: it computes the instantaneous component-resolved gamma-ray signal, while abundance limits, line-of-sight modeling, and detector-level inference are performed externally.