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

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.


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][2][3], direct observation of the Milky Way supermassive central BH [4], and the cosmological and gravitational questions related to primordial BHs (see for example [5][6][7][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 emission of BHs, which however focus on higherdimensional 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 Hawking evaporation radiation into stable or long-lived particles of 4−dimensional BHs mass distributions and its evolution in time.
This document constitutes the manual of BlackHawk v1.0 and is organized as follows: Section 2 is a brief overview of BHs and Hawking evaporation physics, Section 3 presents the structure and file content of the code, and the compilation and run instructions. Section 4 describes the input parameters needed to run BlackHawk, Section 5 gives a detailed description of all the routines written in the code. Section 6 follows the normal execution of BlackHawk programs and gives examples of screen output. Section 7 presents the format of the data files generated by a run along with examples, Section 8 gives an estimation of the memory usage and Section 9 provides instructions for the users on how to modify the code.

Physics of Hawking evaporation
In this section we give a short overview of the main physical aspects of Hawking evaporation. This concerns BHs of primordial origin (Primordial Black Holes, denoted as PBHs), as well as any other BHs.
In the following, all formulas are in natural units where = c = k B = G = 1, unless stated otherwise.

Testing BH 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 Hawking emission 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]. 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.

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 for example [12]): where n ≈ 1.3 and R c is measured using the Cosmic Microwave Background (CMB) to be R c = (24.0 ± 1.2) × 10 −10 at the scale k 0 = 0.002 Mpc −1 . The comoving number density of BHs resulting from this power spectrum is obtained in [12] through peaktheory: where: and: (2.4) in which • H 0 = 67.8 km·s −1 ·Mpc −1 is the current Hubble parameter [13].
• Ω m = 0.308 is the matter mass fraction in the Universe [13].
• g * eq = 3.36 is the number of relativistic energy degrees of freedom (dof) at radiation-matter equality [12].
• g * = 106.75 is the number of relativistic energy dof at the time of BH formation (here the end of the inflation) [13]; • ζ th = 0.7 parametrizes the direct collapse of a density fluctuation into a BH [12].

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: where A is the amplitude, M c is the position of the peak and σ is its width.

Power-law distribution
The power-law distribution [11] is a less refined version of Eq. (2.2). It also derives from scale-invariant primordial density fluctuations and is given by: where γ = −2w/(1 + w) and w is defined through the equation of state of the dominating energy in the Universe at the epoch of BH formation such as P = wρ.

Critical collapse distribution
The critical collapse distribution [11] derives from a Dirac power spectrum for primordial density fluctuations. It is defined as: where A is an amplitude factor and M f an upper cut-off.

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 cm 3 .

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 [14] that BH horizons emit elementary particles as blackbodies with a temperature linked to their mass M through 4 : The number of particles emitted per units of time and energy is: where the sum is over the number of quantum dof (see Table 2 in Appendix C) and the ± are for fermions and bosons, respectively. The factor Γ 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: To obtain instantaneous quantities for a single BH of mass M 0 , one just needs to take: The greybody factors describe the probability that a spherical wave representing 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 µ: 12) in the Schwarzschild metric: where ∆(r) ≡ r 2 h(r), K(r) ≡ r 2 E 2 and E is the particle frequency (or equivalently its energy). In this equation, the separation constant λ ls ≡ 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 σ(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 [18]: The method used in BlackHawk to compute those greybody factors is described in Appendix B.1.

Kerr Black Holes
Kerr Black Holes are an extension of the Schwarzschild ones with an additional parameter: their spin a ≡ J/M ∈ [0, M] (in the following we will denote the reduced spin parameter by a * ≡ a/M ∈ [0, 1]) where J is the BH angular momentum. These rotating BHs could gain their spin through their formation mechanism [19], accretion [20] or merging process [21]. They are axially symmetric and require a specific treatment.
The temperature of a rotating BH is given by: where r + ≡ M + √ M 2 − a 2 is the Kerr external radius. The Teukolsky equation (2.14) has to be modified with ∆(r) ≡ r 2 − 2Mr + a 2 and K(r) ≡ (r 2 + a 2 )E 2 + am, where m is the projection of the angular momentum l. The separation constant λ 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 γ = a * ME, as given in [22] 5 .
The number of particles emitted per units of time and energy is now: where E ′ ≡ E − mΩ and Ω ≡ a * /(2r + ) is the angular velocity at the horizon [22]. The method used to compute these greybody factors in BlackHawk is also described in Appendix B.1.

Schwarzschild Black Hole
Once the greybody factors are known, it is possible to integrate Eq. (2.9) to obtain a differential equation for the mass loss of a BH through Hawking evaporation [23]: The factor f (M) accounts for the number of quantum dof that a BH of mass M can emit. It is obtained through [22]: The computation of the f (M) factor in BlackHawk is described in Appendix B.2.

Kerr Black Hole
For Kerr BH, 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. The equation for the factor f (M, a * ) becomes: 20) and the differential equation describing the angular momentum J is [22] Once the f (M, a * ) and g(M, a * ) factors are obtained, the evolution of a * is straightforwardly obtained through: The computation of the f (M, a * ) and g(M, a * ) factors in BlackHawk are described in Appendix B.2.

Hadronization
The elementary particles emitted by BHs are not the final products of the Hawking emission. 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 [24] and PYTHIA [25] 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 −8 s. These particles are listed in the Table 2 of Appendix C.
The time-dependent comoving density of Hawking secondary particle i emitted by a distribution of BHs per units of time and energy is computed with the integral: where the sum is taken over Hawking primary particles, and Section B.3 describes how hadronization tables dN i j (E ′ , E) have been computed to transform the primary spectra into secondary spectra in BlackHawk.

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 and Windows (using Cygwin64).

Main directory
The main directory contains: • the source codes BlackHawk_inst.c and BlackHawk_tot.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.

src/ subfolder
This folder contains: • a header file include.h containing the declaration of all routines along with the parameter structure struct param (see Section 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, hadro_pythia.c, hadro_pythianew.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.

results/ subfolder
This folder is designed to receive subfolders of data generated by running the BlackHawk code (see Section 7).

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

scripts/ subfolder
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 [26][27][28]. These scripts can be used to generate the needed tables. They are accompanied by README.txt files explaining how to use them.

Compilation
The compilation of BlackHawk has been tested on Linux and Windows (using Cygwin64) distributions. The code is written in C99 standard. To compile the code, simply cd into the main directory and type 7 : make BlackHawk_* where * denotes tot or inst. This will create a library file libblackhawk.a and an executable file 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 8 : where parameter_file is the name of a parameter file. To compile only the library, just cd into the main directory and type: make 7 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. 8 In case of memory problem at execution, increasing the stack size with the command ulimit -s unlimited can help solving the problem.

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

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: Most routines described in Section 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.

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 (full-_output = 1) or not (full_output = 0). It can be useful for debugging the code or seeing 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).

BH spectrum parameters
The second set of parameters defines the quantities used to compute the BH density distribution (see Section 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 Section 9). 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 ≈ 2×10 −5 g < Mmin, Mmax, where M p is the Planck mass. For a mass distribution, one must have Mmin < Mmax. If they are not compatible with boundaries of the mass distribution, the computation will stop (see below).
• amplitude is the amplitude A present in Eqs. (2.5), (2.6) and (2.7). It is the normalization of the corresponding BH distribution and thus strictly positive.
• variance is the variance σ in the log-normal distribution of Eq. (2.5). It has to be strictly positive.
• crit_mass is the characteristic mass M c in Eq. (2.5) and M f in Eq. (2.7). It has to be strictly positive.
• eq_state defines the equation of state w (see Section 2.1.3).
• table is the name of a user-defined BH distribution table. It has to be a string with any file extension.

BH evolution parameters
The next set of parameters defines the quantities used to compute the BH evolution (see Section 5.3): • tmin is the initial integration time of the evolution of BH, in seconds. It can have any positive value.
• 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 Section 5.3). It is fixed to 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 Section 9).
• 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 Section 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 Section 9).

Primary spectrum parameters
This set of parameters defines the quantities related to the primary Hawking spectra (see Section 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 ± 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 Section 9).
• grav determines whether the emission of gravitons by BH will be taken into account (grav = 1) or not (grav = 0).
• nb_gamma_a and nb_gamma_x are respectively the number of spins and x ≡ 2 × E × M tabulated in the greybody factor tables. They should not be modified unless the corresponding tables are recomputed (see Section 9).

Hadronization parameters
This last set of parameters defines the quantities used during the hadronization (see Section 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.
• hadronization_choice determines which hadronization tables will be used to compute the secondary spectra (see Section 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 Section 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 Section 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 Section 9).

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).

General routines
There are 4 general routines in the BlackHawk code. The principal ones are the main routines, described in Section 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 Section 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 Section 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 be performed.
• int memory_estimation(struct param *parameters, int session): this routine gives a rough estimate of the usage of both RAM and disk space (see Section 8). If the user decides to cancel the run the value 0 is returned, otherwise it is 1. The output is given in MB.

BH spectrum routines
There are 4 routines contributing to the BH initial spectrum computation (see Section 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 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.

BH evolution routines
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 −3 , 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: • void write_life_evolutions(double **life_masses, double **life_--spins, double **life_times, int *evolution_length, struct param *parameters): this routine writes the BH time-dependent masses until full evaporation in the file life_evolutions.txt, saved in destination_folder/ (see Section 7.1). The results are converted from GeV to CGS units.

Primary spectra routines
There are 5 routines contributing to the computation of the primary Hawking spectra (see Section 2.2): • 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 Section 7.2). The results are converted from GeV to CGS units.

Secondary spectrum routines
There are 9 routines contributing to the computation of the secondary Hawking spectra (see Section 2.4): • • void total_spectra(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 steptime 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 time step.
• 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 *_primary_spectrum.txt and *_secondary_spectrum.txt files. • 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[][].

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 use in the same destination_folder/ because the output files will not enter in conflict (see Section 7). We will now describe the structure of the main routines together with screen output examples.

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 −10 defines the precision up to which two double numbers are considered as equal.
• Mp ≡ 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 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 Section 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 if of the form 9  The subsequent execution steps depend on the program. Output examples are given in the mode full_output = 0.

BlackHawk_tot: Full 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

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 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 Section 5.5).
This is the end of the execution of BlackHawk_inst. The expected output is of the form:

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.

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 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 ).
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.

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 the routine 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 the routine write_instantaneous_hadronized_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.

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.

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 Section 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): Using the parameters of Appendix D.1, the arrays occupy at most ∼ 10 MB.

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 × 3 × BHnumber.
Using the parameters of Appendix D.1, the total written disk space is ∼ 35 kB.

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

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 src/tables/ subfolder; • 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 subfolder scripts/ together with README files.

Using another BHs initial distribution
The user may be interested in testing its own BHs 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 subfolder src/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 spectrum_choice to 5, putting the table file in the subfolder 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.

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 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 (for example the file names or columns in the writing routines), • eventually add its(their) contribution(s) to the secondary spectra.

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 (for example the file names or columns in the writing routines), • add the corresponding contribution(s) to the routine contribution_instantaneous.

Conclusion
BlackHawk is the first public code generating both primary and secondary Hawking evaporation spectra for any 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 evolution is obtained by considering the energy loss via Hawking radiation and the 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 evaporation on observable quantities and thus to disqualify or set constraints on cosmological models implying the formation of black holes, as well as to test the Hawking radiation assumptions and study black hole general properties.

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 = c = k B = G = 1 to CGS or SI.

A.1 Energy
The energy conversion from GeV to Joule is:

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

A.3 Time
The dimensional link between energy and time is [t] = [ /E], and the conversion from GeV to seconds is:

A.4 Distance
The dimensional link between energy and distance is [l] = [ c/E], and the conversion from GeV to meters is:

A.5 Temperature
The dimensional link between energy and temperature is [T ] = [E/k B ], and the conversion from GeV to Kelvins is:

B Computation of the tables B.1 Greybody factors
Chandrasekhar and Detweiler have shown that the Teukolsky equation can be reduced to a wave equation for Kerr Black Holes [29][30][31][32]. It is indeed difficult to find short-range potentials allowing for precise numerical computation. They give the form of such potentials in [29,30] for spin 2, [31] for spins 0 and 1 and [32] for spin 1/2, and found necessary to define a modified Eddington-Finkelstein radial coordinate r * by the equation: where ρ(r) 2 ≡ r 2 +α 2 and α 2 ≡ 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: Unfortunately, the inverse of this equation has to be found numerically and is generally difficult to determine with accurate precision. As boundary conditions, we use a purely outgoing wave. The solution at the horizon has the form: At infinity, the solution has the form: The Schrödinger-like wave equation is for all spins: The method to transform Eq. (2.14) into this simple wave equation was proposed in the Chandrasekhar-Detweiler papers [29][30][31][32]. The potentials are 10 : 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) = (4ν + 6)ρ 2 + 8νr 2 − 6r 2 + 36Mr − 12a 2 , (B.12) q ′ ∆−q∆ ′ = −2(r−M)νρ 4 +2ρ 2 (2νr∆−3M(r 2 +a 2 )+6ra 2 )+12r∆(Mr−a 2 ) , (B.13) β ± = ±3α 2 , (B.14) where ν ≡ λ 2 lm + 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 2 div = −α 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.5) up to slightly before the divergence (e.g. r div − ǫ). At this point, the behaviour of the potential V s is known, and Eq. (B.5) is simplified. Since the form of the function ψ s can be obtained, by continuity of the function R s of Eq. (2.14) one can extrapolate this form up to slightly after the divergence (e.g. r div + ǫ) 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 − β ± ∆ term. For this extra divergence, we try to integrate for a given spin with one of the potentials (e.g. κ + , β + ), and in case of problem we try with the other potentials (e.g. κ + , β − ), as it seems that at least one of the four potential does not generate any divergence.
The greybody factor is given by the transmission coefficient of the wave from the horizon to the infinity: Practically, we compute the value of the single dof emissivities: for some values of a * and for a range of 0.01 < x ≡ 2Er 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: 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 [18] in the Schwarzschild case and of [33] in the Kerr case. The Mathematica scripts spin_*.m, the fitting script exploitation.m as well as a C formating script formating.c and a README are provided in the subfolder: scripts/greybody_scripts/greybody_factors/

B.2 Evolution tables
To compute the integrals of Eqs. (2.20) and (2.21), 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 [18] for example), thus the integral does not need to be computed over all energies, but a restrained set 10 −5 × T < E < 10 5 × 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). Masses are given in GeV (corresponding to grams) and f (M, a * ) and g(M, a * ) are in GeV 4 (corresponding to g 3 ·s −1 and g 2 ·GeV·s −1 respectively). We have checked that the value of f (M, 0) is consistent with that of [23] in the Schwarzschild case and of [34] in the Kerr case.
The C script fM.c used to compute the tables and a README are provided in the subfolder: scripts/greybody_scripts/fM/

B.3 Hadronization
Two particle physics codes have been used to compute hadronization tables: HERWIG [24] and PYTHIA [25]. In both cases, the strategy is to generate the output of a collision (for example e + + e − → u + u → ...), 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). It 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 particle PYTHIA (new) 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 highenergy 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 formating.c used to format the hadronization tables and README files are provided in the subfolders: scripts/pythia_scripts/ scripts/herwig_scripts/ scripts/pythia_scripts_new/ Please contact one of the authors if you have issues using these scripts.

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

D.1 Parameters
An example of parameters.txt file is given here: