# BlackHawk: A public code for calculating the Hawking evaporation spectra of any black hole distribution

Alexandre Arbey, Jérémy Auffinger

► To cite this version:

Alexandre Arbey, Jérémy Auffinger. BlackHawk: A public code for calculating the Hawking evaporation spectra of any black hole distribution. European Physical Journal C: Particles and Fields, 2019, 79 (8), pp.693. 10.1140/epjc/s10052-019-7161-1 . hal-02144274

HAL Id: hal-02144274

https://hal.science/hal-02144274v1

Submitted on 28 May 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# BlackHawk: a public code for calculating the Hawking evaporation spectra of any black hole distribution

Alexandre Arbey<sup>1,2,5,a</sup>, Jérémy Auffinger<sup>1,3,4,b</sup>

<sup>1</sup> Univ Lyon, Univ Lyon 1, CNRS/IN2P3, Institut de Physique Nucléaire de Lyon, UMR5822, 69622 Villeurbanne, France

<sup>2</sup> Theoretical Physics Department, CERN, CH-1211 Geneva 23, Switzerland

<sup>3</sup> Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

<sup>4</sup> Département de Physique, École Normale Supérieure de Lyon, 69342 Lyon, France

<sup>5</sup> Also Institut Universitaire de France, 103 boulevard Saint-Michel, 75005 Paris, France

Received: 22 May 2019 / Accepted: 18 July 2019 / Published online: 19 August 2019

© The Author(s) 2019

**Abstract** We describe BlackHawk, a public C program for calculating the Hawking evaporation spectra of any black hole distribution. This program enables the users to compute the primary and secondary spectra of stable or long-lived particles generated by Hawking radiation of the distribution of black holes, and to study their evolution in time. The physics of Hawking radiation is presented, and the capabilities, features and usage of BlackHawk are described here under the form of a manual. The BlackHawk code can be downloaded from https://blackhawk.hepforge.org.

## Contents

1 Introduction .................................................... 2

2 Physics of Hawking radiation ................................. 2

2.1 Testing Black Holes distributions .................... 2

2.1.1 Peak theory distribution .......................... 2

2.1.2 Log-normal distribution .......................... 3

2.1.3 Power-law distribution .......................... 3

2.1.4 Critical collapse distribution ...................... 3

2.1.5 Dirac distribution .................................. 3

2.2 Hawking radiation ........................................ 3

2.2.1 Schwarzschild Black Holes .......................... 3

2.2.2 Kerr Black Holes .................................... 4

2.2.3 Exotic Black Holes .................................. 4

2.3 Black Hole evolution .................................... 4

2.3.1 Schwarzschild Black Holes .......................... 4

2.3.2 Kerr Black Holes .................................... 5

2.3.3 Exotic Black Holes .................................. 5

2.4 Hadronization ............................................ 5

3 Content and compilation .................................... 5

3.1 Main directory ............................................ 5

3.2 src/ sub-folder ............................................ 5

3.3 results/ sub-folder .................................... 6

3.4 manual/ sub-folder .................................... 6

3.5 scripts/ sub-folder .................................... 6

3.6 Compilation ................................................ 6

4 Input parameters ............................................ 6

4.1 Parameter structure .................................... 6

4.2 General parameters .................................... 7

4.3 Black Hole spectrum parameters ...................... 7

4.4 Black Holes evolution parameters ...................... 7

4.5 Primary spectrum parameters .......................... 8

4.6 Hadronization parameters .............................. 8

5 Routines .................................................... 8

5.1 General routines ............................................ 8

5.2 Black Hole spectrum routines .......................... 9

5.3 Black Holes evolution routines .......................... 9

5.4 Primary spectra routines ................................ 10

5.5 Secondary spectra routines ............................ 11

5.6 Auxiliary routines ........................................ 12

6 Programs ...................................................... 12

6.1 Common features ........................................ 12

6.2 BlackHawk_tot: Time-dependent Hawking spectra .................................................... 13

6.3 BlackHawk_inst: Instantaneous Hawking spectra .................................................... 13

7 Output files .................................................... 14

7.1 BlackHawk_tot ............................................ 14

7.2 BlackHawk_inst ............................................ 15

8 Memory use .................................................... 15

8.1 RAM used .................................................... 15

8.2 Static disk memory used ................................ 16

9 Other applications ............................................ 16

9.1 Computing new numerical tables ...................... 16

<footnote><sup>a</sup> e-mail: alexandre.arbey@ens-lyon.fr</footnote>

<footnote><sup>b</sup> e-mail: j.auffinger@ipnl.in2p3.fr</footnote>
<table><tbody><tr><td>9.2</td><td>Using another Black Hole mass function</td><td>17</td></tr><tr><td>9.3</td><td>Adding primary particles</td><td>17</td></tr><tr><td>9.4</td><td>Adding secondary particles</td><td>17</td></tr><tr><td>9.5</td><td>Other types of Black Holes</td><td>17</td></tr><tr><td>10</td><td>Conclusion</td><td>17</td></tr><tr><td></td><td>Appendix A: Units</td><td>18</td></tr><tr><td></td><td>Appendix A.1: Energy</td><td>18</td></tr><tr><td></td><td>Appendix A.2: Mass</td><td>18</td></tr><tr><td></td><td>Appendix A.3: Time</td><td>18</td></tr><tr><td></td><td>Appendix A.4: Distance</td><td>18</td></tr><tr><td></td><td>Appendix A.5: Temperature</td><td>18</td></tr><tr><td></td><td>Appendix B: Computation of the tables</td><td>18</td></tr><tr><td></td><td>Appendix B.1: Greybody factors</td><td>18</td></tr><tr><td></td><td>Appendix B.2: Evolution tables</td><td>20</td></tr><tr><td></td><td>Appendix B.3: Hadronization</td><td>20</td></tr><tr><td></td><td>Appendix C: Particle information</td><td>20</td></tr><tr><td></td><td>Appendix D: Results</td><td>21</td></tr><tr><td></td><td>Appendix D.1: Parameters</td><td>21</td></tr><tr><td></td><td>Appendix D.2: BlackHawk_tot</td><td>22</td></tr><tr><td></td><td>Appendix D.3: BlackHawk_inst</td><td>24</td></tr><tr><td></td><td>References</td><td>25</td></tr></tbody></table>

# 1 Introduction

Black Holes (BHs) are fundamental objects which are of utmost importance for the understanding of gravitation. With the detection of gravitational waves from mergers of binary BHs [1–3], direct observation of the Milky Way supermassive central BH [4], and the cosmological and gravitational questions related to primordial BHs (PBHs, see for example [5–8]), these compact objects are currently under intense scrutiny. It is therefore important to find methods to characterize their properties, and we present here a program for studying multi-messenger probes of BHs.

Other codes, such as BlackMax [9] and Charybdis [10], have already been released in order to compute the Hawking radiation (HR) of BHs, which however focus on higher-dimensional models of general relativity where the Planck mass is decreased and allow the users to make predictions for generation and evaporation of micro black holes at high-energy colliders.

We present here BlackHawk, which is the first public code for the computation of the time-dependent HR into stable or long-lived particles of 4-dimensional Schwarzschild and Kerr BHs distributed in mass.

This document constitutes the manual of BlackHawk v1.0 and is organized as follows: Sect. 2 is a brief overview of BHs and HR physics, Sect. 3 presents the structure and file content of the code, and the compilation and run instructions, Sect. 4 describes the input parameters needed to run BlackHawk, Sect. 5 gives a detailed description of all the routines written in the code, Sect. 6 follows the normal execution of BlackHawk programs and gives examples of screen output, Sect. 7 presents the format of the data files generated by a run along with examples, Sect. 8 gives an estimation of the memory usage and Sect. 9 provides instructions for the users on how to modify the code.

# 2 Physics of Hawking radiation

In this section we give a short overview of the main physical aspects of HR. This concerns BHs of primordial origin (PBHs), as well as any other BHs.

In the following, all formulas are in natural units where $\hbar = c = k_B = G = 1$, unless stated otherwise.

## 2.1 Testing Black Holes distributions

BlackHawk has been designed to provide tests of compatibility between observations and BH distributions at different main steps of the history of the Universe. For this purpose, it computes the HR of a distribution of BHs, and its evolution in time. The obtained spectra can then be used to check whether the amount of produced particles has an effect on observable cosmological quantities.

The distribution of BHs as a function of their mass is completely model-dependent and recent studies have proven some previously set constraints to be irrelevant [11,12]. BlackHawk can in principle work with any distribution of BHs. Several BH mass functions are already built-in and depend on the details of the BH formation mechanisms. For these built-in mass functions, all BHs are considered to have the same spin.

### 2.1.1 Peak theory distribution

The peak theory distribution is derived from the scale-invariant model, assuming that the power spectrum of the primordial density fluctuations is a power-law (see e.g. [13,14])

$$P(k) = R_c \left(\frac{k}{k_0}\right)^{n-1}, \qquad (1)$$

where $n \approx 1.3$ and $R_c$ is measured using the Cosmic Microwave Background (CMB) to be $R_c = (24.0 \pm 1.2) \times 10^{-10}$ at the scale $k_0 = 0.002 \text{ Mpc}^{-1}$. The comoving number density of PBHs resulting from this power spectrum is obtained in [13] through peak-theory

$$dn \approx \frac{1}{4\pi^2 M} \left(\frac{X(n-1)}{6M}\right)^{3/2} \frac{(n-1)}{2} \nu^4 e^{-\nu^2/2} dM, \quad (2)$$

where

$$\nu(M) = \left( \frac{2(k_0^2 M/X)^{(n-1)/2}}{R_c \Gamma((n-1)/2)} \right)^{1/2} \zeta_{\text{th}}, \quad (3)$$
and

$$X = \frac{4\pi}{3} \left(\frac{8\pi G}{3}\right)^{-1} \left(\frac{H_0^2 \Omega_m}{1 + z_{\text{eq}}} \left(\frac{g_{*\text{eq}}}{g_*}\right)^{1/3}\right)^{1/2}, \quad (4)$$

in which:

- $H_0 = 67.8 \text{ km s}^{-1} \text{ Mpc}^{-1}$ is the current Hubble parameter [15].

- $\Omega_m = 0.308$ is the matter mass fraction in the Universe [15].

- $z_{\text{eq}} = 3200$ is the radiation-matter equality redshift [13].

- $g_{*\text{eq}} = 3.36$ is the number of relativistic energy degrees of freedom (dof) at radiation-matter equality [13].

- $g_* = 106.75$ is the number of relativistic energy dof at the time of PBH formation (here the end of the inflation) [15].

- $\zeta_{\text{th}} = 0.7$ parametrizes the direct collapse of a primordial density fluctuation into a PBH [13].

### 2.1.2 Log-normal distribution

The log-normal distribution [11] is considered to be the general mass function originating from a peak in the power spectrum of primordial fluctuations. It is parametrized through

$$dn = \frac{A}{\sqrt{2\pi}\sigma M^2} \exp\left(-\frac{\ln(M/M_c)^2}{2\sigma^2}\right) dM, \quad (5)$$

where A is the amplitude, $M_c$ is the position of the peak and $\sigma$ is its width. Note that this is a log-normal distribution for the comoving density $Mdn/dM$ and not for the comoving number density $dn/dM$ — the two differing only by a factor of M.

### 2.1.3 Power-law distribution

The power-law distribution [11] is a less refined version of Eq. (2). It also derives from scale-invariant primordial density fluctuations and is given by

$$dn = AM^{\gamma-2}dM, \quad (6)$$

where $\gamma = -2w/(1+w)$ and $w$ is defined through the equation of state of the dominating energy in the Universe at the epoch of PBH formation such as $P = w\rho$.

### 2.1.4 Critical collapse distribution

The critical collapse distribution [11] derives from a Dirac power spectrum for primordial density fluctuations. It is defined as

$$dn = AM^{1.85} \exp\left(-\left(\frac{M}{M_f}\right)^{2.85}\right) dM, \quad (7)$$

where A is an amplitude factor and $M_f$ an upper cut-off.

### 2.1.5 Dirac distribution

The Dirac distribution simulates a Dirac BH mass function. It is useful to perform time-dependent monochromatic analyses and checks for a single BH. It is normalized to 1 BH per comoving $\text{cm}^3$.

## 2.2 Hawking radiation

### 2.2.1 Schwarzschild Black Holes

Schwarzschild Black Holes are the simplest form of BHs. They are spherically symmetric and only described by their mass M. Hawking has shown [16] that BH horizons emit elementary particles as blackbodies with a temperature linked to their mass M through<fnref n="1" />

$$T = \frac{1}{8\pi M}. \quad (8)$$

The number of particles emitted per units of time and energy is

$$\frac{d^2 N}{dt dE} = \sum_{\text{dof}} \frac{\Gamma_s(E, M)/2\pi}{e^{E/T} \pm 1}, \quad (9)$$

where the sum is over the number of quantum dof (see Table 2 in Appendix C) and the $\pm$ are for fermions and bosons, respectively. The factor $\Gamma_s$ is called the greybody factor and is detailed below.

The time-dependent comoving density of Hawking elementary particle i emitted by a distribution of BHs per units of time and energy is then computed through the integral

$$\frac{d^2 n_i}{dt dE} = \int_{M_{\text{min}}}^{M_{\text{max}}} \frac{d^2 N_i}{dt dE} \cdot \frac{dn}{dM} dM. \quad (10)$$

To obtain instantaneous quantities for a single BH of mass $M_0$, one just needs to take

$$\frac{dn}{dM} = \delta(M - M_0). \quad (11)$$

The greybody factors describe the probability that an elementary particle generated by thermal fluctuations of the vacuum at the BH horizon escapes its gravitational well. Starting from Dirac (spin $s = 1/2$) and Proca (integer spin $s$) wave equations for a particle of rest mass $\mu$

$$(i\not\partial - \mu)\psi = 0, \quad (12)$$

$$(\square + \mu^2)\phi = 0, \quad (13)$$

<footnote n="1">¹ We recall that the Newton constant G has been set to 1.</footnote>
in the Schwarzschild metric

$$ds^2 = h(r)dt^2 - h(r)^{-1}dr^2 - r^2(d\theta^2 + \sin(\theta)^2d\phi^2), \quad (14)$$

where $h(r) \equiv 1 - r_H/r$ and $r_H \equiv 2M$ is the Schwarzschild radius, Teukolsky & Press have shown [17, 18] that the wave equation can be separated into a radial equation and an angular equation if the wave is decomposed into spin weighted spherical harmonics $S_{sl}(\theta)$ and a radial component $R_s(r)$. The radial component of the master equation is for all spins $s$ [19]

$$\frac{1}{\Delta^s} \frac{d}{dr} \left( \Delta^{s+1} \frac{dR_s}{dr} \right) + \left( \frac{K^2 + i s (2r - r_H)K}{\Delta} - 4i s E r \right) - \lambda_{sl} - \mu^2 r^2 \right) R_s = 0, \quad (15)$$

where $\Delta(r) \equiv r^2h(r)$, $K(r) \equiv r^2E^2$ and $E$ is the particle frequency (or equivalently its energy). In this equation, the separation constant $\lambda_{ls} \equiv l(l+1) - s(s+1)$ is the eigenvalue of the angular equation, where $l$ denotes the angular momentum of the spherical harmonics.

To obtain the greybody factors, one has to compute the transmission coefficients of the wave between the BH horizon and the spatial infinity. The cross-section $\sigma(E)$ of the spherical wave on the BH is a sum on all spherical modes $l$ obtained through the optical theorem. The greybody factor is finally given by [20]

$$\Gamma(E, M) = \frac{\sigma(E, M)E^2}{\pi}. \quad (16)$$

The method used in BlackHawk to compute those greybody factors is described in Appendix B.1.

### 2.2.2 Kerr Black Holes

Kerr Black Holes are an extension of the Schwarzschild ones with an additional parameter: their spin $a \equiv J/M \in [0, M]$ (in the following we will denote the reduced spin parameter by $a^* \equiv a/M \in [0, 1]$) where $J$ is the BH angular momentum. These rotating BHs could gain their spin through their formation mechanism [21], accretion [22] or merging process [23]. They are axially symmetric and require a specific treatment.

The temperature of a rotating BH is given by [24]

$$T \equiv \frac{1}{2\pi} \left( \frac{r_+ - M}{r_+^2 + a^2} \right), \quad (17)$$

where $r_+ \equiv M + \sqrt{M^2 - a^2}$ is the Kerr external radius. The Teukolsky equation (15) has to be modified with $\Delta(r) \equiv r^2 - 2Mr + a^2$ and $K(r) \equiv (r^2 + a^2)E^2 + am$, where $m$ is the projection of the angular momentum $l$. The separation constant $\lambda_{slm}$, now resulting from the angular solution for

spheroidal harmonics, is more difficult to compute. We will use the 5th order expansion in terms of $\gamma = a^*ME$, as given in [24].<fnref n="2" />

The number of particles emitted per units of time and energy is now

$$\frac{d^2 N}{dt dE} = \sum_{\text{dof.}} \frac{\Gamma_s(E, M, a^*)/2\pi}{e^{E/T} \pm 1}, \quad (18)$$

where $E' \equiv E - m\Omega$ and $\Omega \equiv a^*/(2r_+)$ is the angular velocity at the horizon [24].

The method used to compute these greybody factors in BlackHawk is also described in Appendix B.1.

### 2.2.3 Exotic Black Holes

There are numerous other types of BHs, either in the classical standard cosmological model framework, such as the charged Reissner-Nordström BHs which possess a $U(1)$ electric charge (e.g. [25, 26]), or in alternative models such as (A)ds BHs [27–30], scalar-tensor theories [31–33], higher-dimensional theories [34–38], massive gravity [39, 40], ... These BHs still exhibit a Hawking radiation process in most cases, the two main differences being Hawking temperature and greybody factors. Equations giving these quantities for specific cases can usually be found in the associated literature. Possible implementations of beyond-standard BHs in BlackHawk are described in Sect. 9.5.

## 2.3 Black Hole evolution

### 2.3.1 Schwarzschild Black Holes

Once the greybody factors are known, it is possible to integrate Eq. (9) to obtain a differential equation for the mass loss of a BH through HR [41]

$$\frac{dM}{dt} = -\frac{f(M)}{M^2}. \quad (19)$$

The Page factor $f(M)$ accounts for the number of quantum dofs that a BH of mass $M$ can emit. It is obtained through [41]

$$f(M) = -M^2 \frac{dM}{dt} = M^2 \int_0^{+\infty} \frac{E}{2\pi} \sum_{\text{dof}} \frac{\Gamma_s(E, M)}{e^{E/T} \pm 1} dE. \quad (20)$$

The computation of the $f(M)$ factor in BlackHawk is described in Appendix B.2.

<footnote n="2">² With our conventions we have an opposite sign for $\gamma$ compared to Ref. [24] (all odd-terms in their Appendix A have to be switched in sign).</footnote>
### 2.3.2 Kerr Black Holes

For Kerr BHs, a new phenomenon arises. The rotation of the BH enhances the emission of particles with high angular momentum, and with a projection *m* of that angular momentum aligned with the BH spin, thus effectively extracting angular momentum from the BH [42]. The equation for the Page factor *f*(*M*, *a*<sup>*</sup>) becomes [24,43,44]

$$ 
\begin{aligned}
f(M, a^*) &= -M^2 \frac{dM}{dt} \\
&= M^2 \int_0^{+\infty} \frac{E}{2\pi} \sum_{\text{dof}} \frac{\Gamma_s(E, M, a^*)}{e^{E/T} \pm 1} dE,
\end{aligned} 
\quad (21) 
$$

and the differential equation describing the angular momentum *J* is [24,43,44]<sup>3</sup>

$$ 
\begin{aligned}
g(M, a^*) &= -\frac{M}{a^*} \frac{dJ}{dt} \\
&= -\frac{M}{a^*} \int_0^{+\infty} \sum_{\text{dof}} \frac{m}{2\pi} \frac{\Gamma_s(E, M, a^*)}{e^{E/T} \pm 1} dE.
\end{aligned} 
\quad (22) 
$$

Once the *f*(*M*, *a*<sup>*</sup>) and *g*(*M*, *a*<sup>*</sup>) Page factors are obtained, the evolution of *a*<sup>*</sup> is straightforwardly obtained through

$$ 
\begin{aligned}
\frac{da^*}{dt} &= \frac{d(J/M^2)}{dt} \\
&= \frac{1}{M^2} \frac{dJ}{dt} - 2 \frac{J}{M^3} \frac{dM}{dt} \\
&= a^* \frac{2f(M, a^*) - g(M, a^*)}{M^3}.
\end{aligned} 
\quad (23) 
$$

The computation of the *f*(*M*, *a*<sup>*</sup>) and *g*(*M*, *a*<sup>*</sup>) Page factors in BlackHawk is described in Appendix B.2.

### 2.3.3 Exotic Black Holes

Exotic BHs listed in Sect. 2.2.3 can have a modified evolution as compared to the Schwarzschild and Kerr cases, for two main reasons. First, since their greybody factors and temperature are different, the *f* and *g* parameters are expected to be different as well and the master Eq. (18) will give a different emission rate. Second, these BHs can possess other scalar degrees of freedom, such as a *U* (1) charge (e.g. the electric charge in Reissner–Nordström BHs [25]), which experience a specific evolution. Evolution equations for these additional charges have to be derived, and would be similar to Eqs. (21) and (22). The implementation of beyond-standard BHs in BlackHawk is described in Sect. 9.5.

## 2.4 Hadronization

The elementary particles emitted by BHs are not the final products of the HR. Some of them are unstable, others only exist in hadrons. A particle physics code has to be used in order to evolve the elementary particles into final products. We used HERWIG [45] and PYTHIA [46] for this purpose.

The final particles, hereby denoted as “secondary Hawking particles” (the elementary being the “primary Hawking particles”), depend on the cosmological context in which they are emitted. For Big-Bang Nucleosynthesis (BBN) studies, an estimation of the reaction rates imposes to keep the particles with a lifetime longer than ~ 10<sup>-8</sup> s [47]. These particles are listed in the Table 2 of Appendix C.

The time-dependent comoving density of Hawking secondary particle *j* emitted by a distribution of BHs per units of time and energy is computed with the integral

$$ \frac{d^2 n_j}{dt dE} = \int \sum_i \frac{d^2 n_i}{dt dE'} \cdot \frac{dN_j^i}{dE} dE', \quad (24) $$

where the sum is taken over Hawking primary particles *i*, and Sect. Appendix B.3 describes how hadronization tables $dN_j^i(E', E)$ have been computed to transform the primary spectra into secondary spectra in BlackHawk.

# 3 Content and compilation

This section describes the structure and file content of the code and explains its usage. BlackHawk is written in C and has been tested under Linux, Mac and Windows (using Cygwin64). It can be obtained from

blackhawk.hepforge.org

## 3.1 Main directory

The main directory contains:

- the source codes BlackHawk_*.c containing the main routines,

- a pre-built parameter file parameters.txt,

- a compilation file Makefile,

- a README.txt file containing general information about the code,

- four folders src/, results/, manual/ and scripts/ that are described in the following.

## 3.2 src/ sub-folder

This folder contains:

<footnote n="3"><sup>3</sup> Same remark as above, we have in our conventions an opposite sign for *g*.</footnote>
- a header file `include.h` containing the declaration of all routines along with the parameter structure `struct param` (see Sect. 4.1) and the numerical values of general quantities (units conversion factors, constants, particle masses...),

- ten source files containing the definition of all the BlackHawk routines:

- `evolution.c`,

- `general.c`,

- `hadro_herwig.c`,

- `hadroPythia.c`,

- `hadroPythianew.c`,

- `primary.c`,

- `secondary.c`,

- `spectrum.c`,

- `technical.c`,

- two compilation files `Makefile` and `FlagsForMake`,

- a subfolder `tables/` containing all the numerical tables which will be described in the following.

### 3.3 results/ sub-folder

This folder is designed to receive sub-folders of data generated by running the BlackHawk code (see Sect. 7).

### 3.4 manual/ sub-folder

This folder contains an up-to-date version of the present manual.

### 3.5 scripts/ sub-folder

This folder contains all the scripts used to compute the numerical tables mentioned in the following, as well as visualization scripts and a main program for SuperIso Relic [48–50]. These scripts can be used to generate the needed tables. They are accompanied by `README.txt` files explaining how to use them.

### 3.6 Compilation

The compilation of BlackHawk has been tested on Linux, Mac and Windows (using Cygwin64) distributions. The code is written in C99 standard. To compile the code, simply `cd` into the main directory and type<fnref n="4" />:

```
make BlackHawk_*
```

where `*` denotes `tot` or `inst`. This will create a library file `libblackhawk.a` and an executable `BlackHawk_*.x`. The compiler and compilation flags can be modified in `Makefile` if needed.

To run the code, `cd` to the main directory and type<fnref n="5" />:

```
./BlackHawk_*.x parameter_file
```

where `parameter_file` is the name of a parameter file. To compile only the library, just `cd` into the main directory and type:

```
make
```

## 4 Input parameters

In this section we describe how input parameters are handled in BlackHawk and their meaning.

### 4.1 Parameter structure

The input parameters used by BlackHawk are listed in the `parameters.txt` file. This file can be modified by the user and is saved for each new run of the code in the destination directory. A C structure has been defined in `include.h` to embed all the parameters:

```c
struct param {
    char destination folder[32];
    int full_output;
    int interpolation_method;

    int BHnumber;
    double Mmin;
    double Mmax;
    double a;
    int spectrum-choice;
    double amplitude;
    double variance;
    double crit-mass;
    double eqstate;
    char table[32];

    double tmin;
    int nb_fin_times;
    int limit;
    double Mmin_fM;
    double Mmax_fM;
    double amin_fM;
    double amax_fM;
```

<footnote n="4">⁴ In case of problems of memory size at compilation, editing `src/include.h` and commenting `#define HARDTABLES` can solve the problem at the price of a longer execution time.</footnote>

<footnote n="5">⁵ In case of memory problem at execution, increasing the stack size with the command `ulimit -s unlimited` can help solving the problem.</footnote>
int nb_fM masses;
int nb_fM_a;

int Enumer;
double Emin;
double Emax;
int particle_number;
int grav;
int nb_a;
int nb_x;

int primary_only;
int hadronization-choice;
double Emin_hadro;
double Emax_hadro;
int nb_init_en;
int nb_fin_en;
int nb_init_part;
int nb_fin_part;

];

Most routines described in Sect. 5 will use this structure as an argument in order to have an easy access to the run parameters. Depending on the choices of the parameters, some parameters can be irrelevant for a given run and will therefore not be taken into account, and no error message will be displayed for the irrelevant/unused parameters.

## 4.2 General parameters

The first set of parameters defines the general variables:

- `destination folder` is the name of the output folder that will be created in `results/` to save run data.

- `full_output` determines whether the shell output will be expanded (1) or not (0). It can be useful to debug the code or to see the progress in time-consuming routines.

- `interpolation_method` determines whether the interpolations in the tables are made linearly (interpolation between the tabulated values) or logarithmically (linear interpolation between the decimal logarithm of the tabulated values).

## 4.3 Black Hole spectrum parameters

The second set of parameters defines the quantities used to compute the BH density distribution (see Sect. 5.2):

- `BHnumber` is the number of BH masses that will be simulated. If the parameter `spectrum-choice` is not set to 5, it has to be an integer greater than or equal to 1. If it is equal to 1, the only BH mass will be `Mmin` (see

below). If the parameter `spectrum-choice` is set to 5, it has to be the number of tabulated values in the user-defined BH distribution (see below and Sect. 9.2). It will be automatically set to 1 if `spectrum-choice` is set to 0.

- `Mmin` and `Mmax` are respectively the lowest and highest BH masses that will be simulated. They have to be given in grams and satisfy the condition $M_p \approx 2 \times 10^{-5} \text{ g} < M_{min}$, $M_{max}$, where $M_p$ is the Planck mass. For a mass distribution, one must have $M_{min} < M_{max}$. If they are not compatible with boundaries of the mass distribution, the computation will stop (see below).

- `spectrum-choice` selects the form of the BH mass distribution (see Sect. 2.1). It has to be an integer among 0 (Dirac, mimicking a single BH<fnref n="6" />), 1 (log-normal), 2 (power-law), 3 (critical collapse), 4 (peak theory) and 5 (user-defined distribution, see below and Sect. 9.2).

- `amplitude` is the amplitude $A$ present in Eqs. (5)–(7). It is the normalization of the corresponding BH distribution and thus strictly positive.

- `variance` is the variance $\sigma$ in the log-normal distribution of Eq. (5). It has to be strictly positive.

- `crit-mass` is the characteristic mass $M_c$ in Eq. (5) and $M_f$ in Eq. (7). It has to be strictly positive.

- `eq_state` defines the equation of state $w$ (see Sect. 2.1.3).

- `table` is the name of a user-defined BH distribution table. It has to be a string with any file extension.

## 4.4 Black Holes evolution parameters

The next set of parameters defines the quantities used to compute the BHs evolution (see Sect. 5.3):

- `tmin` is the initial integration time of the evolution of BH, in seconds. It can have any positive value, but we recommend that it is lower than the lifetime of the lightest BH under consideration.

- `nb_fin_times` is the number of final integration times that will be used in the computations. It will be set automatically by the integration procedure.

- `limit` is the iteration limit when computing the time evolution of a single BH (see Sect. 5.3). It is fixed to $\text{limit} = 5000$ even if the effective iteration numbers hardly reach 1000. It should be increased if the integration does not reach the complete evaporation of BHs.

- `Mmin_fM` and `Mmax_fM` are the BH mass boundaries used to compute the $f(M, a^*)$ and $g(M, a^*)$ tables. They should not be modified unless the user recomputes the corresponding tables (see Sect. 9).

<footnote n="6">⁶ That is to say, the emissivities obtained are those of a single BH. This option can be useful to compute known test emissivities of single BHs.</footnote>
- `amin_fM` and `amax_fM` are the BH spin boundaries used to compute the $f(M, a^*)$ and $g(M, a^*)$ tables. They should not be modified unless the user recomputes the corresponding tables (see Sect. 9).

- `nb_fM_masses` and `nb_fM_a` are respectively the number of BH masses and spins tabulated in the $f(M, a^*)$ and $g(M, a^*)$ tables. They should not be modified unless the corresponding tables are recomputed (see Sect. 9).

## 4.5 Primary spectrum parameters

This set of parameters defines the quantities related to the primary Hawking spectra (see Sect. 5.4):

- `Emin` and `Emax` are the minimum and maximum primary particle energies, respectively. They must be compatible with the table boundaries (see below) and satisfy $0 < Emin < Emax$.

- `Enumber` is the number of primary particles energies that will be simulated. It has to be an integer greater than or equal to 2.

- `particle_number` is the number of primary particle types. It is fixed to 15 (photon, gluon, $W^\pm$ boson, $Z^0$ boson, Higgs boson, neutrino, 3 leptons (electron, muon, tau) and 6 quarks (up, down, charm, strange, top, bottom)) and should not be modified unless the user recomputes the primary particle table (see Sect. 9).

- `grav` determines whether the emission of gravitons by BH will be taken into account (`grav = 1`) or not (`grav = 0`).

- `nb_x` and `nbgamma_x` are respectively the number of spins $a^*$ and values of $x \equiv 2 \times E \times M$ tabulated in the greybody factor tables. They should not be modified unless the corresponding tables are recomputed (see Sect. 9).

## 4.6 Hadronization parameters

This last set of parameters defines the quantities used during the hadronization (see Sect. 5.5):

- `primary_only` determines whether the secondary spectra will be computed or not. It has to be an integer between 0 (primary spectra only) and 1 (primary and secondary spectra). In the case where the parameters `Emin` and `Emax` are not compatible with the hadronization table boundaries (see below), a warning will be displayed and extrapolation used.

- `hadronizationChoice` determines which hadronization tables will be used to compute the secondary spectra (see Sect. Appendix B.3). It has to be an integer between 0 (`PYTHIA` tables – Early Universe/BBN epoch), 1 (`HERWIG` tables – Early Universe/BBN epoch) and 2 (new `PYTHIA` tables – present epoch).

- `Emin_hadro` and `Emax_hadro` are the energy boundaries of the hadronization tables. They should not be changed unless the user recomputes the corresponding tables (see Sect. 9).

- `nb_init_en` and `nb_fin_en` are the number of initial and final particle energy entries in the selected hadronization tables, respectively. They should not be modified unless the corresponding tables are recomputed (see Sect. 9).

- `nb_init_part` and `nb_fin_part` are the number of primary and secondary particle types in the selected hadronization tables, respectively. They should not be modified unless the corresponding tables are recomputed (see Sect. 9).

# 5 Routines

Below are listed the main routines defined in `BlackHawk`. To simplify the analytic formulas, all intermediate quantities are in GeV (see Appendix A for conversion rules).

## 5.1 General routines

There are 4 general routines in the `BlackHawk` code. The principal ones are the main routines, described in Sect. 6. The other two are:

- `int read_params(struct param *parameters, char name[], int session)`: this routine reads the file `name`. The parameters are converted from `CGS` units to GeV. The user should respect the original syntax when modifying the parameters (concerning spaces, underscores, ...), except for comments which are preceded by a `#` symbol. It takes a pointer to a `struct param` object (see Sect. 4.1) as an argument and fills it using the file `name`. The argument `session` shows which of the main program has been launched (0 for `BlackHawk_tot`, 1 for `BlackHawk_inst`). If one parameter is not of the type described in Sect. 4 this function will display an error message. Any of these errors will end the `BlackHawk` run. If one parameter is in small contradiction with the others but the computation can still be partly done (e.g. only the primary spectra can be computed with the given parameters) a warning message will be displayed. In such case, the problematic parameters will be set automatically (e.g. `primary_only = 1`) and the computation will proceed.

- `int memory_estimation(struct param *parameters, int session)`: this routine gives a rough estimate of the usage of both RAM and disk
space (see Sect. 8). If the user decides to cancel the run
the value 0 is returned, otherwise it is 1. The output is
given in MB.

## 5.2 Black Hole spectrum routines

There are 4 routines contributing to the BH initial spectrum computation (see Sect. 2.1):

- void read_users_table(double *init_masses, double *init_spins, double *spec_table, struct param *parameters): this routine reads a user-defined BH distribution table in the file given by the parameter table. It fills the arrays init_masses[], init_spins[] and spec_table[] with results converted from CGS units to GeV.

- double nu(double M): this routine takes a BH mass as an argument and computes the dimensionless quantity $\nu(M)$ defined in Eq. (3).

- double n_cov(double M, double *table_masses, double *table codensities, int index, struct param *parameters): this routine takes a BH mass as an argument and computes the comoving density $\frac{dn}{dM}$ defined in Eq. (2) (using the nu routine) or (5) or (6) or (7) (in GeV$^2 \rightarrow$ cm$^{-3}$.g$^{-1}$), depending on the parameter spectrum-choice (see Sect. 4.3). If this parameter is set to 0, a flat distribution is used with only one BH mass, mimicking a Dirac distribution normalized to one BH per comoving cm$^3$.

- void spectrum(double *init_masses, double *init_spins, double *spec_table, double *table_masses, double *table codensities, struct param *parameters): this routine fills the array init_masses[] with BHnumber BH masses logarithmically distributed between Mmin and Mmax. If the parameter BHnumber is set to 1, the only BH initial mass will be Mmin. For each BH mass, it then fills the array init_spins[] with a spin a (the same for each mass) and the array spec_tables[] computing the corresponding comoving densities $dn$ (in GeV$^3 \rightarrow$ cm$^{-3}$) using the n_cov routine where $dM$ is taken around the considered mass. The result is rescaled by a factor $10^{100}$ due to the very small numbers involved in the dimensionless computation.

- void write_spectrum(double *init_masses, double *init_spins, double *spec_table, struct param *parameters): this routine writes the BH initial masses, spins and comoving densities in a file BH_spectrum.txt, saved in destination folder (see Sect. 7.1).

The results are converted from GeV to CGS units.

## 5.3 Black Holes evolution routines

There are 6 routines contributing to the BH time evolution computation (see Sect. 2.3):

- double rplus_BH(double M, double a): this routine gives the external Kerr radius of a rotating BH for a given mass $M$ and spin $a^*$ (see Sect. 2.2.2) (in GeV$^{-1} \rightarrow$ cm);

- double temp_BH(double M, double a): this routine gives the Hawking temperature of a Kerr BH for a given mass $M$ and spin $a^*$ using Eq. (17) (in GeV $\rightarrow$ K).

- void read_fM_table(double **fM_table, double *fM_masses, double *fM_a, struct param *parameters): this routine reads the $f(M, a^*)$ factor (see Eq. (21)) in the table contained in the folder fM_tables/ (see Sect. Appendix B.2). It fills the arrays fM_masses[] (in GeV $\rightarrow$ g), fM_a[] and fM_table[][] (in GeV$^4 \rightarrow$ g$^3$.s$^{-1}$).

- void read_gM_table(double **gM_table, double *fM_masses, double *fM_a, struct param *parameters): this routine reads the $g(M, a^*)$ factor (see Eq. (22)) in the table contained in the folder fM_tables/ (see Sect. Appendix B.2). It fills the arrays fM_masses[] (in GeV $\rightarrow$ g), fM_a[] and gM_table[][] (in GeV$^4 \rightarrow$ g$^2$.GeV.s$^{-1}$).

- double loss_rate_M(double M, double a, double **fM_table, double *fM_masses, double *fM_a, int counter_M, int counter_a, struct param *parameters): this routine computes the quantity $\frac{dM}{dt}$ defined in Eq. (21) (in GeV$^2 \rightarrow$ g.s$^{-1}$).

- double loss_rate_a(double M, double a, double **fM_table, double **gM_table, double *fM_masses, double *fM_a, int counter_M, int counter_a, struct param *parameters): this routine computes the quantity $\frac{da^*}{dt}$ defined in Eq. (23) (in GeV $\rightarrow$ s$^{-1}$).

- void life_evolution(double **life_masses, double **life_spins, double *life_times, double *dts, int *evolution_length, double *init_masses, double *init_spins, double **fM_table, double **gM_table, double *fM_masses, double *fM_a, struct param *parameters): this routine computes the evolution of each of the initial BH masses in init_masses[] and BH spins in destination folder (see Sect. 7.1).
init_spin[s]. The initial time life_times[0] is set to tmin, the initial masses life_masses[i][0] are set to init_masses[i] and the initial spins life_spin[i][0] are set to init_spin[i]. Iteratively, the next masses and spins are estimated using the Euler method

$$M(t + dt) = M(t) + \frac{dM}{dt}dt, \qquad (25)$$

$$a^*(t + dt) = a^*(t) + \frac{da^*}{dt}dt, \qquad (26)$$

where the derivatives are computed using the loss\_rate\_* routines. If one of the relative variations is too large (|dX/X| > 0.1) then the time interval is divided by 2. If all the variations are very small (|dX/X| < 0.001), and if the current timestep is reasonable compared to the current timescale (dt/t ≲ 1) then the time interval is multiplied by 2. Once the dimensionless spin reaches 10⁻³, we stop computing its variation and simply set it to 0, and it does not enter anymore in the adaptive timesteps conditions. This goes on until each mass reaches the Planck mass or the recursion limit limit × BHnumber is attained, in which case the following error is displayed

[life_evolution] : ERROR ITERATION LIMIT REACHED !

This may be a sign that the parameter limit should be increased. The intermediate time intervals dt, times t, masses M and spins a* are stored in the arrays dts[], life_times[] (both in GeV⁻¹ → s), life_masses[][] (in GeV → g) and life_spin[][], respectively. The number of intermediate iterations for each initial mass is stored in the array evolution_length[].

- void write_life_evolutions(double **life_masses, double **life_spin, double **life_times, int *evolution_length, struct param *parameters): this routine writes the BH time-dependent masses and spins until full evaporation in the file life_evolutions.txt, saved in destination folder (see Sect. 7.1). The results are converted from GeV to CGS units.

## 5.4 Primary spectra routines

There are 5 routines contributing to the computation of the primary Hawking spectra (see Sect. 2.2):

- void read_gammas, double *gamma_a, double *gamma_x, struct param *parameters): this routine reads the quantities Γ/(eE/T ± 1), defined

in Eq. (18), in the tables spin\_*.txt in the folder gamma_tables/. It fills the arrays gamma_a[] and gamma_x[] with the tabulated spins a* (dimensionless) and x ≡ ErBH (dimensionless → GeV·cm), respectively. It fills the array gammas[] [] [] with the corresponding dimensionless greybody factors in format [type][spin][x] (see Appendix B.1).

- void read_asymp_fits(double ***fits, struct param *parameters): this routine reads the asymptotic fit parameters for the greybody factors, contained in the tables spin\_\_*\_fits.txt in the folder gamma_tables/. It fills the array fits[] [] [] in format [type][spin][parameters] (see Appendix B.1).

- double dNtdtE(double E, double M, double a, int particle_index, double ***gammas, double *gamma_a, double *gamma_x, double ***fits, double *dof, double *spins, double *masses_primary, int counter_a, int counter_x, struct param *parameters): this routine computes the emission rate d²N/dtdE of the primary particle particle_index (see Eq. (18)), for a given particle energy E, the BH mass M, the BH spin a and the particle informations contained in dof[], spins[] and masses_primary[]. If x ≡ ErBH is in the greybody factor boundaries, the values are interpolated in those tables at position counter_a and counter_x. Otherwise, we use the asymptotic fits tables (see Appendix B.1). The result is dimensionless (→ GeV⁻¹·s⁻¹).

- void instantaneous_primary_spectrum (double **instantaneous_primary_spectra, double *BH_masses, double BH_spin, double *spec_table, double *energies, double ***gammas, double *gamma_a, double *gamma_x, double ***fits, double *dof, double *spins, double *masses_primary, struct param *parameters): this routine computes the instantaneous primary Hawking spectra for a distribution of BHs given by the routine spectrum, namely the quantity d²n/dtdE in Eq. (10) for each primary particle and each energy in energies[], computed with the routine dNtdtE. The results are stored in the array instantaneous_primary_spectra[] [] in format [particle][energy].

- void write_instantaneous_primary_spectra (double **instantaneous_primary_spectra, double *energies, struct param *parameters): this routine writes the instantaneous primary Hawking spectra in a file instantaneous_primary_spectra.txt, saved
in `destination folder`/ (see Sect. 7.2). The results are converted from GeV to CGS units.

## 5.5 Secondary spectra routines

There are 9 routines contributing to the computation of the secondary Hawking spectra (see Sect. 2.4):

- `void convertHadronization_tables` (double ****tables, double *initial_energies, double *final_energies, struct param *parameters):` this routine is auxiliary. It writes hardcoded versions of the hadronization tables (see Appendix B.3) in files `hadronization_tables_*.h` in the `tables/` subfolder in order to accelerate the code execution, while slowing its compilation.

- `void readHadronization_tables` (double ****tables, double *initial_energies, double *final_energies, struct param *parameters):` this routine reads the hadronization table (see Appendix B.3) determined by `hadronization-choice`. If `HARDTABLES` is defined, it uses the table included at compilation using the routines `readHadronization_*,` otherwise it reads the corresponding table in the `tables` subfolder. It fills the arrays `initial_energies[]` and `final_energies[]` with the tabulated primary particles and secondary particles energies (in GeV), respectively, and fills the array `tables[][][]` with the corresponding branching ratios $\frac{dN_j^i}{dE'}$ in Eq. (24) (in GeV$^{-1}$) in format `[secondary particle][initial energy][final energy][primary particle]`.

- `void totalSpectra` (double ***partial_hadronized_spectra, double **partial_primary_spectra, double **partial_integrated_hadronized_spectra, double ****tables, double *initial_energies, double *final_energies, double ***primary_spectra, double *times, double *energies, double *masses_secondary, struct param *parameters):` this routine is a container that uses the “instantaneous” routines to compute the Hawking primary and secondary spectra at each timestep in `times` and writes it directly in the output in order to save RAM memory. To do so, it creates the output files `*_primary_spectrum.txt` and `*_secondary_spectrum.txt` (if `primary_only` is set to 0). Then, it fills the partial arrays `partial_*` with the instantaneous primary spectra, hadronized spectra and integrated spectra. Finally, it calls the routine

`write_lines` to write the partial result in the output before moving to the next timestep.

- `void write_lines` (char **file_names, double **partial_integrated_hadronized_spectra, double time, struct param *parameters):` given a `time` and instantaneous primary and secondary spectra (if `primary_only` is set to 0), this routine writes a new line in the `*_spectrum.txt` files. The arrays `write_*[]` determine whether the values of each particles are written or not, thus potentially saving disc memory. Results are converted from GeV to CGS units (see Sect. 7.1).

- `double contribution_instantaneous` (int j, int counter, int k, double **instantaneous_primary_spectra, double ****tables, double *initial_energies, double *final_energies, int particle_type, int hadronization-choice):` this routine computes the instantaneous integrated of Eq. (24) (in GeV$^{-1} \rightarrow GeV^{-2} \cdot s^{-1}$) for the secondary particle `particle_type`, initial energy $E' = energies[j]$, corresponding tabulated initial energy `initial_energies(counter]` and final energy $E = final_energies[k]$. The sum over channels of production of the secondary particles may depend on the structure of the hadronization tables.

- `void hadronize_instantaneous` (double ***instantaneous_hadronized_spectra, double ****tables, double *initial_energies, double *final_energies, double **instantaneous_primary_spectra, double *energies, struct param *parameters):` this routine computes the instantaneous secondary Hawking spectra for all secondary particles, all initial energies in `energies[]` and all final energies in `final_energies[]`. It fills the array `instantaneous_hadronized_spectra[][][]` using the routine `contribution_instantaneous`, in format `[secondary particle][initial energy][final energy]`. If the initial energy is not in the hadronization tables, the contribution is extrapolated.

- `void integrate_initial_energies_instantaneous` (double ***hadronized emission_spectra, double **integrated_hadronized_spectra, double *energies, double *final_energies, struct param *parameters):` this routine computes the integral Eq. (24) (dimensionless $\rightarrow GeV^{-1} \cdot s^{-1}$) using the `trapeze` routine. The results are stored in the array `instantaneous_integrated_hadronized_spectra[][]` in format `[secondary particle][final energy]`.
- `void add_* _instantaneous (double **instantaneous_primary_spectra, double **instantaneous_integrated_ hadronized_spectra, double *energies, double *final_energies, struct param *parameters)`: these two routines add the contribution of the primary photons/neutrinos to the secondary produced ones. The value in term of final energies is interpolated in the primary spectrum and added to the hadronized spectrum `instantaneous_integrated_hadronized_spectra[][]`.

- `void write_instantaneous_hadronized_spectra(double **instantaneous_integrated_hadronized_spectra, double *hadronized_energies, struct param *parameters)`: this routine writes the instantaneous secondary Hawking spectra in the file `instantaneous_secondary_spectra.txt`, saved in destination folder / (see Sect. 7.2). The results are converted from GeV to CGS units.

## 5.6 Auxiliary routines

8 auxiliary routines are used throughout the code:

- `double trapeze(double x1, double x2, double y1, double y2)`: this routine performs the trapeze integration of a function f that takes values y1 in x1 and y2 in x2 using

$$\int_{x1}^{x2} f(x)dx \approx \frac{1}{2}(x2 - x1) \times (y1 + y2). \quad (27)$$

- `void free* (*)`: these routines perform a proper memory freeing of n-dimensional arrays of various types, by recursively applying the native free routine.

- `int ind_max(double *table, int llength)`: this routine returns the index of the maximum of the array `table[]` of length `llength`.

## 6 Programs

The BlackHawk code is split into two programs, which are presented in this section:

- `BlackHawk_tot`: full time-dependent Hawking spectra;

- `BlackHawk_inst`: instantaneous Hawking spectra.

Once a set of parameters is chosen, the two programs can be launched in the same `destination folder/` because the output files will not enter in conflict (see Sect. 7). We

will now describe the structure of the main routines together with screen output examples.

### 6.1 Common features

When running the BlackHawk code, some routines will be called regardless of the program choice. First, some general quantities are fixed (which are converted into GeV when applicable, see Appendix A):

- `machine_precision = 10⁻¹⁰` defines the precision up to which two double numbers are considered as equal.

- $G = 6.67408 \times 10^{-11} \text{ m}^3 \cdot \text{kg}^{-1} \cdot \text{s}^{-2}$ is the Newton constant in SI units.

- $M_p \equiv G^{-1/2}$ is the Planck mass in the natural system of units.

- `m_*` are the masses of the Standard Model particles (see Table 2 in Appendix C).

- `*_conversion` are the quantities used to convert units from CGS/SI to GeV (see Appendix A).

The code runs in several steps, which are separated on the output screen. A new step starts with:

`[main] : ***** ...`

and ends with:

### DONE

If the `full_output` parameter is set to 1, then more information will be displayed about the progress of the steps. In the case where information appears with the name of another routine inside brackets, it means that an error occurred.

The first common step is the definition and filling of the parameters structure using `read_params`. Then an estimation of the memory that will be used is displayed by `memory_estimation`. The user can choose to go on or to cancel the run (see Sect. 5.1). If no error was found in the input parameters, the output directory `destination folder/` is created. If it already exists, the user has the choice to overwrite the existing data or to stop the execution in order to choose another output folder. For a subsequent data interpretation, the parameters file is copied in the output folder. The expected output at this stage is of the form<fnref n="7" />:

<footnote n="7">⁷ No user checking will be done if `CHECK_USER` is defined to 0 in `include.h`.</footnote>
########## BLACKHAWK v1.0
##########
# HAWKING SPECTRUM
# COMPUTATION DEVICE

[main] : STARTING EXECUTION...
[main] : READING THE RUN PARAMETERS IN
'parameters.txt'... DONE
[main] : ESTIMATION OF THE MEMORY
USE...

Running this session will use at least 101.603 MB of RAM and 3.408 MB of disc memory.

Do you want to continue? (type y or n) y

DONE

[main] : SAVING RUN PARAMETERS... DONE

The subsequent execution steps depend on the program. Output examples are given in the mode `full_output = 0`.

## 6.2 BlackHawk_tot: Time-dependent Hawking spectra

In this program, BlackHawk computes the time-dependent Hawking spectra of a chosen initial distribution of BHs.

BlackHawk will compute the initial distribution of BHs (at `tmin`) using the routine `spectrum` or will read the user-defined BH distribution file `table` with the routine `read_users_table` (depending on the `spectrum-` choice), filling the arrays `init masses`[, `init spins`][] and `spec_table`[]. It writes the results in the output with `write_spectrum` (see Sect. 5.2).

It then reads the $f(M, a^*)$ and $g(M, a^*)$ tables using the `read_fM_table` and `read_gM_table` routines, respectively, filling the arrays `fM_table`[[], `gM_table`[][], `fM masses`[], `fM_a`[], in order to evolve in time each initial BH spin and mass until the Planck mass limit using the routine `life_evolution`. This fills the arrays `life_times`[, `life masses`[][], `life_spin`[][], `dts`[] and `evolution_length`[]. The evolutions in time are written in the output using the routine `write_life_evolutions` (see Sect. 5.3).

Then BlackHawk reads the greybody factor tables using the `read_gamma_tables` routine, filling the arrays `gammas`[][][], `gamma_a`[] and `gamma_x`[], and the fits tables using `read_asymp_fits`, filling the array `fits`[][][]. The common time range `times`[] is filled with the times in `life_times`[] until the evaporation of

the last BH. This time range thus embeds all interesting intermediate evolution timesteps.

If the parameter `primary_only` has been set to 0, BlackHawk reads the suitable hadronization tables (depending on the `hadronization-` choice) with the routine `read_hadronization_tables`, filling the arrays `tables`[][][][], `initial_energies`[] and `final_energies`[]. It uses all these tables to compute the primary and secondary (if `primary_only = 0`) Hawking spectra using the routine `total_spectra`. Due to the large number of intermediate timesteps when a full distribution is considered, we do not perform the full computation in one step in the RAM memory, but rather do it timestep by timestep using the intermediate arrays `partial_primary_spectra`[][], `partial_hadronized_spectra`[][][] and `partial_integrated_hadronized_spectra`[][], and the instantaneous routines `hadronize_instantaneous`, `integrate_initial_energies_instantaneous` and `add_*_instantaneous`. The intermediate results are written in the output thanks to `write_lines` (see Sect. 5.5).

This is the end of the execution of `BlackHawk_tot`. The expected output is of the form:

[main] : COMPUTING THE INITIAL DISTRIBUTION OF BLACK HOLES... DONE
[main] : WRITING INTO FILE 'BH_spectrum.txt'... DONE
[main] : READING EVOLUTION TABLES... DONE
[main] : COMPUTING THE EVOLUTION OF BLACK HOLES... DONE
[main] : WRITING INTO FILE 'life_evolutions.txt'... DONE
[main] : READING GAMMA TABLES... DONE
[main] : READING FIT TABLES... DONE
[main] : READING HADRONIZATION TABLES... DONE
[main] : COMPUTING SPECTRA... DONE
[main] : END OF EXECUTION

## 6.3 BlackHawk_inst: Instantaneous Hawking spectra

In this program, BlackHawk computes the instantaneous Hawking spectra of a distribution of BHs.

First BlackHawk will compute the initial distribution of BHs (at `tmin`) using the routine `spectrum` or it will read the user-defined BH distribution file `table` with the routine `read_users_table` (depending on the `spectrum-` choice), filling the arrays `init masses`[], `init_spin`[] and `spec_table`[]. It then writes
the results in the output with `write_spectrum` (see Sect. 5.2).

Then `BlackHawk` reads the greybody factor tables using the routine `read_tables`, filling the arrays `gammas`[] [] []], `gamma_masses`[] and `gamma_energies`[] and the fit table with the routine `read_asymp_fits`, filling the array `fits`[] [] []], to compute the primary Hawking spectra using the routine `instantaneous_primary_spectrum`, filling the arrays `instantaneous_primary spectra`[] []]. The results are written in the output by the routine `write_instantaneous_primary spectra` (see Sect. 5.4).

If the parameter `primary_only` has been set to 0, `BlackHawk` reads the hadronization tables (depending on the hadronization `choice`) using the routine `read_hadronization_tables`, filling the arrays `tables`[] [] [] []], `initial_energies`[] and `final_energies`[]], and uses them to compute the secondary Hawking spectra using the routine `hadronize_instantaneous`, filling the array `instantaneous_hadronized spectra`[] [] [].

The initial energy dependence of the spectra is integrated out with the routine `integrate_initial_energies_instantaneous`, which fills the array `instantaneous_integrated_hadronized spectra`[] []]. The contributions from primary photons and neutrinos are added to the secondary spectra by the routines `add_*_instantaneous`. The results are written in the output by the routine `write_instantaneous_hadronized_spectra` (see Sect. 5.5).

This is the end of the execution of `BlackHawk_inst`. The expected output is of the form:

```
[main] : COMPUTING THE INITIAL DISTRIBUTION OF BLACK HOLES... DONE
[main] : WRITING INTO FILE 'BH_spectrum.txt'... DONE
[main] : READING GAMMA TABLES... DONE
[main] : READING FITS TABLES... DONE
[main] : COMPUTING PRIMARY SPECTRA... DONE
[main] : WRITING INTO FILE instantaneous_primary_spectra.txt... DONE
[main] : READING HADRONIZATION TABLES ... DONE
[main] : HADRONIZING PARTICLES... DONE
[main] : INTEGRATING OVER INITIAL ENERGIES... DONE
[main] : WRITING INTO FILE instantaneous_secondary_spectra.txt... DONE
[main] : END OF EXECUTION
```

# 7 Output files

As explained in the previous sections, all the output files generated by a run of `BlackHawk` will be stored in a `destination folder`/. In this section we describe the format of these files created by each program. Examples of results can be found in Appendix D. In all the cases, the parameter file `parameters.txt` used for the run is copied in the output folder in order to allow for subsequent data interpretation.

Python vizualisation scripts have been incorporated in the sub-folder `scripts/` in order to plot the data produced by both programs. They come with a file `README.txt` that explains how to configure them. You can of course modify these scripts to your own purpose or use any other plotting program.

## 7.1 BlackHawk_tot

Running `BlackHawk_tot` produces 4 (or 3) types of output files:

- `BH_spectrum.txt`: this file is written by the routine `write_spectrum`. It contains the initial density spectrum of BHs and has 3 columns: the first one is a list of the BHs initial masses (in g), the second one the corresponding list of initial spins (dimensionless) and the third one is the comoving number densities (in $cm^{-3}$).

- `life_evolutions.txt`: this file is written by the routine `write_life_evolutions`. It contains all the integrated timesteps for each initial BH mass. It includes a list of the number of integration timesteps for each initial BH mass. Also it contains a table in which the first column is the time (in s), and each other column is the evolution of the mass of a BH (in g) as a function of time. Finally it includes a table with the same format giving the evolution of the spins (dimensionless).

- `*_primary_spectrum.txt`: these files are written by the routine `write_lines`. They contain the emission rates of each primary particle at each final time and for each simulated initial energy. The first line gives the list of energies (in GeV), the first column gives the list of times (in s), and each further column is the emission rate of the particle per unit energy, time and covolume (in $GeV^{-1}s^{-1}cm^{-3}$).

- `*_secondary_spectrum.txt`: these files are also written by `write_lines`. They contain the emission rates of each secondary particles at each final times and for each simulated final energies. The first line gives the list of energies (in GeV), the first column gives the list of times (in s), and each other column is the emission
rate of the particle per units of energy, time and covolume (in $GeV^{-1}s^{-1}cm^{-3}$). These files will not be generated if the parameter primary_only has been set to 1.

## 7.2 BlackHawk_inst

Running `BlackHawk_inst` produces 3 (or 2) output files:

- `BH_spectrum.txt`: this file is written by the routine `write_spectrum`. It contains the initial density spectrum of BHs, and has 3 columns: the first one is a list of BHs initial masses (in g), the second one the corresponding list of initial spins (dimensionless) and the third one is the comoving number densities (in $cm^{-3}$).

- `instantaneous_primary_spectra.txt`: this file is written by `write_instantaneous_primary_spectra`. It contains the emission rates of the primary particles for each simulated initial energy. The first line is the list of primary particles, the first column is the list of energies (in GeV), and each other column is the emission rate per unit energy and time (in $GeV^{-1}s^{-1}cm^{-3}$).

- `instantaneous_secondary_spectra.txt`: this file is written by `write_instantaneousHadronized_spectra`. It contains the emission rates of the secondary particles for each simulated final energy. The first line is the list of secondary particles, the first column is that of energies, and each other column is the emission rate per unit energy and time (in $GeV^{-1}s^{-1}cm^{-3}$). It will not be generated if the parameter primary_only has been set to 1.

# 8 Memory use

The code `BlackHawk` has been designed to minimize the memory used (both RAM and disk) and the computation time while avoiding excessive approximations. In this Section we give estimates of the memory used by each program.

## 8.1 RAM used

To every array defined in `BlackHawk`, a memory space is allocated with a `malloc` call. This memory is freed at the moment the array stops being necessary for the following part of the run. Then, the RAM used by `BlackHawk` at a given step of a session (corresponding to a paragraph in Sect. 6) can be estimated as a sum over all active arrays at that time. `double` are coded in 8 bytes and `int` in 4 bytes. Memory spaces $M$ are given in bytes. For `BlackHawk_tot` we have:

- step 1 (BH spectrum):

- `init masses [] = 8 x BHnumber`

- `init spins [] = 8 x BHnumber`

- `spec table [] = 8 x BHnumber`

- step 2 (BH evolution):

- `init masses [] = 8 x BHnumber`

- `init spins [] = 8 x BHnumber`

- `spec table [] = 8 x BHnumber`

- `fM_table [] [] = 8 x nb_fM_a x nb_fM_masses`

- `gM_table [] [] = 8 x nb_fM_a x nb_fM_masses`

- `fM masses [] = 8 x nb_fM_masses`

- `fM_a [] = 8 x nb_fM_a`

- `life masses [] [] = 8 x BHnumber^2 x limit`

- `life spins [] [] = 8 x BHnumber^2 x limit`

- `life times [] = 8 x BHnumber x limit`

- `dts [] = 8 x BHnumber x limit`

- `evolution length [] = 4 x BHnumber`

- step 3 (primary and secondary spectra):

- `spec table [] = 8 x BHnumber`

- `life masses [] [] = 8 x BHnumber^2 x limit`

- `life spins [] [] = 8 x BHnumber^2 x limit`

- `life times [] = 8 x BHnumber x limit`

- `dts [] = 8 x BHnumber x limit`

- `evolution length [] = 4 x BHnumber`

- `gammas [] [] [] = 8 x 4 x nb{gamma_a} x nb{gamma_x}`

- `gamma_a [] = 8 x nb{gamma_a}`

- `gamma_x [] = 8 x nb{gamma_x}`

- `fits [] [] [] = 8 x 4 x nb{gamma_a} x 7`

- `dof [] = 8 x (particle_number + grav)`

- `spins [] = 8 x (particle_number + grav)`

- `masses_primary [] = 8 x (particle_number + grav)`

- `times [] ≈ 8 x limit x BHnumber`

- `energies [] = 8 x Enumer`

- `tables [] [] [] [] = 8 x nb_fin_part x nb_init_en x nb_fin_en x nb_fin_part`

- `initial energies [] = 8 x nb_init_en`

- `final energies [] = 8 x nb_fin_en`

- `partial_hadronized_spectra [] [] [] = 8 x nb_fin_part x Enumer x nb_fin_en`

- `partial_primary_spectra [] [] = 8 x (particle_number + grav) x Enumer`

- `partial_integrated_hadronized_spectra [] [] = 8 x nb_fin_part x nb_fin_en`

- `masses_secondary [] = 8 x nb_fin_part`
Using the parameters of Appendix D.1, the arrays occupy at most ~ 150 MB. For `BlackHawk_inst` we have:

- step 1 (BH spectrum):

- BH masses [] = 8 x BHnumber

- BH spins [] = 8 x BHnumber

- spec table [] = 8 x BHnumber

- step 2 (primary spectra):

- BH masses [] = 8 x BHnumber

- BH spins [] = 8 x BHnumber

- spec table [] = 8 x BHnumber

- gammas [] [] [] = 8 x 4 x nb Gamma_a x nb Gamma_x

- gamma_a [] = 8 x nb Gamma_a

- gamma_x [] = 8 x nb Gamma_x

- fits [] [] [] = 8 x 4 x nb Gamma_a x 7

- dof [] = 8 x (particle_number + grav)

- spins [] = 8 x (particle_number + grav)

- masses primary [] = 8 x (particle_number + grav)

- instantaneous primary spectra [] [] = 8 x (particle_number + grav) x Enumer

- energies [] = 8 x Enumer

- step 3 (during hadronization):

- instantaneous primary spectra [] [] = 8 x (particle_number + grav) x Enumer

- energies [] = 8 x Enumer

- tables [] [] [] [] = 8 x nb Fin_part x nb Init_en x nb Fin_en x nb Fin_part

- initial energies [] = 8 x nb Init_en

- final energies [] = 8 x nb Fin_en

- masses secondary [] = 8 x nb Fin_part

- instantaneous hadronized spectra [] [] = 8 x nb Fin_part x Enumer x nb Fin_en

- step 3 bis (during integration):

- instantaneous primary spectra [] [] = 8 x (particle_number + grav) x Enumer

- energies [] = 8 x Enumer

- initial energies [] = 8 x nb Init_en

- final energies [] = 8 x nb Fin_en

- instantaneous hadronized spectra [] [] = 8 x nb Fin_times x Enumer x nb Fin_en

- instantaneous integrated hadronized spectra [] [] = 8 x nb Fin_part x nb Fin_en

Using the parameters of Appendix D.1, the arrays occupy at most ~ 10 MB.

## 8.2 Static disk memory used

The output generated is written in `.txt` files using a precision of 5 significant digits. Adding the exponent and the coma, we get to 12 characters per written number, which is 12 bytes. For `BlackHawk_tot` we have:

- file BH spectrum.txt: $M = 12 \times 3 \times BHnumber$

- file life evolutions.txt: $M \approx 4 \times 3 \times BHnumber + 12 \times 2 \times BHnumber^2 \times 1000$ where an average number of 1000 iterations for the mass integration of BHs has been assumed.

- files * primary spectrum.txt: $M = 12 \times (particle_number + grav) \times Enumer \times 1000 \times BHnumber$ where an average number of 1000 iterations for the mass integration of BHs has been assumed.

- files * secondary spectrum.txt: $M = 12 \times nb Fin_part \times nb Fin_en \times 1000 \times BHnumber$ where an average number of 1000 iterations for the mass integration of BHs has been assumed.

Using the parameters of Appendix D.1, the total written disk space is ~ 230 MB. For `BlackHawk_inst` we have:

- file BH spectrum.txt: $M = 12 \times 3 \times BHnumber$

- file instantaneous primary spectra.txt: $M = 12 \times Enumer \times (particle_number + grav)$

- file instantaneous secondary spectra.txt: $M = 12 \times nb Fin_en \times nb Fin_part$

Using the parameters of Appendix D.1, the total written disk space is ~ 35 kB.

# 9 Other applications

In this Section we present some hints on how to modify `BlackHawk`. Most of these modifications will require additions in the file `parameters.txt` and thus a modification of the routine `read_params` and of the structure `struct param`.

## 9.1 Computing new numerical tables

The user may be interested in recomputing the tables described in Appendix B, either to have more entries or to compute them with different methods for comparison. The easiest way to add tables in `BlackHawk` would be:

- authorize the corresponding "choice" parameters to have other integer values;

- put the new tables in a new directory in the `tables/` sub-folder;
- add a switch into the tables reading routines;

- make sure that the way tables are used in the routines will be compatible with the format of the new ones.

All the scripts used to compute the current tables are included in *BlackHawk* in the sub-folder `scripts/` together with `README.txt` files.

## 9.2 Using another Black Hole mass function

The user may be interested in testing its own BH distribution. Here are the main steps to add a pre-built distribution:

- add a “choice” parameter to the `struct param` choosing the distribution,

- add the corresponding analytical formula to the routine `n_cov` or tabulated values in the sub-folder `tables/`,

- modify the parameter `tmin` if the distribution is valid at a different initial time.

Providing a tabulated initial distribution to *BlackHawk* is done by switching the parameter `spectrumansto` to 5, putting the table file in the sub-folder `users_spectra/` and giving its full file name (including the extension) to the parameter `table`. The format has to be:

- three same-length columns, the first one for BHs masses $M$, the second one for BHs spins $a^*$ and the third one for the comoving number densities $dn(M)$ (with $dM$ taken around $M$),

- masses and densities in CGS units ($g$ and $cm^{-3}$ respectively), spins in dimensionless form,

- numbers in standard scientific notation,

- no additional text.

## 9.3 Adding primary particles

If the user wants to add hypothetical primary Hawking particles, the following steps have to be undertaken:

- enhance the parameter `particle_number` or add the new particle(s) with a switch similar to the one of the graviton,

- recompute the $f(M, a^*)$ and $g(M, a^*)$ tables to account for this(ese) new emission(s),

- if the spin(s) of the new particle(s) is(are) not among the greybody factor tables, compute the new ones,

- add the new particle(s) to all the fixed length arrays of particle types (e.g. the file names or columns in the writing routines),

- eventually add its(their) contribution(s) to the secondary spectra.

## 9.4 Adding secondary particles

In order to add secondary Hawking particles to the code, one has to:

- recompute the hadronization tables to take new branching ratios into account,

- add the new particle(s) to all the fixed length arrays of particle types (e.g. the file names or columns in the writing routines),

- add the corresponding contribution(s) to the routine `contribution instantanous`.

## 9.5 Other types of Black Holes

If the user wants to compute the Hawking emission of BHs different from the Schwarzschild or Kerr ones, several ingredients are needed:

- add a switch to the parameter file to select amongst the new types of BHs,

- modify/add the Hawking temperature function `temp_BH` for these BHs,

- modify/add evolution routines `loss_rate_*` and `life_evolution` (e.g. for charged BHs a routine `loss_rate_Q` for the evolution of the charge parameter $Q$),

- compute the corresponding $f, g$ and eventually new evolution parameters tables and add the corresponding reading routines (e.g. for charged BHs a routine `read_hM_table` to read the $q(M, a^*, Q)$ table where $q$ would describe the evolution of the electric charge $Q$),

- compute the new greybody factors tables and update the corresponding reading and interpolating routines `read, read_gammas, read_gammas_fits` and `dNtdE`.

Depending on the complexity of the BH model, the user may need to implement some or all of the above modifications.

# 10 Conclusion

*BlackHawk* is the first public code generating both primary and secondary Hawking radiation spectra for any mass distribution of Schwarzschild and Kerr Black Holes, and their evolution in time. The primary spectra are obtained using greybody factors, and the secondary ones result from the decay and hadronization of the primary particles. The Black Hole and spectrum evolutions are obtained by considering the energy loss via Hawking radiation and the subsequent modification of the temperature of the Black Hole. *BlackHawk* is
designed in a user-friendly way and modifications can be easily implemented. The prime application is to study the effects of particles generated by Hawking radiation on observable quantities and thus to disqualify or set constraints on cosmological models involving the formation of Black Holes, as well as to test the Hawking radiation assumptions and study Black Hole general properties.

**Acknowledgements** We gratefully acknowledge helpful exchanges with P. Richardson in particular on the hadronization procedure and the HERWIG code. We are also thankful to J. Silk for many constructive discussions, to P. Skands for help with PYTHIA and hadronization, and to G. Robbins for the interface with SuperIso Relic. The authors thank the CERN theory group for its hospitality during which part of this work was done.

**Data Availability Statement** This manuscript has associated data in a data repository. [Authors' comment: The code can be downloaded from https://blackhawk.hepforge.org.]

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP³.

## Appendix A: Units

The BlackHawk code uses the GeV unit internally in order to have simpler analytical expressions. However, to make the user interface more accessible, the input parameters as well as the output files are in CGS units. We provide below unit conversions from the natural system of units where $\hbar = c = k_B = G = 1$ to CGS or SI.

### Appendix A.1: Energy

The energy conversion from GeV to Joule is:

$$E_J = 1.602176565 \times 10^{-10} E_{\text{GeV}}. \qquad (A.1)$$

### Appendix A.2: Mass

The dimensional link between energy and mass is $[m] = [E/c^2]$, and the conversion from GeV to grams is:

$$m_g = 5.60958884 \times 10^{23} m_{\text{GeV}}. \qquad (A.2)$$

### Appendix A.3: Time

The dimensional link between energy and time is $[t] = [\hbar/E]$, and the conversion from GeV to seconds is:

$$t_s = 1.519267407 \times 10^{24} t_{\text{GeV}^{-1}}. \qquad (A.3)$$

## Appendix A.4: Distance

The dimensional link between energy and distance is $[l] = [\hbar c/E]$, and the conversion from GeV to meters is:

$$l_{\text{cm}} = 5.06773058 \times 10^{13} l_{\text{GeV}^{-1}}. \qquad (A.4)$$

## Appendix A.5: Temperature

The dimensional link between energy and temperature is $[T] = [E/k_B]$, and the conversion from GeV to Kelvins is:

$$T_K = 8.61733063 \times 10^{-14} T_{\text{GeV}}. \qquad (A.5)$$

## Appendix B: Computation of the tables

### Appendix B.1: Greybody factors

Chandrasekhar and Detweiler have shown that the Teukolsky equation can be reduced to a wave equation for Kerr Black Holes [51–54]. It is indeed difficult to find short-range potentials allowing for precise numerical computation. They give the form of such potentials in [51,52] for spin 2, [53] for spins 0 and 1 and [54] for spin 1/2, and find necessary to define a modified Eddington-Finkelstein radial coordinate $r^*$ by

$$\frac{dr^*}{dr} = \frac{\rho^2}{\Delta}, \qquad (B.6)$$

where $\rho(r)^2 \equiv r^2 + \alpha^2$ and $\alpha^2 \equiv a^2 + am/E$, $a$ being the BH spin and $m$ the projection of the angular momentum $l$. This equation can be integrated to give

$$r^*(r) = r + \frac{r_H r_+ + am/E}{r_+ - r_-} \ln \left( \frac{r}{r_+} - 1 \right) - \frac{r_H r_- + am/E}{r_+ - r_-} \ln \left( \frac{r}{r_-} - 1 \right). \qquad (B.7)$$

Unfortunately, the inverse of this equation has to be found numerically and is generally difficult to determine with accurate precision.

As boundary conditions for the wave equation, we use a purely outgoing wave. The solution at the horizon has the form

$$\psi_s = e^{-iEr^*}. \qquad (B.8)$$

At infinity, the solution has the form

$$\psi_s = A_{\text{in}}e^{iEr^*} + A_{\text{out}}e^{-iEr^*}, \qquad (B.9)$$

The Schrödinger-like wave equation is for all spins

$$\frac{d^2\psi_s}{dr^{*2}} + (E^2 - V_s(r^*))\psi_s = 0. \qquad (B.10)$$
The method to transform Eq. (15) into this simple wave equation was proposed in the Chandrasekhar & Detweiler papers [51–54]. The potentials are<fnref n="8" />

$$V_0(r) = \frac{\Delta}{\rho^4} \left( \lambda_{0lm} + \frac{\Delta + 2r(r-M)}{\rho^2} - \frac{3r^2\Delta}{\rho^4} \right), \quad (B.11)$$

$$V_{1/2,\pm}(r) = (\lambda_{1/2lm} + 1) \frac{\Delta}{\rho^4} \mp \frac{\sqrt{(\lambda_{1/2,l,m} + 1)\Delta}}{\rho^4} \times \left( (r-M) - \frac{2r\Delta}{\rho^2} \right), \quad (B.12)$$

$$V_{1,\pm}(r) = \frac{\Delta}{\rho^4} \left( (\lambda_{1lm} + 2) - \alpha^2 \frac{\Delta}{\rho^4} \mp i\alpha\rho^2 \frac{d}{dr} \left( \frac{\Delta}{\rho^4} \right) \right), \quad (B.13)$$

$$V_2(r) = \frac{\Delta}{\rho^8} \left( q - \frac{\rho^2}{(q - \beta\Delta)^2} \left( (q - \beta\Delta)(\rho^2\Delta q'' - 2\rho^2q - 2r(q'\Delta - q\Delta') + \rho^2(\kappa\rho^2 - q' + \beta\Delta')) \times (q'\Delta - q\Delta') \right) \right). \quad (B.14)$$

The different potentials for a given spin lead to the same results. In the potential for spin 2 particles, the following quantities appear

$$q(r) = \nu\rho^4 + 3\rho^2(r^2 - a^2) - 3r^2\Delta, \quad (B.15)$$

$$q'(r) = r \left( (4\nu + 6)\rho^2 - 6(r^2 - 3Mr + 2a^2) \right), \quad (B.16)$$

$$q''(r) = (4\nu + 6)\rho^2 + 8\nu r^2 - 6r^2 + 36Mr - 12a^2, \quad (B.17)$$

$$q'\Delta - q\Delta' = -2(r - M)\nu\rho^4 + 2\rho^2(2\nu r\Delta - 3M(r^2 + a^2) + 6ra^2) + 12r\Delta(Mr - a^2), \quad (B.18)$$

$$\beta_{\pm} = \pm 3\alpha^2, \quad (B.19)$$

$$\kappa_{\pm} = \pm \sqrt{36M^2 - 2\nu(\alpha^2(5\nu + 6) - 12a^2) + 2\beta\nu(\nu + 2)}, \quad (B.20)$$

$$q - \beta_{+}\Delta = \rho^2(v\rho^2 + 6Mr - 6a^2), \quad (B.21)$$

$$q - \beta_{-}\Delta = \nu\rho^4 + 6r^2(\alpha^2 - a^2) + 6Mr(r^2 - \alpha^2), \quad (B.22)$$

where $\nu \equiv \lambda_{2lm} + 4$.

In the Schwarzschild limit ($a = 0$), we recover the Regge-Wheeler potentials. The angular momentum projection $m$

only appears multiplied by $a$, which simplifies the calculation since only one common value for all $m$ has to be chosen once $l$ is fixed.

The $r^*$ variable change used in these potentials leads to divergences in the potentials, when $r_{\text{div}}^2 = -\alpha^2$. This can happen for sufficiently low energies and high (negative) angular momentum projections, and it corresponds to the *superradiance* regime. As discussed in the Chandrasekhar-Detweiler papers, the technique to avoid this divergence is to integrate Eq. (B.10) up to slightly before the divergence (e.g. $r_{\text{div}} - \epsilon$). At this point, the asymptotic behaviour of the potential $V_s$ is known, and Eq. (B.10) is simplified. Since the form of the function $\psi_s$ can be obtained, by continuity of the function $R_s$ of Eq. (15) one can extrapolate this form up to slightly *after* the divergence (e.g. $r_{\text{div}} + \epsilon$) and continue the integration.

Another difficulty which can arise is the fact that there can be an additional divergence in the spin 2 potential because of the $q - \beta_{\pm}\Delta$ term. For this extra divergence, we try to integrate the wave equation with one of the potentials (e.g. $\kappa_+$, $\beta_+$), and in case of problem we try with the other potentials (e.g. $\kappa_+$, $\beta_-$), as it seems that at least one of the four potentials does not generate any divergence.

The greybody factor is given by the transmission coefficient of the wave from the horizon to space infinity

$$\Gamma_{slm} \equiv T_{slm} = \frac{1}{|A_{\text{out}}|^2}. \quad (B.23)$$

Practically, we compute the value of the single dof emissivities

$$Q_s \equiv \sum_{l,m} \frac{\Gamma_{slm}}{(e^{E/T} \pm 1)}, \quad (B.24)$$

for some values of $a^*$ and for a range of $0.01 < x \equiv 2Er_{\text{BH}} < 5$ (dimensionless), since we can show that these are the only relevant parameters for massless particles. For $x$ out of this range, we have found easier to find empiric asymptotic forms of the emissivities. At low energies, we have for all spins

$$\log_{10}(Q_s) \approx a_{1,s} \log_{10}(x) + a_{2,s}, \quad (B.25)$$

and at high energies

$$\log_{10}(Q_s) \approx a_{3,s}x + a_{4,s} + a_{5,s}\cos(a_{7,s}x) + a_{6,s}\sin(a_{7,s}x). \quad (B.26)$$

We fitted the computed emissivities to find the values of the parameters $a_{i,s}$. We checked that they agree with the asymptotic limits of [20] in the Schwarzschild case and of [43] in the Kerr case.

The Mathematica scripts `spin_*` .m, the fitting script `exploitation.m` as well as a C formatting script `formatting.c` and a `README.txt` are provided in the sub-folder:

<footnote n="8">⁸ We found that the spin 0 potential had a missing “r” in [53].</footnote>
scripts/greybody scripts/greybody_factors/

Please contact one of the authors if you have issues using
these scripts.

## Appendix B.2: Evolution tables

To compute the integrals of Eqs. (21) and (22), we use the greybody factor tables and the fits computed in Appendix B.1. The peak of Hawking emission lies around the BH temperature (see [20] for example), thus the integral does not need to be computed over all energies, but a restrained set $10^{-5} \times T < E < 10^{5} \times T$ is sufficient. The domains of integration are segmented over logarithmically distributed energies, and computed for masses between $M_P$ to $10^{46}$ GeV (~$10^{-5} - 10^{22}$ g). In the limit where neutrinos are massless, $f(M, a^*)$ and $g(M, a^*)$ are not expected to change for masses higher than $10^{22}$ g, the tables can therefore be extended manually without any new computation. Masses are given in GeV (corresponding to grams) and $f(M, a^*)$ and $g(M, a^*)$ are in GeV$^4$ (corresponding to $g^3 \cdot s^{-1}$ and $g^2 \cdot \text{GeV} \cdot s^{-1}$, respectively). We have checked that the value of $f(M, 0)$ is consistent with that of [41] in the Schwarzschild case and that the values of $f(M, a^*)$ and $g(M, a^*)$ are consistent with [55] in the Kerr case.

The C script `fM.c` used to compute the tables and a `README.txt` are provided in the subfolder:

scripts/greybody scripts/fM/

Please contact one of the authors if you have issues using
these scripts.

## Appendix B.3: Hadronization

Two particle physics codes have been used to compute hadronization tables: HERWIG [45] and PYTHIA [46]. In both cases, the strategy is to generate the output of a collision (for example $e^+e^- \to u+\bar{u} \to ...$), and then to count the number of final particles (here denoted as dots) normalized by the number (here 2) of initial particles (here $u$, see Table 2 in Appendix C) satisfying the desired stability criterion: Table 3 for Early Universe/BBN particles (PYTHIA and HERWIG tables) and Table 4 for present epoch particles

(PYTHIA “new” tables). This gives the number of secondary Hawking particles of each type that a primary particle will generate.

To build the PYTHIA and HERWIG tables, we have simulated for each channel listed in Table 1, $10^5$ events for initial energies $E'$ (half of the center of mass energy) logarithmically distributed between 5 GeV and $10^5$ GeV (PYTHIA and PYTHIA “new”) or between 25 GeV and $10^5$ GeV (HERWIG). Then, the final particles have been listed as a function of their final energy $E$, into a range of $10^{-6}$ GeV to $10^5$ GeV and the counts have been averaged over the number or simulated events. This gives the dimensionless quantity $dN_j^i(E', E)$ of Eq. (24). Energies are given in GeV.

The branching ratios $e^\pm \to \gamma\gamma \to \dots$ and $e^\pm \to \nu\bar{\nu} \to \dots$ have not been computed. The contribution from the primary photons and neutrinos is directly added to the secondary spectra with a branching ratio of 1.

fOR initial energies lower than the cutoff of the computed tables, branching ratios from the lowest relevant initial energy will be extrapolated at lower energies once shifted to the considered energy, taking into account that no emission can arise below the rest mass of the final particles. The same kind of extrapolation is used with high-energy hadronization. There is however no guarantee that the extrapolations remain valid far beyond the cutoff energies.

The PYTHIA (new) and HERWIG scripts used to run the particle physics codes, as well as the C scripts `formatting.c` used to format the hadronization tables and `README.txt` files are provided in the sub-folders:

scripts/pythia scripts/

scripts/herwig scripts/

scripts/pythia scripts_new/

Please contact one of the authors if you have issues using
these scripts.

## Appendix C: Particle information

In this Appendix we give the particle properties taken from [15] used in BlackHawk and summed in Tables 2, 3 and 4.

<table><caption>**Table 1** List of the channels used to compute the hadronization tables</caption><thead><tr><th>Particle</th><th>PYTHIA (new)</th><th>HERWIG</th></tr></thead><tbody><tr><td>Gluons</td><td>e<sup>+</sup>e<sup>-</sup> → h<sup>0</sup> → g<span style="text-decoration: overline">g</span></td><td>e<sup>+</sup>e<sup>-</sup> → h<sup>0</sup> → g<span style="text-decoration: overline">g</span></td></tr><tr><td>Higgs boson</td><td>e<sup>+</sup>e<sup>-</sup> → h<sup>0</sup></td><td>e<sup>+</sup>e<sup>-</sup> → h<sup>0</sup></td></tr><tr><td>W bosons</td><td>e<sup>+</sup>e<sup>-</sup> → Z<sup>0</sup>/γ<sup>*</sup> → W<sup>+</sup>W<sup>-</sup></td><td>e<sup>+</sup>e<sup>-</sup> → Z<sup>0</sup>/γ<sup>*</sup> → W<sup>+</sup>W<sup>-</sup></td></tr><tr><td>Z boson</td><td>e<sup>+</sup>e<sup>-</sup> → h<sup>0</sup> → Z<sup>0</sup>Z<sup>0</sup></td><td>e<sup>+</sup>e<sup>-</sup> → Z<sup>0</sup>/γ<sup>*</sup> → Z<sup>0</sup>Z<sup>0</sup></td></tr><tr><td>Leptons</td><td>e<sup>+</sup>e<sup>-</sup> → h<sup>0</sup> → l<sup>+</sup>l<sup>-</sup></td><td>e<sup>+</sup>e<sup>-</sup> → Z<sup>0</sup>/γ<sup>*</sup> → l<sup>+</sup>l<sup>-</sup></td></tr><tr><td>Quarks</td><td>e<sup>+</sup>e<sup>-</sup> → Z<sup>0</sup>/γ<sup>*</sup> → q<span style="text-decoration: overline">q</span></td><td>e<sup>+</sup>e<sup>-</sup> → Z<sup>0</sup>/γ<sup>*</sup> → q<span style="text-decoration: overline">q</span></td></tr></tbody></table>
<table><caption>**Table 2** Properties of the elementary particles of the Standard Model, in addition to the graviton [15]. The number of quantum dof is the product of the family, antiparticle, the colour and the helicity multiplicities. Neutrinos are here considered massless. In the code, gluons have been assigned an effective mass to account for the QCD energy scale $\Lambda \approx 200$ MeV</caption><thead><tr><th>Particle</th><th>Symbol</th><th>Mass (GeV/c<sup>2</sup>)</th><th>Spin</th><th>Quantum dof</th></tr></thead><tbody><tr><td>Higgs boson</td><td>h<sup>0</sup></td><td>1.2503 &times; 10<sup>2</sup></td><td>0</td><td>1</td></tr><tr><td>Photon</td><td>&gamma;</td><td>0</td><td>1</td><td>2</td></tr><tr><td>Gluons</td><td>g</td><td>0</td><td>1</td><td>16</td></tr><tr><td>W bosons</td><td>W<sup>&pm;</sup></td><td>8.0403 &times; 10<sup>1</sup></td><td>1</td><td>6</td></tr><tr><td>Z boson</td><td>Z<sup>0</sup></td><td>9.11876 &times; 10<sup>1</sup></td><td>1</td><td>3</td></tr><tr><td>Neutrinos</td><td>&nu;<sub>e,&mu;&tau;</sub>, &overline{&nu;<sub>e,&mu;&tau;</sub>}</td><td>0</td><td>1/2</td><td>6</td></tr><tr><td>Electron</td><td>e<sup>&pm;</sup></td><td>5.109989461 &times; 10<sup>-4</sup></td><td>1/2</td><td>4</td></tr><tr><td>Muon</td><td>&mu;<sup>&pm;</sup></td><td>1.056583745 &times; 10<sup>-1</sup></td><td>1/2</td><td>4</td></tr><tr><td>Tau</td><td>&tau;<sup>&pm;</sup></td><td>1.77686</td><td>1/2</td><td>4</td></tr><tr><td>Up quark</td><td>u, &overline{u}</td><td>2.2 &times; 10<sup>-3</sup></td><td>1/2</td><td>12</td></tr><tr><td>Down quark</td><td>d, &overline{d}</td><td>4.7 &times; 10<sup>-3</sup></td><td>1/2</td><td>12</td></tr><tr><td>Charm quark</td><td>c, &overline{c}</td><td>1.27</td><td>1/2</td><td>12</td></tr><tr><td>Strange quark</td><td>s, &overline{s}</td><td>9.6 &times; 10<sup>-2</sup></td><td>1/2</td><td>12</td></tr><tr><td>Top quark</td><td>t, &overline{t}</td><td>1.7321 &times; 10<sup>2</sup></td><td>1/2</td><td>12</td></tr><tr><td>Bottom quark</td><td>b, &overline{b}</td><td>4.18</td><td>1/2</td><td>12</td></tr><tr><td>Graviton</td><td>G</td><td>0</td><td>2</td><td>2</td></tr></tbody></table>

<table><caption>**Table 3** Particles with a lifetime longer than $10^{-8}$ s, relevant for Early Universe/BBN studies [15,47] and used to compute the hadronization tables in `pythia_tables/` and `herwig_tables/`</caption><thead><tr><th>Particle</th><th>Symbol</th><th>Lifetime (s)</th></tr></thead><tbody><tr><td>Photon</td><td>&gamma;</td><td>&infin;</td></tr><tr><td>Electron</td><td>e<sup>&pm;</sup></td><td>&infin;</td></tr><tr><td>Muon</td><td>&mu;<sup>&pm;</sup></td><td>(2.1969811 &pm; 0.0000022) &times; 10<sup>-6</sup></td></tr><tr><td>Neutrinos</td><td>&nu;<sub>e,&mu;&tau;</sub>, &overline{&nu;<sub>e,&mu;&tau;</sub>}</td><td>&infin;</td></tr><tr><td>Charged pions</td><td>&pi;<sup>&pm;</sup></td><td>(2.6033 &pm; 0.0005) &times; 10<sup>-8</sup></td></tr><tr><td>Neutral “long” kaon</td><td>K<sup>0</sup><sub>L</sub></td><td>(5.099 &pm; 0.021) &times; 10<sup>-8</sup></td></tr><tr><td>Charged kaons</td><td>K<sup>&pm;</sup></td><td>(1.2379 &pm; 0.0021) &times; 10<sup>-8</sup></td></tr><tr><td>Proton</td><td>p, &overline{p}</td><td>&infin;</td></tr><tr><td>Neutron</td><td>n, &overline{n}</td><td>880.2 &pm; 1</td></tr></tbody></table>

<table><caption>**Table 4** Stable particles, relevant for evaporating BH in the present Universe and used to compute the hadronization table in `pythia_tables_new/`</caption><thead><tr><th>Particle</th><th>Symbol</th><th>Lifetime (s)</th></tr></thead><tbody><tr><td>Photon</td><td>&gamma;</td><td>&infin;</td></tr><tr><td>Electron</td><td>e<sup>&pm;</sup></td><td>&infin;</td></tr><tr><td>Neutrinos</td><td>&nu;<sub>e,&mu;&tau;</sub>, &overline{&nu;<sub>e,&mu;&tau;</sub>}</td><td>&infin;</td></tr><tr><td>Proton</td><td>p, &overline{p}</td><td>&infin;</td></tr></tbody></table>

## Appendix D: Results

The results in the output files are given in CGS units.

### Appendix D.1: Parameters

An example of `parameters.txt` file is given here:
```
destination folder = test
full_output = 0
```

`interpolation_method = 1`

BHnumber = 10

Mmin = 1.e+15

Mmax = 1.e+17

a = 0.

`spectrum.choice = 1`

`amplitude_lognormal = 1.e-30`

`variance_lognormal = 1.`

`critPhiMass_lognormal = 1.e+10`

`amplitude_powerlaw = 1.`

`eqStatePowerlaw = 0.3333`

`amplitude_critical_collapse = 1.`

`critPhiMass_critical_collapse = 1.e+10`
```
table = table.txt

tmin = 1.e-5
limit = 5000

Mmin_fM = 2.177e-5
Mmax_fM = 1.783e+22
amin_fM = 0.
amax_fM = 0.9999
nb_fM大人 = 100
nb_fM_a = 50

Enumber = 100
Emin = 1.e-6
Emax = 100000.
particle_number = 15
grav = 1

nb_a = 50
nbgamma_x = 200

primary_only = 0

hadronizationChoice = 2

Emin_hadro_PYTHIA = 5.
Emax_hadro_PYTHIA = 100000.
nb_init_en_PYTHIA = 250
nb_fi_n_en_PYTHIA = 500
nb_init_part_PYTHIA = 14
nb_fi_n_part_PYTHIA = 11

Emin_hadro_HERWIG = 25.
Emax_hadro_HERWIG = 100000.
nb_init_en_HERWIG = 100
nb_fi_n_en_HERWIG = 100
nb_init_part_HERWIG = 14
nb_fi_n_part_HERWIG = 11
```

```
Emin_hadro_PYTHIA_new = 5.
Emax_hadro_PYTHIA_new = 100000.
nb_init_en_PYTHIA_new = 250
nb_fi_n_en_PYTHIA_new = 500
nb_init_part_PYTHIA_new = 14
nb_fi_n_part_PYTHIA_new = 6
```

## Appendix D.2: BlackHawk_tot

When running BlackHawk_tot with the parameters of Appendix D.1, the output file BH_spectrum.txt is generated in the folder test/ and should look like Table 5.

The output file life_evolutions.txt should look like Table 6.

The output file photon_primary_spectrum.txt should look like Table 7.

The output file photon_secondary_spectrum.txt should look like Table 8.

<table><caption>**Table 5** Initial BH comoving number density as a function of their mass</caption><thead><tr><td>Mass</td><td>Spin</td><td>Density</td></tr></thead><tbody><tr><td>1.00000e+015</td><td>0.00000e+000</td><td>4.43371e-034</td></tr><tr><td>1.66810e+015</td><td>0.00000e+000</td><td>6.44560e-037</td></tr><tr><td>2.78256e+015</td><td>0.00000e+000</td><td>7.21193e-040</td></tr><tr><td>4.64159e+015</td><td>0.00000e+000</td><td>6.21057e-043</td></tr><tr><td>7.74264e+015</td><td>0.00000e+000</td><td>4.11627e-046</td></tr><tr><td>1.29155e+016</td><td>0.00000e+000</td><td>2.09975e-049</td></tr><tr><td>2.15443e+016</td><td>0.00000e+000</td><td>8.24375e-053</td></tr><tr><td>3.59381e+016</td><td>0.00000e+000</td><td>2.49100e-056</td></tr><tr><td>5.99484e+016</td><td>0.00000e+000</td><td>5.79314e-060</td></tr><tr><td>1.00000e+017</td><td>0.00000e+000</td><td>1.03692e-063</td></tr></tbody></table>
<table><caption>**Table 6** Evolution of the BH masses and spins as functions of time</caption><thead><tr><th colspan="10">Number of timesteps for each BH initial mass</th></tr><tr><th>838</th><th>1596</th><th>2449</th><th>3181</th><th>4006</th><th>4866</th><th>5779</th><th>6639</th><th>7574</th><th>8395</th></tr><tr><th>Time</th><th colspan="9">Masses</th></tr></thead><tbody><tr><td>1.00000e-005</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>2.00000e-005</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>4.00000e-005</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>8.00000e-005</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1.60000e-004</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>3.20000e-004</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>6.40000e-004</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1.28000e-003</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>2.56000e-003</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>5.12000e-003</td><td>1.00000e+015</td><td>1.66810e+015</td><td>2.78256e+015</td><td>4.64159e+015</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><th>Time</th><th colspan="9">Spins</th></tr><tr><td>1.00000e-005</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>2.00000e-005</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

<table><caption>**Table 7** Hawking primary spectrum as a function of time</caption><thead><tr><th>Time/energy</th><th>1.00000e-006</th><th>1.29155e-006</th><th>1.66810e-006</th><th>2.15443e-006</th><th>...</th></tr></thead><tbody><tr><td>1.00000e-005</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>2.00000e-005</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>4.00000e-005</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>8.00000e-005</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>1.60000e-004</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>3.20000e-004</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>6.40000e-004</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>1.28000e-003</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>2.56000e-003</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>5.12000e-003</td><td>1.65033e-026</td><td>3.48856e-026</td><td>7.37429e-026</td><td>1.55882e-025</td><td>...</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>
<table><caption>**Table 8** Hawking secondary spectrum as a function of time</caption><thead><tr><th>Time/energy</th><th>1.00000e-006</th><th>1.05200e-006</th><th>1.10700e-006</th><th>1.16400e-006</th><th>...</th></tr></thead><tbody><tr><td>1.00000e-005</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>2.00000e-005</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>4.00000e-005</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>8.00000e-005</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>1.60000e-004</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>3.20000e-004</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>6.40000e-004</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>1.28000e-003</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>2.56000e-003</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>5.12000e-003</td><td>2.34734e-010</td><td>2.27114e-010</td><td>2.21976e-010</td><td>2.06182e-010</td><td>...</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

## Appendix D.3: BlackHawk_inst

When running `BlackHawk_inst` with the parameters of Appendix D.1, the output file `instantaneous_primary_spectra.txt` is generated in the folder `test/` and it should look like Table 9.

The output file `instantaneous_secondary_spectra.txt` should look like Table 10.

<table><caption>**Table 9** Hawking primary spectra for each particle types</caption><thead><tr><th>Energy/particle</th><th>Photon</th><th>Gluons</th><th>Higgs</th><th>W+-</th><th>...</th></tr></thead><tbody><tr><td>1.00000e-006</td><td>1.65033e-026</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>1.29155e-006</td><td>3.48856e-026</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>1.66810e-006</td><td>7.37429e-026</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>2.15443e-006</td><td>1.55882e-025</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>2.78256e-006</td><td>3.29511e-025</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>3.59381e-006</td><td>6.96536e-025</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>4.64159e-006</td><td>1.47237e-024</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>5.99484e-006</td><td>3.11238e-024</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>7.74264e-006</td><td>6.57911e-024</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>1.00000e-005</td><td>1.39073e-023</td><td>0.00000e+000</td><td>0.00000e+000</td><td>0.00000e+000</td><td>...</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

<figure><img src="image_3.png" /></figure>
<table><caption>**Table 10** Hawking secondary spectra for each particle type</caption><thead><tr><th>Energy/particle</th><th>Photon</th><th>Electron</th><th>nu_e</th><th>nu_mu</th><th>...</th></tr></thead><tbody><tr><td>1.00000e-006</td><td>2.34734e-010</td><td>0.00000e+000</td><td>7.60472e-011</td><td>6.42832e-015</td><td>...</td></tr><tr><td>1.05200e-006</td><td>2.27114e-010</td><td>0.00000e+000</td><td>7.84847e-011</td><td>5.29763e-015</td><td>...</td></tr><tr><td>1.10700e-006</td><td>2.21976e-010</td><td>0.00000e+000</td><td>8.64882e-011</td><td>5.65187e-015</td><td>...</td></tr><tr><td>1.16400e-006</td><td>2.06182e-010</td><td>0.00000e+000</td><td>9.49170e-011</td><td>6.81580e-015</td><td>...</td></tr><tr><td>1.22500e-006</td><td>2.05790e-010</td><td>0.00000e+000</td><td>8.94352e-011</td><td>1.15943e-014</td><td>...</td></tr><tr><td>1.28900e-006</td><td>1.84274e-010</td><td>0.00000e+000</td><td>9.97086e-011</td><td>1.42860e-014</td><td>...</td></tr><tr><td>1.35600e-006</td><td>1.78922e-010</td><td>0.00000e+000</td><td>9.95837e-011</td><td>1.27595e-014</td><td>...</td></tr><tr><td>1.42700e-006</td><td>1.73629e-010</td><td>0.00000e+000</td><td>1.08663e-010</td><td>1.34540e-014</td><td>...</td></tr><tr><td>1.50100e-006</td><td>1.64730e-010</td><td>0.00000e+000</td><td>1.17571e-010</td><td>1.59443e-014</td><td>...</td></tr><tr><td>1.57900e-006</td><td>1.59870e-010</td><td>0.00000e+000</td><td>1.12879e-010</td><td>2.51215e-014</td><td>...</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

References

1. B.P. Abbott et al., Phys. Rev. Lett. **116**(6), 061102 (2016). https://doi.org/10.1103/PhysRevLett.116.061102

2. B.P. Abbott et al., Astrophys. J. **848**(2), L12 (2017). https://doi.org/10.3847/2041-8213/aa91c9

3. LIGO Scientific, Virgo Collaboration, B.P. Abbott et al. GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs (2018). arXiv:1811.12907

4. K. Akiyama et al., Astrophys. J. **875**(1), L1 (2019). https://doi.org/10.3847/2041-8213/ab0ec7

5. B.J. Carr, K. Kohri, Y. Sendouda, J. Yokoyama, Phys. Rev. D **81**, 104019 (2010). https://doi.org/10.1103/PhysRevD.81.104019

6. B. Carr, F. Kuhnel, M. Sandstad, Phys. Rev. D **94**(8), 083504 (2016). https://doi.org/10.1103/PhysRevD.94.083504

7. L. Barack et al. Class. Quant. Grav. **36**(14), 143001 (2019). arXiv:1806.05195 [gr-qc]

8. A. Kashlinsky et al. Electromagnetic probes of primordial black holes as dark matter (2019). arXiv:1903.04424

9. D.C. Dai, C. Issever, E. Rizvi, G. Starkman, D. Stojkovic, J. Tseng. Comput. Phys. Commun. **236**, 285–301 (2019). arXiv:0902.3577 [hep-ph]

10. J.A. Frost, J.R. Gaunt, M.O.P. Sampaio, M. Casals, S.R. Dolan, M.A. Parker, B.R. Webber, J. High Energy Phys. **10**, 014 (2009). https://doi.org/10.1088/1126-6708/2009/10/014

11. B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, H. Veermäe, Phys. Rev. D **96**(2), 023514 (2017). https://doi.org/10.1103/PhysRevD.96.023514

12. A. Katz, J. Kopp, S. Sibiryakov, W. Xue, J. Cosmol. Astropart. Phys. **2018**(12), 005 (2018). https://doi.org/10.1088/1475-7516/2018/12/005

13. H. Tashiro, N. Sugiyama, Phys. Rev. D **78**(2), 023004 (2008). https://doi.org/10.1103/PhysRevD.78.023004

14. C. Germani, I. Musco, Phys. Rev. Lett. **122**(14), 141302 (2019). https://doi.org/10.1103/PhysRevLett.122.141302

15. M. Tanabashi et al., Phys. Rev. D **98**(3), 030001 (2018). https://doi.org/10.1103/PhysRevD.98.030001

16. S.W. Hawking, Commun. Math. Phys. **43**, 199 (1975). https://doi.org/10.1007/BF02345020

17. S.A. Teukolsky, Astrophys. Phys. J. **185**, 635 (1973). https://doi.org/10.1086/152444

18. S.A. Teukolsky, W.H. Press, Astrophys. Phys. J. **193**, 443 (1974). https://doi.org/10.1086/153180

19. C.M. Harris, P. Kanti, J. High Energy Phys. **10**, 014 (2003). https://doi.org/10.1088/1126-6708/2003/10/014

20. J.H. MacGibbon, B.R. Webber, Phys. Rev. D **41**, 3052 (1990). https://doi.org/10.1103/PhysRevD.41.3052

21. E. Cotner, A. Kusenko, Phys. Rev. D **96**(10), 103002 (2017). https://doi.org/10.1103/PhysRevD.96.103002

22. K.S. Thorne, Astrophys. J. **191**, 507 (1974). https://doi.org/10.1086/152991

23. M. Kesden, G. Lockhart, E.S. Phinney, Phys. Rev. D **82**(12), 124045 (2010). https://doi.org/10.1103/PhysRevD.82.124045

24. R. Dong, W.H. Kinney, D. Stojkovic, J. Cosmol. Astropart. Phys. **10**, 034 (2016). https://doi.org/10.1088/1475-7516/2016/10/034

25. D.N. Page, Phys. Rev. D **16**, 2402 (1977). https://doi.org/10.1103/PhysRevD.16.2402

26. J. Ahmed, K. Saifullah, Eur. Phys. J. C **78**(4), 316 (2018). https://doi.org/10.1140/epjc/s10052-018-5800-6

27. J.V. Rocha, J. High Energy Phys. **2009**(8), 027 (2009). https://doi.org/10.1088/1126-6708/2009/08/027

28. P. Kanti, T. Pappas, Phys. Rev. D **96**(2), 024038 (2017). https://doi.org/10.1103/PhysRevD.96.024038

29. C.Y. Zhang, P.C. Li, B. Chen, Phys. Rev. D **97**(4), 044013 (2018). https://doi.org/10.1103/PhysRevD.97.044013

30. P.C. Li, C.Y. Zhang, Phys. Rev. D **99**(2), 024030 (2019). https://doi.org/10.1103/PhysRevD.99.024030

31. L.C.B. Crispino, A. Higuchi, E.S. Oliveira, J.V. Rocha, Phys. Rev. D **87**(10), 104034 (2013). https://doi.org/10.1103/PhysRevD.87.104034

32. R. Dong, J. Sakstein, D. Stojkovic, Phys. Rev. D **96**(6), 064048 (2017). https://doi.org/10.1103/PhysRevD.96.064048

33. C. Charmousis, M. Crisostomi, D. Langlois, K. Noui, arXiv e-prints. arXiv:1907.02924 (2019)

34. C.M. Harris, P. Kanti, J. High Energy Phys. **2003**(10), 014 (2003). https://doi.org/10.1088/1126-6708/2003/10/014

35. G. Duffy, C.M. Harris, P. Kanti, E. Winstanley, J. High Energy Phys. **2005**(9), 049 (2005). https://doi.org/10.1088/1126-6708/2005/09/049

36. M. Casals, P. Kanti, E. Winstanley, J. High Energy Phys. **2006**(2), 051 (2006). https://doi.org/10.1088/1126-6708/2006/02/051

37. M. Casals, S. Dolan, P. Kanti, E. Winstanley, J. High Energy Phys. **2007**(3), 019 (2007). https://doi.org/10.1088/1126-6708/2007/03/019

38. Y.H. Hyun, Y. Kim, S.C. Park, J. High Energy Phys. **2019**(6), 41 (2019). https://doi.org/10.1007/JHEP06(2019)041

39. R. Dong, D. Stojkovic, Phys. Rev. D **92**(8), 084045 (2015). https://doi.org/10.1103/PhysRevD.92.084045
40. P. Boonserm, T. Ngampitipan, P. Wongjun, Eur. Phys. J. C **78**(6), 492 (2018). https://doi.org/10.1140/epjc/s10052-018-5975-x

41. J.H. MacGibbon, Phys. Rev. D **44**, 376 (1991). https://doi.org/10.1103/PhysRevD.44.376

42. D.C. Dai, D. Stojkovic, J. High Energy Phys. **2010**, 16 (2010). https://doi.org/10.1007/JHEP08(2010)016

43. D.N. Page, Phys. Rev. D **13**, 198 (1976). https://doi.org/10.1103/PhysRevD.13.198

44. D.N. Page, Phys. Rev. D **14**, 3260 (1976). https://doi.org/10.1103/PhysRevD.14.3260

45. J. Bellm et al., Eur. Phys. J. C **76**(4), 196 (2016). https://doi.org/10.1140/epjc/s10052-016-4018-8

46. T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C.O. Rasmussen, P.Z. Skands, Comput. Phys. Commun. **191**, 159 (2015). https://doi.org/10.1016/j.cpc.2015.01.024

47. K. Kohri, J. Yokoyama, Phys. Rev. D **61**(2), 023501 (2000). https://doi.org/10.1103/PhysRevD.61.023501

48. A. Arbey, F. Mahmoudi, Comput. Phys. Commun. **181**, 1277 (2010). https://doi.org/10.1016/j.cpc.2010.03.010

49. A. Arbey, F. Mahmoudi, Comput. Phys. Commun. **182**, 1582 (2011). https://doi.org/10.1016/j.cpc.2011.03.019

50. A. Arbey, F. Mahmoudi, G. Robbins, Comput. Phys. Commun. **239**, 238 (2019). https://doi.org/10.1016/j.cpc.2019.01.014

51. S. Chandrasekhar, S. Detweiler, Proc. R. Soc. Lond. Ser. A **345**, 145 (1975). https://doi.org/10.1098/rspa.1975.0130

52. S. Chandrasekhar, S. Detweiler, Proc. R. Soc. Lond. Ser. A **350**, 165 (1976). https://doi.org/10.1098/rspa.1976.0101

53. S. Chandrasekhar, Proc. R. Soc. Lond. Ser. A **348**, 39 (1976). https://doi.org/10.1098/rspa.1976.0022

54. S. Chandrasekhar, S. Detweiler, Proc. R. Soc. Lond. Ser. A **352**, 325 (1977). https://doi.org/10.1098/rspa.1977.0002

55. B.E. Taylor, C.M. Chambers, W.A. Hiscock, Phys. Rev. D **58**(4), 044012 (1998). https://doi.org/10.1103/PhysRevD.58.044012