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

TL;DR

BlackHawk v2.0 is a publicly available C program that enables detailed calculations of Hawking radiation spectra from any black hole distribution, facilitating studies of black hole evaporation and particle emission.

Contribution

The paper introduces BlackHawk v2.0, a comprehensive and user-friendly tool for computing Hawking radiation spectra for diverse black hole populations, with new features and improved capabilities.

Findings

Enables calculation of primary and secondary particle spectra from black holes.

Supports evolution of black hole distributions over time.

Provides a manual and examples for user guidance.

Abstract

We describe BlackHawk, a public C program for calculating the Hawking evaporation spectra of any black hole distribution. This program allows 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. This is the BlackHawk v2.0 manual, which is available on the BlackHawk webpage http://blackhawk.hepforge.org/. A brief release note summarizing the new aspects of BlackHawk v2.0 as well as illustrating examples can be found in https://arxiv.org/abs/2108.02737.

Tables

Table 1. Table 1: Equation of state for different cosmological fluids.

Fluid	Equation of state
Matter	$w=0$
Radiation	$w=1/3$
Cosmological constant	$w=-1$
Curvature	$w=-1/3$

Table 2. Table 2: Horizon temperatures for different BH metrics.

Metric	Temperature
Schwarzschild	$T_S=\frac{1}{4\pi r_S}$
Reissner–Nordström	$T_Q=\frac{r_{+}-r_{-}}{4\pi r_{+}^2}$
Higher-dimensional	$T_n=\frac{n+1}{4\pi r_H}$
Polymerized	$T_{\mathrm{LQG}}=\frac{r_{+}^2\left(r_{+}-r_{-}\right)}{4\pi \left(r_{+}^4+a_0^2\right)}$
Kerr	$T_K=\frac{1}{2\pi }\left(\frac{r_{+}-M}{r_{+}^2+a^2}\right)$

Table 3. Table 3: Properties of the elementary particles of the Standard Model, in addition to the graviton  [ 35 ] . 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 Λ ≈ 200 Λ 200 \Lambda\approx 200\, MeV.

particle	symbol	mass (GeV/c2)	spin	quantum dof $g_i$
Higgs boson	$h^0$	$1.2503\times {10}^2$	$0$	$1$
photon	$\gamma$	$0$	$1$	$2$
gluons	$g$	$0$	$1$	$16$
W bosons	$W^{\pm }$	$8.0403\times {10}^1$	$1$	$6$
Z boson	$Z^0$	$9.11876\times {10}^1$	$1$	$3$
neutrinos	${\nu }_{e,\mu ,\tau },{\overline{\nu }}_{e,\mu ,\tau }$	$0$	$1/2$	$6$
electron	$e^{\pm }$	$5.109989461\times {10}^{-4}$	$1/2$	$4$
muon	${\mu }^{\pm }$	$1.056583745\times {10}^{-1}$	$1/2$	$4$
tau	${\tau }^{\pm }$	$1.77686$	$1/2$	$4$
up quark	$u,\overline{u}$	$2.2\times {10}^{-3}$	$1/2$	$12$
down quark	$d,\overline{d}$	$4.7\times {10}^{-3}$	$1/2$	$12$
charm quark	$c,\overline{c}$	$1.27$	$1/2$	$12$
strange quark	$s,\overline{s}$	$9.6\times {10}^{-2}$	$1/2$	$12$
top quark	$t,\overline{t}$	$1.7321\times {10}^2$	$1/2$	$12$
bottom quark	$b,\overline{b}$	$4.18$	$1/2$	$12$
graviton	$G$	$0$	$2$	$2$

Table 4. Table 4: Particles with a lifetime longer than 10 − 8 superscript 10 8 10^{-8}\, s  [ 35 ] , relevant for BBN studies and used to compute the hadronization tables in pythia_tables/ and herwig_tables/ .

particle	symbol	mass (GeV/c2)	lifetime (s)
photon	$\gamma$	$0$	$\infty$
electron	$e^{\pm }$	$5.109989461\times {10}^{-4}$	$\infty$
muon	${\mu }^{\pm }$	$1.056583745\times {10}^{-1}$	$\left(2.1969811\pm 0.0000022\right)\times {10}^{-6}$
neutrinos	${\nu }_{e,\mu ,\tau },{\overline{\nu }}_{e,\mu ,\tau }$	0	$\infty$
charged pions	${\pi }^{\pm }$	$1.3957018\times {10}^{-1}$	$\left(2.6033\pm 0.0005\right)\times {10}^{-8}$
neutral “long” kaon	$K_L^0$	$4.977\times {10}^{-1}$	$\left(5.099\pm 0.021\right)\times {10}^{-8}$
charged kaons	$K^{\pm }$	$4.937\times {10}^{-1}$	$\left(1.2379\pm 0.0021\right)\times {10}^{-8}$
proton	p,$\overline{p}$	$9.38272\times {10}^{-1}$	$\infty$
neutron	n,$\overline{n}$	$9.395654\times {10}^{-1}$	$880.2\pm 1$

Table 5. Table 5: Stable particles  [ 35 ] , relevant for evaporating BHs in the low redshift universe and used to compute the hadronization table in pythia_tables_new/ and hazma_tables/ .

particle	symbol	mass (GeV/c2)	lifetime (s)
photon	$\gamma$	0	$\infty$
electron	$e^{\pm }$	$5.109989461\times {10}^{-4}$	$\infty$
neutrinos	${\nu }_{e,\mu ,\tau },{\overline{\nu }}_{e,\mu ,\tau }$	0	$\infty$
proton	p,$\overline{p}$	$9.38272\times {10}^{-1}$	$\infty$

Table 6. Table 6: List of the channels used to compute the hadronization tables for the Monte-Carlo codes.

particle	PYTHIA (new)	HERWIG
gluons	$e^{+}{e}^{-}\to h^0\to g\overline{g}$	$e^{+}{e}^{-}\to h^0\to g\overline{g}$
Higgs boson	$e^{+}{e}^{-}\to h^0$	$e^{+}{e}^{-}\to h^0$
W bosons	$e^{+}{e}^{-}\to Z^0/{\gamma }^{\ast }\to W^{+}{W}^{-}$	$e^{+}{e}^{-}\to Z^0/{\gamma }^{\ast }\to W^{+}{W}^{-}$
Z boson	$e^{+}{e}^{-}\to h^0\to Z^0{Z}^0$	$e^{+}{e}^{-}\to Z^0/{\gamma }^{\ast }\to Z^0{Z}^0$
leptons	$e^{+}{e}^{-}\to h^0\to l^{+}{l}^{-}$	$e^{+}{e}^{-}\to Z^0/{\gamma }^{\ast }\to l^{+}{l}^{-}$
quarks	$e^{+}{e}^{-}\to Z^0/{\gamma }^{\ast }\to q\overline{q}$	$e^{+}{e}^{-}\to Z^0/{\gamma }^{\ast }\to q\overline{q}$

Equations

$$t_{\mathrm{f}}\sim \frac{M}{\gamma },$$

$$n_{\mathrm{BH}}={\int }_{M_{\min }}^{M_{\max }}\frac{\mathrm{d}n}{\mathrm{d}M}\mathrm{d}M,$$

$${\rho }_{\mathrm{BH}}={\int }_{M_{\min }}^{M_{\max }}M\frac{\mathrm{d}n}{\mathrm{d}M}\mathrm{d}M.$$

$$P(k)=R_{\mathrm{c}}{\left(\frac{k}{k_0}\right)}^{n-1},$$

$$\mathrm{d}n\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}\mathrm{d}M,$$

$$\nu (M)\equiv {\left(\frac{2(k_0^{2}M/X{)}^{(n-1)/2}}{R_{\mathrm{c}}\Gamma ((n-1)/2)}\right)}^{1/2}{\zeta }_{\mathrm{th}},X\equiv \frac{4\pi }{3}{\left(\frac{8\pi G}{3}\right)}^{-1}{\left\{\frac{H_0^2{\Omega }_{\mathrm{m}}}{1+z_{\mathrm{eq}}}{\left(\frac{g_{\ast ,\mathrm{eq}}}{g_{\ast }}\right)}^{1/3}\right\}}^{1/2},$$

$$\mathrm{d}n=\frac{A}{\sqrt{2\pi }\sigma M^2}\exp \left(-\frac{\ln (M/M_{\mathrm{c}}{)}^2}{2{\sigma }^2}\right)\mathrm{d}M,$$

$$\mathrm{d}n=\frac{A}{\sqrt{2\pi }\sigma M}\exp \left(-\frac{\ln (M/M_{\mathrm{c}}{)}^2}{2{\sigma }^2}\right)\mathrm{d}M.$$

$$\mathrm{d}n=A{M}^{\gamma -2}\mathrm{d}M,$$

$$\mathrm{d}n=A{M}^{1.85}\exp \left(-{\left(\frac{M}{M_{\mathrm{f}}}\right)}^{2.85}\right)\mathrm{d}M,$$

$$\mathrm{d}n=\frac{A}{M_{\max }-M_{\min }}\mathrm{d}M,$$

$$\mathrm{d}n=A\delta (M-M_{\mathrm{c}}).$$

$$a^{\ast 2}+Q^{\ast 2}<1,$$

$${\tilde{n}}_{\mathrm{BH}}={\int }_0^1\frac{\mathrm{d}\tilde{n}}{\mathrm{d}x}\mathrm{d}x=1.$$

$$\mathrm{d}\tilde{n}=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(x-x_{\mathrm{c}}{)}^2}{2{\sigma }^2}\right)\mathrm{d}x,$$

$$\mathrm{d}\tilde{n}=\frac{\mathrm{d}x}{x_{\max }-x_{\min }},$$

$$\mathrm{d}\tilde{n}=\delta (x-x_{\mathrm{c}}).$$

$$\mathrm{d}{s}^2=-G(r)\mathrm{d}{t}^2+\frac{1}{F(r)}\mathrm{d}{r}^2+H(r)\mathrm{d}{\Omega }^2,$$

$$F(r),G(r)\underset{r\to +\infty }{⟶}1,H(r)\underset{r\to +\infty }{\sim }{r}^2,$$

$$F(r)=G(r)=1-\frac{r_{\mathrm{H}}}{r},H(r)=r^2.$$

$$F(r)=G(r)=1-\frac{r_{\mathrm{S}}}{r}+\frac{r_{\mathrm{Q}}^2}{r^2},H(r)=r^2,$$

$$r_{\mathrm{H},\mathrm{C}}\equiv r_{\pm }\equiv r_{\mathrm{S}}\frac{1\pm \sqrt{1-4r_{\mathrm{Q}}^2/r_{\mathrm{S}}^2}}{2}.$$

$$F(r)=G(r)=1-{\left(\frac{r_{\mathrm{H}}}{r}\right)}^{n+1},H(r)=r^2,$$

$$r_{\rm H}=\dfrac{1}{\sqrt{\pi}M_{*}}\left(\dfrac{M}{M_{*}}\right)^{1/(n+1)}\left(\dfrac{8\Gamma\big{(}(n+3)/2\big{)}}{n+2}\right)^{1/(n+1)},$$

$$F(r)=\frac{(r-r_{+})(r-r_{-})r^4}{(r+r_{\ast }{)}^2(r^4+a_0^2)},G(r)=\frac{(r-r_{+})(r-r_{-})(r+r_{\ast }{)}^2}{r^4+a_0^2},H(r)=r^2+\frac{a_0^2}{r^2},$$

$$r_{\mathrm{H}}\equiv r_{+}\equiv 2m(\epsilon ),r_{\mathrm{C}}\equiv r_{-}\equiv 2m(\epsilon )P(\epsilon {)}^2,$$

$$P(\epsilon )=\frac{\sqrt{1+{\epsilon }^2}-1}{\sqrt{1+{\epsilon }^2}+1},$$

$$M=m(\epsilon )(1+P(\epsilon ){)}^2.$$

$${\rm d}s^{2}=\big{(}{\rm d}t-a\sin^{2}(\theta){\rm d}\phi\big{)}^{2}\,\frac{\Delta}{\Sigma}-\left(\frac{{\rm d}r^{2}}{\Delta}+{\rm d}\theta^{2}\right)\Sigma-\big{(}(r^{2}+a^{2}){\rm d}\phi-a{\rm d}t\big{)}^{2}\,\frac{\sin^{2}(\theta)}{\Sigma}\,,$$

$$r_{\mathrm{H},\mathrm{C}}\equiv r_{\pm }\equiv r_{\mathrm{S}}\frac{1\pm \sqrt{1-a^{\ast 2}}}{2}.$$

Full text

1. I Introduction

2. II Hawking radiation

3. II.1 Black hole distributions

4. II.1.1 Mass distributions

5. II.1.2 Secondary parameter distribution

6. II.2 Hawking evaporation physics

7. II.2.1 Spherically symmetric and static black holes

8. II.2.2 Rotating black holes

9. II.2.3 Particle creation by black holes

10. II.2.4 Black hole evolution

11. II.2.5 Hadronization and decay

12. III Content and compilation

13. III.1 Content

14. III.2 Compilation

15. IV Input parameters

16. IV.1 Parameter structure

17. IV.2 General parameters

18. IV.3 BH spectrum parameters

19. IV.4 BH evolution parameters

20. IV.5 Primary spectrum parameters

21. IV.6 Hadronization parameters

22. IV.7 Tables parameters

23. V Routines

24. V.1 General routines

25. V.2 BH initial distribution routines

26. V.3 BH evolution routines

27. V.4 Primary spectra routines

28. V.5 Secondary spectrum routines

29. V.6 Auxiliary routines

30. V.7 General feature

31. VI BlackHawk programs

32. VI.1 Common features

33. VI.2 BlackHawk_tot: Full time-dependent Hawking spectra

34. VI.3 BlackHawk_inst: Instantaneous Hawking spectra

35. VII Output files

36. VII.1 BlackHawk_tot program

37. VII.2 BlackHawk_inst program

38. VIII Memory usage

39. VIII.1 RAM used

40. VIII.2 Static disk memory used

41. IX Other applications

42. IX.1 Computing new numerical tables

43. IX.2 Using another BH mass and spin distribution

44. IX.3 Adding primary particles

45. IX.4 Adding secondary particles

46. IX.5 Other types of black holes

47. X Conclusion

48. A Units

49. A.1 Energy

50. A.2 Mass

51. A.3 Time

52. A.4 Distance

53. A.5 Temperature

54. A.6 Electric charge

55. B Particle information

56. C Numerical tables

57. C.1 Greybody factors

58. C.1.1 General framework

59. C.1.2 Kerr black holes

60. C.1.3 Spherically symmetric and static black holes

61. C.1.4 Boundary conditions

62. C.2 Page factors

63. C.3 Hadronization tables

64. D Example runs

65. D.1 Built-in parameters.txt file

66. D.2 BlackHawk_tot output files

67. D.3 BlackHawk_inst output files

I Introduction

Black holes (BHs) may be the most promising astrophysical objects in the study of beyond general relativity (GR) and beyond Standard Model (SM) physics. BHs of stellar mass $M\sim 5-100\,M_{\odot}$ have been observed indirectly from X-ray emission of binaries for quite some time (see e.g. [2] and references therein). More recently, gravitational wave (GW) instruments such as LIGO/VIRGO have detected the perturbation of spacetime generated by the merging of stellar mass BHs, confirming GR predictions with great accuracy [3, 4, 5, 6]. Future instruments may improve over these catalogs, giving access to relevant BH population statistics and to an extended mass and binary separation ranges with unprecedented precision (these include KAGRA [7], the Einstein Telescope [8] and LISA [9]). Supermassive BHs with $M\gg 10^{6}\,M_{\odot}$ are believed to lie at the centre of every big galaxy and have been detected by radio interferometers; such as the central BH of M87 (see [10], and subsequent papers for details). There are no observation in the intermediate mass range yet, but GW interferometers could detect $M\sim 10^{2}-10^{3}\,M_{\odot}$ in the near future [11]. Interestingly, the first BH population statistics extracted from these catalogs are not easily explained by stellar evolution. Thus, some claim that part of the observed stellar mass BHs are in fact of primordial origin, that is to say BHs formed in the early universe through the collapse of overdensities (for recent reviews on primordial BHs — PBHs — see e.g. [12] and [13]). This channel of formation seems to be consistent with the LIGO/VIRGO observations if dark matter (DM) is composed at most of $\sim 1\%$ of stellar mass PBHs [12, 14]. Future observations such as the precise measure of the spin of the binary BHs components [15], the stochastic GW background [16] or detection of BHs with masses outside the stellar evolution predictions [17, 18] could point toward a primordial origin for some of those; while it is still not clear how supermassive BHs can be present at such high redshifts: they may have been formed from the hierarchical mergers of massive PBH seeds [19]. On the other extreme of the mass range, PBHs could be as small as the Planck mass $M\sim 10^{-5}\,$g, even if some inflation models place lower bounds of $M\gtrsim 10^{-1}\,$g [20, 12]. Their detection represents a challenge of modern astrophysics (see the review [12]).

Hawking proved that BHs emit radiation as blackbodies, therefore slowly losing mass during a process called Hawking evaporation [21]. One major feature of that process is that Hawking radiation (HR) gets more and more energetic as the BH gets lighter, reaching the Planck energy at the end of the BH evaporation. Of course, BHs with initial mass $M\lesssim M_{\odot}$ (below the Oppenheimer–Volkoff limit) are presumably of primordial origin, but even stellar mass BHs would evaporate away within an incredibly long time. Tiny PBHs of $\sim 10^{15}-10^{18}\,$g emit gamma rays that could be detected in Earth of space based observatories. PBHs even smaller, with initial mass $M\lesssim 10^{15}\,$g, would have already disappeared and can not represent a fraction of DM today. Their evaporation in the early universe could have left imprints in cosmological observables [12]. Larger PBHs with $M\gtrsim 10^{18}\,$g still emit HR but with a power that is so small that detection is considered out of reach for now (see however the proposition [22]). Historically, Hawking radiation has been considered as a mean of detection of tiny BHs evaporating instantaneously in high energy colliders such as the LHC, within models where extra spatial dimensions reduce the Planck mass and thus the energy density required to form a BH. The public codes BlackMax [23] and Charybdis [24] were published in that context: they compute the spectrum of particles generated by the evaporation of higher-dimensional (spinning) BHs in a LHC-like detector. As HR is a purely gravitational phenomenon, all particles of the (beyond the) SM spectrum are emitted. These particles are subsequently subject to (beyond the) SM interactions such as hadronization, inter-conversion and decay, before reaching some detector. Hence, particle physics codes such as PYTHIA [25], HERWIG [26] or Hazma [27] have been used to compute the secondary spectra of particles resulting from this primary emission. There is therefore a strong link between HR and particle physics, particularly in the QCD domain [28, 29]. However, before BlackHawk was released, there was no publicly available code to compute the Hawking radiation spectrum of a distribution of BHs in the universe, e.g. to derive precise constraints on the DM fraction under the form of PBHs or to predict detection signals of HR in future instruments; people usually relied on analytical approximations to compute the spectra. This was the motivation that drove the publication of the code. Moreover, since its first release, BlackHawk has received many modifications to embed interesting features: distributions of spinning BHs, emission of DM on top of the SM particles, HR rates in modified BH solutions, etc. Some of these modifications were already present in the previous version of the manual on arXiv [1905.04268v2], but most of them have been added in the recent update presented in the BlackHawk v2.0 release note [1]. The latter also contains a complete overview of the publications realised since BlackHawk v1.0 was released. The present manual is associated with BlackHawk v2.0, which is publicly available on the webpage:

http://blackhawk.hepforge.org/

It is structured as follows: in Section II we review the physics and main formulas of Hawking radiation, in Section III we provide the necessary compilation information, Section IV describes the parameters of a typical run of the code, Section V details the different routines, the two main BlackHawk programs are presented in Section VI, the output files are explained in Section VII, we briefly estimate the memory usage of a run in Section VIII, we outline some possible applications of the code in Section IX and we conclude in Section X.

II Hawking radiation

In this Section we review the physics of Hawking radiation. We present the formulas and results in their most general forms, and use the natural system of units $G=\hbar=k_{\rm B}=c=4\pi\varepsilon_{0}=1$ unless stated otherwise.

II.1 Black hole 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, such as the cosmic microwave background (CMB), extragalactic gamma-ray background (EGRB), the Big-Bang nucleosynthesis (BBN) yields, etc.

The distribution of BHs as a function of their mass or any other parameter (spin, charge) is completely model dependent. In BlackHawk, we are agnostic about the channel of formation of BHs and thus provide the possibility to use any distribution. For primordial BHs formed in the radiation era after inflation, there is typically a general relationship between the PBH formation time $t_{\rm f}$ and its initial mass [12, 30]

[TABLE]

where $\gamma\simeq 0.2$ encodes the details of the gravitational collapse of an overdense region. Of course, before Hawking radiation becomes relevant and depending on the plasma environment of PBHs, accretion (or mergers) could change this initial time-mass relationship [31]. Usually, a monochromatic (Dirac delta) mass/spin function is used to constrain the abundance of BHs. However, recent studies have shown that the generalisation to e.g. extended mass functions is not trivial [32]. The cleanest way consists in recomputing the constraints taking into account the extended distribution, rather than extrapolating from the monochromatic results. BlackHawk can in principle work with any distribution of BHs. Several BH mass and secondary BH parameter distributions are already built-in and listed in the following. It is also possible that the user provides its own distribution for BH masses and their secondary parameters thanks to a .txt file.

II.1.1 Mass distributions

Here we list the BlackHawk built-in mass distributions ${\rm d}n/{\rm d}M$ such that the integral

[TABLE]

gives the total number of comoving BHs per unit volume. To obtain the density of those BHs, one simply considers the integral

[TABLE]

In the following formulas, the normalization factor $A$ which appears for each distribution has a mass unit that depends on the details of the expression of ${\rm d}n/{\rm d}M$.

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. [33, 34])

[TABLE]

where the spectral index is $n\approx 1.3$ and $R_{\rm c}$ is measured using the CMB to be $R_{\rm c}=(24.0\pm 1.2)\times 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 [33] through peak-theory

[TABLE]

where

[TABLE]

in which:

• •

$H_{0}=67.8\,$km$\cdot$s${}^{-1}\cdot$Mpc*-1* is the Hubble parameter [35],

• •

$\Omega_{\rm m}=0.308$ is the matter energy density in the universe [35],

• •

$z_{\rm eq}=3200$ is the radiation-matter equality redshift [33],

• •

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

• •

$g_{*}=106.75$ is the number of (SM) relativistic energy dof at the time of BH formation, here the end of the inflation [35],

• •

$\zeta_{\rm th}=0.7$ is the overdensity threshold that provokes the direct collapse of a density fluctuation into a BH [33].

Log-normal distributions

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

[TABLE]

where $A$ is the amplitude, $M_{\rm c}$ is the position of the peak and $\sigma$ is its width. Note that this is a log-normal distribution for the comoving density $M{\rm d}n/{\rm d}M$ and not for the comoving number density ${\rm d}n/{\rm d}M$ – which differs only by a factor of $M$ and a different normalization $A$

[TABLE]

These distributions reduce to the monochromatic distribution in the limit $\sigma\rightarrow 0$, which can be used as a consistency check.

Power-law distribution

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

[TABLE]

where $\gamma\equiv-2w/(1+w)$ and $w$ is defined through the equation of state of the energy density that dominates the universe at the epoch of BH formation such that $P=w\rho$, see Table 1 below.

Critical collapse distribution

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

[TABLE]

where $A$ is an amplitude factor and $M_{\rm f}$ an upper cut-off.

Uniform distribution

The uniform distribution is a toy-model defined as

[TABLE]

where $A$ is an amplitude factor and $M_{\rm min,max}$ are the distribution bounds.

Dirac distribution

The Dirac distribution simulates a single-valued BH mass function $M=M_{\rm c}$

[TABLE]

It is useful to perform monochromatic analyses and consistency checks for a single BH. It is by default normalized to 1 BH per comoving cm3, any additional amplitude factor $A$ must be applied afterwards.

II.1.2 Secondary parameter distribution

Here we list the secondary BH parameter distributions. The no-hair theorem states that in general relativity, BHs are exclusively described by the set $\{M,a^{*},Q^{*}\}$ where $a^{*}\equiv a/M\equiv J/M^{2}$ is the adimensioned BH angular momentum and $Q^{*}\equiv Q/M$ is the adimensioned BH electric charge. The three parameters are constrained by

[TABLE]

otherwise the interior Cauchy horizon would collapse with the BH horizon and the singularity would be exposed, violating the GR singularity conjecture. Therefore, in GR we have two possible secondary BH parameters, namely the charge $Q^{*}$ and the angular momentum $a^{*}$. As is well known, GR can not be the fundamental description of gravity, as it is unable to provide a quantum mechanics framework for the latter. Lots of alternative models have been proposed, most of which include additional secondary parameters, e.g. the number of extra dimensions $n$ in higher-dimensional theories of gravity or the polymerization parameter $\varepsilon$ in polymerization models deriving from loop quantum gravity. Obviously, the two last examples of secondary parameters are constants common to all the BHs of a given model — the only physical distribution of those is a Dirac function. However, realistic models of BH formation should include an extended distribution of both angular momentum and electric charge, whose functional form depends on the BH formation mechanism. For example, most models predict that BHs form in electrically neutral environments, either from the collapse of stars or in the primordial universe. A phase transition in the early universe could however produce regions of non-zero electric charge that collapse to form charged PBHs. On the other hand, there are numerous models describing the spin distribution of BHs [36], which could originate from accretion [37], hierarchical mergers [38, 39], close encounters in BHs clusters [40, 41] or the initial angular momentum or non-sphericity of the collapsing region [42]. Interestingly, the spin of BHs could act as a distinguishing rule between primordial or stellar origin (e.g. [15, 43]), or increase the detection possibilities (e.g. [44, 45, 46]). Inside BlackHawk, we have built some benchmark distributions ${\rm d}\tilde{n}/{\rm d}x$ for both the angular momentum and the charge $x=a^{*},Q^{*}$, with an integral normalized to unity

[TABLE]

Gaussian distribution

The Gaussian distribution is a generalisation of the monochromatic distribution, which reduces to the latter in the limit $\sigma\rightarrow 0$

[TABLE]

where $\sigma$ is the standard deviation and $x_{\rm c}$ a characteristic value of $x$.

Uniform distribution

This is the same distribution as for the mass, but normalized to unity

[TABLE]

where $x_{\rm min,max}$ are the distribution bounds.

Dirac distribution

This is the same distribution as for the mass, and it also applies for all the possible single-valued secondary parameters $x=x_{\rm c}$

[TABLE]

II.2 Hawking evaporation physics

In this Section we review the basic formulas for the Hawking evaporation process. Hawking radiation is a semi-classical phenomenon arising from the mix between classical general relativity and aspects of quantum mechanics. Random quantum fluctuations at the horizon of BHs create pairs of particles and anti-particles. It is then possible that (anti-)particles escape to spatial infinity, giving rise to a net flux of outgoing radiation denoted as Hawking radiation. This process is qualitatively and quantitatively linked to the details of the BH solution metric; there are many such solutions of the Einstein equations of general relativity, not all of which exhibit Hawking radiation. For the spherically symmetric and static BHs, we rely on the two companion papers [47, 48] (see relevant extended bibliography there).

II.2.1 Spherically symmetric and static black holes

The general metric of a spherically symmetric and static BH solution is in Boyer–Lindquist coordinates

[TABLE]

where ${\rm d}\Omega^{2}\equiv{\rm d}\theta^{2}+\sin(\theta){\rm d}\varphi^{2}$ is the 4-dimensional 2-sphere element. We also require that this solution is asymptotically flat, meaning that

[TABLE]

and that $F(r)$ exhibits at least one pole at $r=r_{\rm H}$, the BH horizon. Examples used in BlackHawk are listed below.

Schwarzschild black hole

This is the simplest example of a BH solution, describing the spacetime curvature around a non-rotating, uncharged mass $M$ enclosed by a horizon at $r_{\rm H}=2M\equiv r_{\rm S}$. The metric coefficients are

[TABLE]

It is an example of the common $tr$-symmetric metrics with $F(r)=G(r)$ and $H(r)=r^{2}$. All other BH solutions embed this case as a limiting behaviour.

Charged black hole

Reissner–Nordström BHs describe the spacetime around a mass $M$ with charge $Q^{*}\equiv Q/M$. The $tr$-symmetric metric coefficients are

[TABLE]

where $r_{\rm Q}^{2}=Q^{2}$ in our system of units. The Hawking horizon $r_{\rm H}$ and Cauchy horizon $r_{\rm C}$ are the poles of $F(r)$, respectively

[TABLE]

We recover the Schwarzschild solution in the limit $Q^{*}\rightarrow 0$ (in that case, $r_{\rm Q}=0=r_{-}$). The opposite limit $Q^{*}\rightarrow 1$ is called “extremal” BH and $Q^{*}=1$ is not physical.

Higher-dimensional black hole

Models of general relativity with higher dimensions (more than 3 spatial dimensions) are not excluded yet by particle physics experiments or cosmological observations if the additional spatial dimensions are small enough and compactified [49]. The $tr$-symmetric metric coefficients are

[TABLE]

where $n=0,1,\dots$ is the number of extra spatial dimensions and

[TABLE]

is the BH horizon with $\Gamma$ the Euler gamma function. In these theories, the Planck mass is rescaled by $M_{\rm P}^{2}\sim M_{*}^{n+2}R^{n}$ where $R$ is the typical size of the extra dimensions. The 4-dimensional Schwarzschild case corresponds to $n=0$.

Polymerized black hole

Polymerized BHs are inspired by the quantization procedure in loop quantum gravity. They are part of a more general class of BH solutions that solve the central singularity problem, as the metric coefficients

[TABLE]

do not diverge at $r\rightarrow 0$. The horizon and Cauchy radii are

[TABLE]

where $\varepsilon$ is the polymerization parameter which describes the scale of the quantum corrections, $P(\varepsilon)$ is the polymerization function

[TABLE]

and the Arnowitt–Deser–Misner mass is given by

[TABLE]

The parameter $a_{0}$ is the area gap in loop quantum gravity. We recover the Schwarzschild case in the limits $\varepsilon,a_{0}\rightarrow 0$.

II.2.2 Rotating black holes

Realistic BH formation channels predict a non-zero angular momentum. The metric of a (rotating) Kerr BH is given by

[TABLE]

where $\Sigma(r)\equiv r^{2}+a^{2}\cos^{2}(\theta)$ and $\Delta(r)\equiv r^{2}-2Mr+a^{2}$. The horizon and Cauchy radii are

[TABLE]

It is not trivial at all to mix this BH solution with the (beyond GR) solutions listed above. Thus, in BlackHawk, we keep them separated and do not consider e.g. Kerr–Newman rotating and charged BHs. We recover the Schwarzschild case in the limit $a^{*}\rightarrow 0$, while the other limit $a^{*}\rightarrow 1$ is called “extremal” (the case $a^{*}=1$ is not physical).

II.2.3 Particle creation by black holes

Black hole temperature

Once the metric of spacetime around a BH is known, Hawking radiation rates are computed as follows. First, the deformation of spacetime is associated with an effective horizon temperature

[TABLE]

where $\kappa$ is the surface gravity of the BH. In the case of a spherically symmetric and static metric, it is given by [48]

[TABLE]

For a Kerr BH, it is

[TABLE]

We can point out that the temperature of an “extremal” BH with $a^{*}=1$, $Q^{*}=1$ or $\varepsilon\rightarrow+\infty$ is 0, which is thermodynamically not physical. The temperatures associated with the different metrics listed above are summarized in Table 2.

Emission rates

Then, the emission rate of some (beyond the) SM degrees of freedom $i$ per unit time and energy is given by the master equation [21]

[TABLE]

In this equation, we identify the Boltzmann statistics factor with the BH temperature in the denominator, with $\pm 1$ for fermions and bosons respectively. The particle energy $E$ has to be corrected for couplings between the emitted particle and the BH; in BlackHawk the only correction is brought by Kerr BHs for which the horizon rotation gives

[TABLE]

We do not consider inside BlackHawk the electric interaction between a Reissner–Nordström BH and a charged particle. In these equations, we have decomposed the wave function of the particle fields on the base of the spin-weighted spheroidal harmonics of spin $s_{i}$, and $l,m$ are respectively the particle angular momentum $l\geq s_{i}$ and its projection on the BH axis $m=-l,\dots,+l$. The last factor we need to determine is the greybody factor $\Gamma_{s_{i}lm}(E,M,x_{j})$. This factor encodes the probability that a (anti-)particle generated at the BH horizon escapes to spatial infinity. It depends on the particle energy $E$, the BH mass $M$, and the set of secondary parameters for the BH metric $x_{j}$ (e.g. the BH spin $a^{*}$). It depends a priori on the particle rest mass $\mu_{i}$, but for simplicity we apply the kinematic condition $E>\mu_{i}$ as a simple cut-off in the emission rates of massive particles. The computation of the greybody factors inside BlackHawk is detailed in Appendix C.1. In order to obtain the emission rate of some particle $i$, we sum over the field angular momentum $l,m$, helicity $g_{\rm helicity}^{i}$, color $g_{\rm color}^{i}$ and antiparticle $g_{\rm anti}^{i}$ multiplicities

[TABLE]

In the case of a spherically symmetric metric, the sum over $m$ reduces to a factor $2l+1$ as all angular momentum projections give the same contribution to the total emission rate. The fields we consider in BlackHawk are the SM fields (photon, gluons, Higgs boson, W± boson, Z0 boson, 3 neutrinos, electron, muon, tau, up, down, charm, strange, top, bottom quarks) plus possibly the graviton and some additional degree of freedom like a DM particle or dark radiation (DR). The multiplicities $g_{i}$, as well as the particles lifetimes and rest masses are summarized in Tables 3, 4 and 5 in Appendix B. To obtain the total primary spectrum of particle $i$ for an extended mass and/or secondary parameter $x$ distribution, we perform the integrals

[TABLE]

II.2.4 Black hole evolution

Once we know the emission rates of all the particles of the (beyond the) SM spectrum, we can determine the evolution of the BH driven by Hawking evaporation. To do that, we need to obtain the differential equations of evolution of the BH mass and the possible secondary parameter. The mass loss rate is [50, 51, 52]

[TABLE]

and the BH angular momentum loss rate is

[TABLE]

We could, on the same model, define the BH charge loss rate, but as we have not computed the emission of charged particles we can not get the full sum on the SM fields $i$. For higher-dimensional BHs, we have computed the emission rates on the 4-dimensional brane, but not in the bulk of the extra dimensions, where particles could acquire Kaluza-Klein piles of massive states. Thus, the sum is also incomplete in this case and we limit ourselves to the instantaneous emission. Finally, for polymerized BHs, only the factor $f(M,\varepsilon,a_{0})$ is relevant as the secondary parameters $\varepsilon,a_{0}$ are constant over the BH evolution. As a convention, we can take $g=0$ in this case. Then, we use Eqs. (38) and (39) to get the differential equations for the BH mass and spin [50, 51, 52]

[TABLE]

The computation of the $f(M,x_{j})$ and $g(M,a^{*})$ Page factors in BlackHawk is described in Appendix C.2.

II.2.5 Hadronization and decay

The elementary particles emitted by BHs are not the final products of Hawking emission. Some of them are unstable, others only exist in composite states (hadrons). A particle physics code has to be used in order to evolve the elementary particles into final products. We used PYTHIA [25], HERWIG [26] and Hazma [27] for this purpose.

The final particles, hereafter denoted as “secondary Hawking particles” (the elementary being the “primary Hawking particles”), depend on the cosmological context in which they are emitted. For BBN studies, an estimation of the reaction rates imposes to consider only the particles with a lifetime longer than $\sim\!\!10^{-8}\,$s, listed in Table 4 of Appendix B. For low redshift cosmological studies, one considers only the particles with an infinite lifetime, listed in Table 5. The emission rate of secondary particle $j$ emitted by a distribution of BHs per units of time and energy is computed with the integral

[TABLE]

where the sum is taken over Hawking primary particles $i$. To obtain the total secondary spectrum of particle $j$ for an extended distribution of mass and/or secondary parameter $x$ we perform the integrals

[TABLE]

Appendix C.3 describes how hadronization tables of the branching ratios ${\rm d}N_{j}^{i}(E^{\prime},E)$ have been computed to transform the primary spectra into secondary spectra inside BlackHawk.

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

III.1 Content

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 IV.1) and the numerical values of general quantities (unit conversion factors, constants, particle masses, etc),

• •

source files containing the definition of all the BlackHawk routines (evolution.c, general.c, hadro_herwig.c, hadro_pythia.c, hadro_pythianew.c, hadro_hazma.c, primary.c, secondary.c, spectrum.c, technical.c),

• •

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

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 [53, 54, 55]. These scripts can be used to generate modified numerical tables for example (see Appendices C.1, C.2 and C.3). They are accompanied by README.txt files explaining how to use them.

III.2 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 type333In 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.:

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 type444In case of memory problem at execution, increasing the stack size with the command ulimit -s unlimited can help solving the problem.:

./BlackHawk_*.x parameter_file

where parameter_file is the name of a parameter file (e.g. parameters.txt for the pre-built one). To compile only the library, just cd into the main directory and type:

make

IV Input parameters

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

IV.1 Parameter structure

The input parameters used by BlackHawk are listed in a parameter file (e.g. parameters.txt for the pre-built one). This file can be modified by the user and is saved for each new run of the code in the results/ directory. A C structure has been defined in include.h to embed all the parameters (i.e. input parameters and default run parameters):

struct param $\{$

char destination_folder[200];

int full_output;

int interpolation_method;

int BH_number;

double Mmin;

double Mmax;

int metric;

int param_number;

double amin;

double amax;

double Qmin;

double Qmax;

double epsilon_LQG;

double a0_LQG;

double n;

double M_star;

int spectrum_choice;

int spectrum_choice_param;

double amplitude;

double stand_dev;

double crit_mass;

double eqstate;

double stand_dev_param;

double mean_param;

char table[32];

int tmin_manual;

double tmin;

int limit;

int nb_fin_times;

int BH_remnant;

double M_remnant;

int E_number;

double Emin;

double Emax;

int particle_number;

int grav;

int add_DM;

double m_DM;

double spin_DM;

double dof_DM;

int primary_only;

int hadronization_choice;

double Mmin_fM;

double Mmax_fM;

int nb_fM_masses;

int nb_fM_param;

double param_min;

double param_max;

int nb_gamma_param;

int nb_gamma_x;

int nb_gamma_spins;

int nb_gamma_fits;

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 Section V 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 of them 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. Error messages are displayed if some parameter is out of its allowed range, automatically cancelling the run.

IV.2 General parameters

This 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 to debug the code or to follow the progress in time-consuming routines. It also determines whether some interactive pre-run checks are done with the user. We thus recommend that this parameter is always set to 1.

• •

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

IV.3 BH spectrum parameters

This set of parameters defines the quantities used to compute the BH density distribution:

• •

BH_number is the number of BH masses that will be simulated. If the parameter spectrum_choice is not set to $-1$, 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 $-1$, it has to be the number of tabulated values in the user-defined BH distribution. It will be automatically set to 1 if spectrum_choice is set to 0 (Dirac distribution).

• •

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_{\rm p}\approx 2\times 10^{-5}\,$g $<$ Mmin, Mmax, where $M_{\rm p}$ is the Planck mass. For an extended mass distribution, one must have Mmin $<$ Mmax. If they are not compatible with boundaries of the mass distribution, the computation will stop.

• •

metric switches between the different BH solutions available in BlackHawk: 0 (Kerr BH), 1 (polymerized BH), 2 (Reissner–Nordström BH), 3 (higher-dimensional BH).

• •

param_number is the number of BH secondary parameters that will be simulated. It is relevant only if metric is 0 or 2, otherwise there is only one possible value of the secondary parameters. If the parameter spectrum_choice is not set to $-1$, it has to be an integer greater than or equal to 1. If it is equal to 1, the only BH secondary parameter will be amin (Kerr BHs), Qmin (Reissner–Nordström BHs) or the single-valued secondary parameter in the case of polymerized or higher-dimensional BHs. If the parameter spectrum_choice is set to $-1$, it has to be the number of tabulated values in the user-defined BH distribution. It will be automatically set to 1 if spectrum_choice_param is set to 0 or if metric is set to 1 or 2.

• •

amin and amax are respectively the lowest and highest BH dimensionless reduced spins that will be simulated. They have to satisfy the condition $0\leq$ amin, amax $<1$. For an extended spin distribution, one must have amin $<$ amax. If they are not compatible with boundaries of the spin distribution, the computation will stop. The paradigm is the same with Qmin and Qmax for the dimensionless BH electric charge.

• •

epsilon_LQG is the dimensionless quantization parameter $\varepsilon$ for polymerized BHs. It has to be a positive number $0\leq\varepsilon\leq 100$; a0_LQG is the area gap $a_{0}$ in loop quantum gravity, it has to be a positive number, but only the two values 0 and $0.11$ have been used to compute the numerical tables.

• •

n is the number of extra spatial dimensions $n$ for higher-dimensional BHs, it has to be an integer $0\leq n\leq 6$; M_star is the rescaled Planck mass in Eq. (24), only the value $M_{*}=1$ has been used in the numerical tables.

• •

spectrum_choice selects the form of the BH mass distribution, it has to be an integer among 0 (Dirac), 1 (log-normal for the mass density), 11 (log-normal for the number density), 2 (power-law), 3 (critical collapse), 4 (peak theory), 5 (uniform) and $-1$ (user defined distribution, provided as a .txt file in the subfolder /src/tables/users_spectra/).

• •

spectrum_choice_param selects the form of the BH secondary parameter distribution. It has to be an integer among 0 (Dirac), 1 (uniform) and 2 (Gaussian). If the parameter spectrum_choice has been set to $-1$, then the spin distribution is also the user-defined one and this parameter is irrelevant. If the BH metric is set to polymerized or higher-dimensional BHs, then the secondary parameter is automatically single-valued.

• •

amplitude is the amplitude $A$ present in Eqs. (7), (8), (9), (10) and (11). It is the normalization of the corresponding BH distribution and thus strictly positive, its unit depends on the distribution chosen, but masses are in grams and densities in centimetres*-3*.

• •

stand_dev is the dimensionless standard deviation $\sigma$ in the log-normal distributions of Eqs. (7) and (8). It has to be strictly positive.

• •

crit_mass is the characteristic mass $M_{\rm c}$ in Eqs. (7) and (8) and $M_{\rm f}$ in Eq. (10). It has to be strictly positive and given in grams.

• •

eq_state defines the equation of state $w$ for the power-law mass distribution of Eq. (9), it is dimensionless.

• •

stand_dev_param is the dimensionless standard deviation $\sigma$ in the Gaussian distribution (15) of the secondary parameter; it has to be strictly positive.

• •

mean_param is the characteristic dimensionless secondary parameter $x_{\rm c}$ in the Gaussian distribution (15) or in the uniform distribution (16), it has to be strictly positive.

• •

table is the name of a user-defined BH distribution table. It has to be a string with any file extension. This file has to have the same format as the BlackHawk generated files, that is to say a two-entry table with lines corresponding to the BH mass and rows to the BH secondary parameters.

If the parameters spectrum_choice and spectrum_choice_param are both set to 0, then the distribution is monochromatic in both mass and secondary parameter, mimicking a single BH — that is to say, the emissivities obtained are those of a single BH. This option can be useful to compute and compare known test emissivities of single BHs.

IV.4 BH evolution parameters

This set of parameters defines the quantities used to compute the BH evolution:

• •

tmin_manual switches between automatically set $t_{\rm min}$ for the formation time of BHs (tmin_manual = 0) following Eq. (1) and manually set $t_{\rm min}$ (tmin_manual = 1). For now, if the BH distribution is extended in mass, then all BHs are formed at the same time — meaning, the formation time is common to all the initial BH masses. This should have negligible impact on the emission rates.

• •

tmin is the initial integration time of the evolution of BH, in seconds, if tmin_manual has been set to 1. It can have any positive value.

• •

limit is the iteration limit when computing the time evolution of a single BH. It is fixed to limit = 5000 even if the effective iteration number hardly reaches 1000. It should be increased if the integration does not reach the complete evaporation of BHs or if the following error appears:

[life_evolution] : ERROR ITERATION LIMIT REACHED !

• •

nb_fin_times is the number of final integration times of the computations. It is set automatically by the integration procedure.

• •

BH_remnant switches between two cases: either the BH evaporation continues until the Planck mass is reached, and at this moment we are agnostic about what happens but we believe that Hawking radiation becomes irrelevant (BH_remnant = 0); or the BH evaporation stops before reaching the Planck mass, leaving a massive remnant (BH_remnant = 1).

• •

M_remnant is the BH remnant mass in grams if BH_remnant = 1.

IV.5 Primary spectrum parameters

This set of parameters defines the quantities related to the primary Hawking spectra:

• •

E_number is the number of primary particle energies that will be simulated. It has to be an integer greater than or equal to 2.

• •

Emin and Emax are the minimum and maximum primary particle energies, respectively. They must be compatible with the numerical tables boundaries and satisfy $0<\texttt{Emin}<\texttt{Emax}$.

• •

particle_number is the number of primary particle types. It is fixed to 15 (photon, gluon, W± boson, Z0 boson, Higgs boson, neutrino, 3 leptons (electron, muon, tau) and 6 quarks (up, down, charm, strange, top, bottom)), but is increased in the code if graviton ($+1$) or DM emission ($+1$) is activated.

• •

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

• •

add_DM swithces between DM emission (add_DM = 1) or not (add_DM = 0). It is automatically set to 1 if hadronization_choice is set to 3 (Hazma tables), as in this case the DM slot will be used to compute the emission of primary pions, thus forbidding a mixed scenario “DM + Hazma hadronization”.

• •

m_DM, spin_DM and dof_DM are respectively the additional DM mass in GeV, DM spin (among 0, 1, 2 and $1/2$, while the value $3/2$ is only available for Schwarzschild BHs as these are the only spin $3/2$ greybody factors we have computed) and DM number of dof (any positive integer). The evolution numerical tables $f(M,x_{j})$ and $g(M,a^{*})$ have been computed for dof_DM = 1 only, but the change should be negligible if dof_DM remains much smaller than the number of dof available at some BH mass $M$ ($g=106.75$ at $M\ll 10^{10}\,$g for the full SM, $g=7.25$ at $M\gg 10^{19}\,$g for the photon and 3 massless Majorana neutrinos); the user can recompute these tables anyway if needed.

IV.6 Hadronization parameters

This set of parameters defines the quantities used during hadronization:

• •

primary_only determines whether the secondary spectra will be computed or not. It has to be an integer between 0 (primary and secondary spectra) and 1 (primary spectra only). In the case where the parameters Emin and Emax are not compatible with the hadronization tables boundaries, a warning will be displayed and extrapolation used (hadronization_choice = 0,1,2), or the computation will stop altogether (hadronization_choice = 3).

• •

hadronization_choice determines which hadronization tables will be used to compute the secondary spectra. It has to be an integer between 0 (PYTHIA tables, early Universe/BBN epoch), 1 (HERWIG tables, early Universe/BBN epoch), 2 (PYTHIA tables, present epoch) and 3 (Hazma tables, present epoch).

IV.7 Tables parameters

This set of parameters is quite special: it regroups the numerical table parameters, which are fixed when the tables are computed and should not be modified by the user. Thus, they are set in separated infos.txt files. Their list is:

• •

Mmin_fM, Mmax_fM, nb_fM_masses and nb_fM_param describe the format of the evolution tables $f(M,x_{j})$ and $g(M,a^{*})$ (see Appendix C.2). They are read in the file:

src/tables/fM_tables/infos.txt

• •

param_min, param_max, nb_gamma_param, nb_gamma_x, nb_gamma_spins and nb_gamma_fits describe the format of the greybody factor tables (see Appendix C.1). They are read in the file:

src/tables/gamma_tables/infos.txt

• •

Emin_hadro, Emax_hadro, nb_init_en, nb_fin_en, nb_init_part and nb_fin_part describe the format of the hadronization tables (see Appendix C.3). They are read in the file:

src/tables/hadronization_tables/infos.txt

V Routines

In this Section, we have listed all the routines defined in BlackHawk. To simplify the analytic formulas, all intermediate quantities are in GeV (see Appendix A for conversion rules).

V.1 General routines

There are seven general routines in the BlackHawk code. The principal ones are the two main routines, described in Section VI. The other five are:

• •

int read_params(struct param *parameters, char name[], int session)

This routine reads the file name and the numerical tables information files thanks to the read_*_infos routines. 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 and are skipped automatically. It takes a pointer to a struct param object (see Section IV.1) as an argument and fills it using the file name. The argument session keeps track of which one of the main programs has been launched (0 for BlackHawk_tot, 1 for BlackHawk_inst). If one parameter is not of the type described in Section IV this function will display an error message of the kind:

[read_params] : ERROR ... !

sometimes accompanied with an explanation message. Any of these errors will make the routine return the value 0 and end the BlackHawk run. If one parameter is in 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, of the kind:

[read_params] : WARNING ... !

In such a case, the problematic parameters will be set automatically (e.g. primary_only = 1) and the computation will be performed.

• •

int read_fM_infos(struct param *parameters)

This routine reads the $f(M,x_{j})$ and $g(M,a^{*})$ table information in the file:

src/tables/fM_tables/infos.txt

It returns 1 if all parameters are correct and 0 if there is an error, displaying an error message in this case:

[read_fM_infos] : ERROR ... !

• •

int read_gamma_infos(struct param *parameters)

This routine reads the greybody factor table information in the file:

src/tables/gamma_tables/infos.txt

It returns 1 if all parameters are correct and 0 if there is an error, displaying a message in this case:

[read_gamma_infos] : ERROR ... !

• •

int read_hadronization_infos(struct param *parameters)

This routine reads the hadronization table information in the file:

src/tables/hadronization_tables/infos.txt

It returns 1 if all parameters are correct and 0 if there is an error, displaying a message in this case:

[read_hadronization_infos] : ERROR ... !

• •

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 VIII). If the parameter full_output is set to 1, the user has the choice to cancel the run and the value 0 is returned, otherwise it is 1. The output is given in MB.

V.2 BH initial distribution routines

There are six routines contributing to the BH initial distribution computation:

• •

void read_users_table(double *init_masses, double *init_params, double **spec_table, struct

param *parameters)

This routine reads a user-defined BH distribution table in the file given by the parameter table, if the parameter spectrum_choice is set to $-1$. It fills the arrays init_masses[], init_params[] 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. (6).

• •

double M_dist(double M, struct param *parameters)

This routine takes a BH mass as an argument and computes the comoving density ${\rm d}n/{\rm d}M$ defined in Eq. (5) (using the nu routine) or Eqs. (7), (8), (9), (10) or (11) (in GeV2 $\rightarrow$ cm${}^{-3}\cdot$g*-1*), depending on the parameter spectrum_choice. If this parameter is set to 0, a flat distribution is used with only one BH mass.

• •

double param_dist(double param, struct param *parameters)

This routine takes a BH secondary parameter as an argument ($x=a^{*}$ or $x=Q^{*}$) and computes the fraction ${\rm d}\tilde{n}/{\rm d}x$ defined in Eqs. (15) or (16) (dimensionless), depending on the parameter spectrum_choice_param. If this parameter is set to 0, a Dirac distribution is used with only one $x$ value.

• •

void spectrum(double *init_masses, double *init_spins, double **spec_table, struct param

*parameters)

This routine fills the array init_masses[] with BH_number BH masses logarithmically distributed between Mmin and Mmax. If the parameter BH_number is set to 1, the only BH initial mass will be Mmin. It then fills the array init_params[] with secondary parameters linearly distributed between $a_{\rm min},a_{\rm max}$ or $Q_{\rm min},Q_{\rm max}$ (for Kerr and Reissner–Norström BHs respectively), and finally fills the array spec_tables[][] computing the corresponding comoving densities ${\rm d}n(M,x)$ (in GeV3 $\rightarrow$ cm*-3*) using the M_dist and param_dist routines where ${\rm d}M$ is taken around the considered mass and ${\rm d}x$ around the considered secondary parameter. The result is rescaled by a factor $10^{100}$ due to the very small numbers involved.

• •

void write_spectrum(double *init_masses, double *init_spins, double **spec_table, struct param *parameters)

This routine writes the BH initial masses, secondary parameters and comoving densities in the file:

results/destination_folder/BH_spectrum.txt

detailed in Section VII. The results are converted from GeV to CGS units.

V.3 BH evolution routines

There are nineteen routines contributing to the BH time evolution computation:

• •

double rplus_Kerr(double M, double a)

This routine gives the external Kerr radius $r_{+}\equiv r_{\rm H}$ of Eq. (30) of a rotating BH for a given mass $M$ and reduced spin $a^{*}$ (in GeV*-1* $\rightarrow$ cm).

• •

double temp_Kerr(double M, double a)

This routine gives the Hawking temperature of a Kerr BH for a given mass $M$ and reduced spin $a^{*}$, see Table 2 (in GeV $\rightarrow$ K).

• •

double P_LQG(double epsilon)

This routine computes the dimensionless polymerization function $P(\varepsilon)$ defined in Eq. (27).

• •

double m_LQG(double M, double epsilon)

This routine computes the effective mass $m(\varepsilon)$ of a polymerized BH with ADM mass $M$, defined in Eq. (28) (in GeV $\rightarrow$ g).

• •

double temp_LQG(double M, double epsilon, double a0)

This routine gives the temperature of a polymerized BH of mass $M$, polymerization factor $\varepsilon$ and area gap $a_{0}$, see Table 2 (in GeV $\rightarrow$ K).

• •

double rplus_charged(double M, double Q) and rminus_charged(double M, double Q)

These routines compute the external horizon and Cauchy radii $r_{\rm H}\equiv r_{+}$ and $r_{\rm C}\equiv r_{-}$ of a Reissner–Nordström BH of mass $M$ and adimensioned charge $Q^{*}$, defined in Eq. (22) (in GeV*-1* $\rightarrow$ cm).

• •

double temp_charged(double M, double Q)

This routine gives the temperature of a Reissner–Nordström BH of mass $M$ and adimensioned charge $Q^{*}$, see Table 2 (in GeV $\rightarrow$ K).

• •

double rH_higher(double M, double n, double M_star)

This routine computes the horizon radius of a higher-dimensional BH of mass $M$ and extra dimensions number $n$ with rescaled Planck mass $M_{*}$, defined in Eq. (24) (in GeV*-1* $\rightarrow$ cm).

• •

double temp_higher(double M, double n, double M_star)

This routine gives the temperature of a higher-dimensional BH of mass $M$ and extra dimensions number $n$ with rescaled Planck mass $M_{*}$, see Table 2 (in GeV $\rightarrow$ K).

• •

void read_fM_table(double **fM_table, double *fM_masses, double *fM_param, struct param

*parameters)

This routine reads the $f(M,x_{j})$ factors defined in Eq. (38) in the relevant table contained in the folder:

src/tables/fM_tables/

depending on the parameter metric. It fills the arrays fM_masses[] (in GeV $\rightarrow$ g), fM_param[] (dimensionless) and fM_table[][] with format [mass][param] (in GeV4 $\rightarrow$ g${}^{3}\cdot$s*-1*).

• •

void read_gM_table(double **gM_table, double *fM_masses, double *fM_param, struct param

*parameters)

This routine reads the $g(M,a^{*})$ factor defined in Eq. (39) in the relevant table contained in the folder:

src/tables/fM_tables/

It fills the arrays fM_masses[] (in GeV $\rightarrow$ g), fM_param[] and gM_table[][] with format [mass][param] (in GeV4 $\rightarrow$ g${}^{2}\cdot$GeV$\cdot$s*-1*).

• •

double loss_rate_M(double M, double param, double **fM_table, double *fM_masses, double

*fM_param, int counter_M, int counter_param, struct param *parameters)

This routine computes the quantity ${\rm d}M/{\rm d}t$ defined in Eq. (40a) (in GeV2 $\rightarrow$ g$\cdot$s*-1*).

• •

double loss_rate_param(double M, double param, double **fM_table, double **gM_table, double

*fM_masses, double *fM_param, int counter_M, int counter_param, struct param *parameters)

This routine computes the quantity ${\rm d}a^{*}/{\rm d}t$ defined in Eq. (40b) (in GeV $\rightarrow$ s*-1*).

• •

evolution_times(double *init_masses, double *init_params, double **evol_times, double

**fM_table, double **gM_table, double *fM_masses, double *fM_a, struct param *parameters)

This routine computes the lifetime of all BHs with masses in init_masses[] and secondary parameters in init_params[] and stores the results in evol_times[][] with format [mass][param] (in GeV*-1* $\rightarrow$ s).

• •

sort_lifetimes(double **evol_times, double *sorted_times, int **rank, struct param

*parameters)

This routine sorts the lifetimes of the BHs computed by evolution_times from shortest to longest using the routine sort_fusion and stores the corresponding ranks in the table rank[][] with format [mass][param].

• •

void life_evolution(double ***life_masses, double ***life_params, double *life_times, double *dts, double *init_masses, double *init_params, int **rank, double **fM_table, double

**gM_table, double *fM_masses, double *fM_param, struct param *parameters)

This routine computes the evolution of each of the initial BH masses in init_masses[] and BH secondary parameters in init_params[]. The initial time life_times[0] is set to tmin, the initial masses life_masses[i][j][0] are set to init_masses[i] and the initial secondary parameters life_params[i][j][0] are set to init_params[j]. Iteratively, the next masses and secondary parameters are estimated using the simple Euler method

[TABLE]

where the derivatives are computed using the loss_rate_* routines. If the mass or the secondary parameter relative variations of the currently interesting BH are too large ($|{\rm d}X/X|>0.1$) then the time interval is divided by 2. If all the variations are very small ($|{\rm d}X/X|<0.001$), and if the current timestep is reasonable compared to the current timescale (${\rm d}t/t\lesssim 1$) then the time interval is multiplied by 2. If the dimensionless spin (Kerr BHs) reaches $10^{-3}$, we stop computing its variation and simply set it to 0, and it does not enter any more in the adaptive timesteps conditions. This goes on until each BH reaches the Planck mass or the remnant mass. As BH evolution is very steep towards the end of evaporation, the BH that determines the time interval is determined by following the order given by the array rank[][] which contains the information about the chronological order in which the BHs will evaporate. If the recursion limit limit $\times$ BH_number is reached, 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 ${\rm d}t$, times $t$, masses $M$ and secondary parameters $x$ are stored in the arrays dts[], life_times[] (both in GeV*-1* $\rightarrow$ s), life_masses[][][] (in GeV $\rightarrow$ g) and life_params[][][] (dimensionless), respectively, with format [mass][param][time] for the two last ones.

• •

void write_life_evolutions(double ***life_masses, double ***life_params, double *life_times, struct param *parameters)

This routine writes the BH time-dependent masses and secondary parameters until the end of evaporation in the file:

results/destination_folder/life_evolutions.txt

detailed in Section VII. If the BH distribution is extended in both mass and secondary parameter, we advise the user to deactivate this writing because of the extensive memory space used in this case (see Section VIII). The results are converted from GeV to CGS units.

V.4 Primary spectra routines

There are five routines contributing to the computation of the primary Hawking spectra:

• •

void read_gamma_tables(double ***gammas, double *gamma_param, double *gamma_x, struct param

*parameters)

This routine reads the greybody factors $\Gamma_{s_{i}lm}(E,M,x_{j})/(e^{E^{\prime}/T}-(-1)^{2s_{i}})$ defined in Eq. (34), in the files:

src/tables/gamma_tables/[/]spin_.txt

It fills the arrays gamma_param[] and gamma_x[] with the tabulated secondary parameters (dimensionless) and $x\equiv Er_{\rm S}$ (dimensionless $\rightarrow$ GeV$\cdot$cm), respectively. It fills the array gammas[][][] with the corresponding dimensionless greybody factors in format [type][param][$x$].

• •

void read_asymp_fits(double ***fits, struct param *parameters)

This routine reads the asymptotic fit parameters for the greybody factors, contained in the files:

src/tables/gamma_tables/[/]spin__fits.txt

It fills the array fits[][][] with format [type][param][coefs].

• •

double dNdtdE(double E, double M, double param, int particle_index, double ***gammas, double *gamma_param, double *gamma_x, double ***fits, double *dof, double *spins, double

*masses_primary, int counter_param, int counter_x, struct param *parameters)

This routine computes the emission rate ${\rm d}^{2}N_{i}/{\rm d}t{\rm d}E$, defined in Eq. (36) of the primary particle particle_index, for a given particle energy $E$, BH mass $M$, secondary parameter and with the primary particle information contained in the arrays dof[], spins[] and masses_primary[]. If $x\equiv Er_{\rm S}$ is inside the greybody factor tables boundaries, the values are interpolated in these tables at position counter_param and counter_x. Otherwise, it uses the asymptotic high and low energy fit tables. If the energy condition $E>\mu$ is broken, where $\mu$ is the particle rest mass, it returns 0. The result is dimensionless ($\rightarrow$ GeV${}^{-1}\cdot$s*-1*).

• •

void instantaneous_primary_spectrum(double **instantaneous_primary_spectra, double *BH_masses, double *BH_params, double **spec_table, double *energies, double ***gammas, double

*gamma_param, 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 quantities ${\rm d}^{2}n_{i}/{\rm d}t{\rm d}E$ defined in Eq. (37) for each primary particle $i$ and each primary energy in energies[], thanks to the routine dNdtdE. The results are stored in the array instantaneous_primary_spectra[][] with 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 the file:

results/destination_folder/instantaneous_primary_spectra.txt

detailed in Section VII. The results are converted from GeV to CGS units.

V.5 Secondary spectrum routines

There are thirteen routines contributing to the computation of the secondary Hawking spectra:

• •

void convert_hadronization_tables(double ****tables, double *initial_energies, double

*final_energies, struct param *parameters)

This routine is auxiliary. It writes hardcoded versions of the hadronization tables in files:

src/tables/hadronization_tables/hadronization_tables_*.h

in order to accelerate the code execution (while slightly slowing its compilation).

• •

void read_hadronization_tables(double ****tables, double *initial_energies, double

*final_energies, struct param *parameters)

This routine reads the relevant hadronization tables in the files:

src/tables/hadronization_tables/_tables/.txt

depending on the parameter hadronization_choice. If HARDTABLES is defined, it uses the hardcoded tables included at compilation thanks to the routines read_hadronization_, which is faster. 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 branching ratios ${\rm d}N_{j}^{i}/{\rm d}E$ defined in Eq. (41) (in GeV-1*) with format [secondary particle][initial energy][final energy][primary particle].

• •

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 timestep in times and writes it directly in the output in order to save RAM memory. To do so, it creates the output files (one per particle)

results/destination_folder/*_primary_spectrum.txt

and

results/destination_folder/*_secondary_spectrum.txt

if primary_only is set to 0. It reads the writing instructions using read_writing_instructions. Then, it fills the partial arrays partial_* with the instantaneous primary spectra, hadronized spectra and integrated spectra at each intermediate time and calls the routine write_lines to write the partial result in the output before moving to the next timestep.

• •

void read_writing_instructions(int *write_primary, int *write_secondary, struct param

*parameters)

This routine reads the writing information contained in the files

src/tables/write_*.txt

depending on the hadronization_choice, and stores it into the arrays write_primary[] and write_secondary[].

• •

void write_lines(char **file_names, double **partial_integrated_hadronized_spectra, int

*write_primary, int *write_secondary, double time, struct param *parameters)

This routine writes, given a time $t$ and instantaneous primary and secondary spectra (if primary_only is set to 0), a new line in the files:

results/destination_folder/*_primary_spectrum.txt

and

results/destination_folder/*_secondary_spectrum.txt

detailed in Section VII. The arrays write_*[] determine whether the values for each particle are written or not, thus potentially saving disk memory. Results are converted from GeV to CGS units.

• •

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 one part $i$ of the instantaneous integrand of Eq. (41) (in GeV*-1* $\rightarrow$ GeV${}^{-2}\cdot$s*-1*) for the secondary particle particle_type, initial energy $E^{\prime}=$ 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, that is to say the integrand of Eq. (41) 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, with format [secondary particle][initial energy][final energy]. If the initial energy is not in the hadronization tables, the contribution is extrapolated (this is not the case for hadronization_choice = 3).

• •

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 of Eq. (42) (dimensionless $\rightarrow$ GeV${}^{-1}\cdot$s*-1*) using the trapeze routine. The results are stored in the array instantaneous_integrated_hadronized_spectra[][] with 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 three routines add the contributions of the primary photons, neutrinos and electrons to the secondary produced ones. These (trivial) branching ratios are not included in the computation of the hadronization tables with PYTHIA or HERWIG. The value in term of final energies is interpolated in the primary spectrum and added to the hadronized spectrum instantaneous_integrated_hadronized_spectra[][].

• •

void add_FSR_instantaneous(double **instantaneous_primary_spectra, double

**instantaneous_integrated_hadronized_spectra, double *energies, double *final_energies,

double *masses_primary, struct param *parameters)

This routine adds the final state radiation from charged primary particles (electrons, muons, pions) to the secondary photon spectrum in the case of hadronization_choice = 3, see Appendix C.3.

• •

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:

results/destination_folder/instantaneous_secondary_spectra.txt

detailed in Section VII. The results are converted from GeV to CGS units.

V.6 Auxiliary routines

Fourteen 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 through

[TABLE]

• •

void free1D_double(double *array), void free1D_int(int *array), void free2D_int(int **array, int l_1stD), void free2D_double(double **array, int l_1stD), void free2D_char(char **array, int l_1stD), void free3D_double(double ***array, int l_1stD, int l_2ndD), void free3D_int(int

***array, int l_1stD, int l_2ndD) and void free4D_double(double

****array, int l_1stD, int l_2ndD, int l_3rdD)

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.

• •

void fusion(double *table, int start1, int end1, int end2), void sort_fusion_bis(double *table, int start, int end) and void sort_fusion(double *table, int llength)

These 3 routines perform a fusion sorting of the table table of length llength.

V.7 General feature

We have defined arrays compute[] that contain either 0’s or 1’s for all the primary and secondary particles and decide whether the contribution of the corresponding particle will be taken into account in the computation of the Hawking spectra. For example, if a “0” is set for the primary photon computation, then its primary spectrum will be all 0’s and thus will not contribute to the secondary spectra. If a “0” is set for the secondary photon computation, then the secondary particles generated by primary photons are not taken into account in the secondary spectra (but the primary spectra of photons is still computed). These arrays are by default all 1’s and should be modified only if the user wants to test specific scenarios or correct a bug.

VI BlackHawk 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. We will now describe the structure of the main routines together with screen output examples.

VI.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^{-10}$ defines the precision up to which two double numbers are considered as equal,

• •

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

• •

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

• •

m_* are the masses of the SM elementary and composite particles (see Tables 3, 4 and 5 in Appendix B),

• •

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

The code further works 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, more information will be displayed about the progress of the steps and an interactive pre-run check will be performed (overwrite existing data? check for memory usage?). In the case where information appears with the name of another routine inside brackets, it means that an error occurred, or that a warning must be displayed due to an unusual use of the code.

The first common step is the definition and filling of the parameters structure using read_params. If full_output is set to 1, 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. If no error is found in the input parameters, the output directory:

results/destination_folder/

is created. If it already exists and if full_output is set to 1, the user has the choice to overwrite the existing data or to stop the execution in order to choose another output folder, otherwise data will be overwritten automatically. For a subsequent data interpretation, the parameters file is copied in the output folder:

results/destination_folder/parameters.txt

The expected output at this stage if of the form (full_output = 1, program BlackHawk_inst):

             ############################

             #          BLACKHAWK v2.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 121.447 MB of RAM and 0.068 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. In the following, screen output examples are given in the mode full_output = 0. The results associated with the built-in parameters.txt file are given in Appendix D and should be used as checks at the first BlackHawk use.

VI.2 BlackHawk_tot: Full time-dependent Hawking spectra

In this program, BlackHawk computes the time-dependent Hawking spectra of a chosen initial distribution of BHs, that is the quantities ${\rm d}^{2}n_{i}/{\rm d}t{\rm d}E$ at all time intervals $t$ from the initial time $t_{\rm min}$ to the end of evaporation.

BlackHawk computes the initial distribution of BHs at $t_{\rm min}$ using the routine spectrum or reads the user-defined BH distribution file table with the routine read_users_table (depending on the spectrum_choice), filling the arrays init_masses[], init_params[] and spec_table[][]. It writes the results in the output with write_spectrum.

It then reads the relevant $f(M,x_{j})$ (and $g(M,a^{*})$ for Kerr BHs) tables, depending on the metric choice, using the read_fM_table (and read_gM_table) routines, filling the arrays fM_table[][], gM_table[][], fM_masses[] and fM_param[], in order to evolve in time each initial BH mass and secondary parameter down to the Planck or remnant mass using the routine life_evolution. This fills the arrays life_times[], life_masses[][][], life_params[][][] and dts[]. The evolutions in time are written in the output using the routine write_life_evolutions.

Then BlackHawk reads the relevant greybody factor tables with the read_gamma_tables routine, depending on the metric choice, filling the arrays gammas[][][], gamma_param[] and gamma_x[], and the fit tables using read_asymp_fits, filling the array fits[][][].

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 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 after 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 using write_lines.

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

VI.3 BlackHawk_inst: Instantaneous Hawking spectra

In this program, BlackHawk computes the instantaneous Hawking spectra of a distribution of BHs. That is, the quantities ${\rm d}^{2}n_{i}/{\rm d}t{\rm d}E$ at the initial masses and secondary parameters, without time evolution.

First BlackHawk computes the initial distribution of BHs using the routine spectrum or it reads the user-defined BH distribution file table with the routine read_users_table (depending on the spectrum_choice), filling the arrays init_masses[], init_params[] and spec_table[][]. It writes the results in the output with write_spectrum.

Then BlackHawk reads the relevant greybody factor tables (depending on the metric choice) using the routine read_gamma_tables, filling the arrays gammas[][][], gamma_masses[] and gamma_x[] and the fit table with the routine read_asymp_fits, filling the array fits[][][]. It computes the primary Hawking spectra using the routine instantaneous_primary_spectrum, filling the arrays instantaneous_primary_spectra[][]. The results are written in the output by write_instantaneous_primary_spectra.

If the parameter primary_only has been set to 0, BlackHawk reads the relevant 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 integrate_initial_energies_instantaneous, which fills the array instantaneous_integrated_hadronized_spectra[][]. The contributions from primary photons, neutrinos and electrons 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.

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 FIT 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

VII Output files

All the output files generated by a run of BlackHawk will be stored in the folder:

results/destination_folder/

In this Section we describe the format of the files created by each program. Examples of results can be found in Appendix D. In all 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 provided in the folder:

scripts/visualization_scripts/*.py

in order to plot the data produced by both programs. They come with a README that explains how to configure them. The user can of course modify these scripts or use any other plotting program.

VII.1 BlackHawk_tot program

Running BlackHawk_tot produces five (or four if primary_only is set to 1) 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 is 2-dimensional: the first column is a list of the BH initial masses (in g), the first line a list of the initial secondary parameters (dimensionless) and the bulk 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 and secondary parameter. It includes the total number of integration timesteps in the description line. It also contains BH_number tables in which the first column gives the time (in s), and each other pair of columns is the evolution of the mass (in g) and secondary parameter (dimensionless) of BHs as a function of time for a fixed initial mass.

• •

dts.txt: this file is also written by write_life_evolutions. It contains the integrated timesteps $t$ (first column) and the corresponding time intervals ${\rm d}t$ (second column), both in s.

• •

_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 ${\rm d}^{2}n_{i}/{\rm d}t{\rm d}E$ of Eq. (37) of the particle per unit energy, time and covolume (in GeV${}^{-1}\cdot$s${}^{-1}\cdot$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 ${\rm d}^{2}n_{j}/{\rm d}t{\rm d}E$ of Eq. (42) of the particle per units of energy, time and covolume (in GeV${}^{-1}\cdot$s${}^{-1}\cdot$cm-3*). These files will not be generated if the parameter primary_only has been set to 1.

VII.2 BlackHawk_inst program

Running BlackHawk_inst produces three (or two if primary_only is set to 1) output files:

• •

BH_spectrum.txt: this file is the same as for BlackHawk_tot.

• •

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 ${\rm d}^{2}n_{i}/{\rm d}t{\rm d}E$ of Eq. (37) per unit energy and time (in GeV${}^{-1}\cdot$s${}^{-1}\cdot$cm*-3*).

• •

instantaneous_secondary_spectra.txt: this file is written by 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 ${\rm d}^{2}n_{j}/{\rm d}t{\rm d}E$ of Eq. (42) per unit energy and time (in GeV${}^{-1}\cdot$s${}^{-1}\cdot$cm*-3*). It will not be generated if the parameter primary_only has been set to 1.

VIII Memory usage

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.

VIII.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 subsequent part of the run. Then, the RAM used by BlackHawk at a given step of a session can be estimated as a sum over all active arrays at that time. double numbers are coded in 8 bytes and int in 4 bytes. Memory spaces $\mathcal{M}$ are given in bytes. We assume that $\langle N\rangle\sim 1000$ timesteps are necessary for the integration of one BH mass and secondary parameter. For BlackHawk_tot we have:

• •

step 1 (BH spectrum):

• –

$\texttt{init\_masses[]}=8\times\texttt{BH\_number}$,

• –

$\texttt{init\_params[]}=8\times\texttt{BH\_number}$,

• –

$\texttt{spec\_table[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}$,

• •

step 2 (BH evolution):

• –

$\texttt{init\_masses[]}=8\times\texttt{BH\_number}$,

• –

$\texttt{init\_params[]}=8\times\texttt{BH\_number}$,

• –

$\texttt{spec\_table[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}$,

• –

$\texttt{fM\_table[][]}=8\times\texttt{nb\_fM\_param}\times\texttt{nb\_fM\_masses}$,

• –

$\texttt{gM\_table[][]}=8\times\texttt{nb\_fM\_param}\times\texttt{nb\_fM\_masses}$ if metric = 0,

• –

$\texttt{fM\_masses[]}=8\times\texttt{nb\_fM\_masses}$,

• –

$\texttt{fM\_param[]}=8\times\texttt{nb\_fM\_param}$,

• –

$\texttt{life\_masses[][][]}=8\times\texttt{BH\_number}^{2}\times\texttt{param\_number}^{2}\times\texttt{limit}$,

• –

$\texttt{life\_params[][][]}=8\times\texttt{BH\_number}^{2}\times\texttt{param\_number}^{2}\times\texttt{limit}$,

• –

$\texttt{life\_times[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}\times\texttt{limit}$,

• –

$\texttt{dts[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}\times\texttt{limit}$,

• •

step 3 (primary and secondary spectra):

• –

$\texttt{spec\_table[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}$,

• –

$\texttt{life\_masses[][][]}=8\times\texttt{BH\_number}^{2}\times\texttt{param\_number}^{2}\times\texttt{limit}$,

• –

$\texttt{life\_params[][][]}=8\times\texttt{BH\_number}^{2}\times\texttt{param\_number}^{2}\times\texttt{limit}$,

• –

$\texttt{life\_times[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}\times\texttt{limit}$,

• –

$\texttt{dts[]}=8\times\texttt{BH\_number}\times\texttt{param\_number}\times\texttt{limit}$,

• –

$\texttt{gammas[][][]}=8\times 4\times\texttt{nb\_gamma\_param}\times\texttt{nb\_gamma\_x}$,

• –

$\texttt{gamma\_param[]}=8\times\texttt{nb\_gamma\_param}$,

• –

$\texttt{gamma\_x[]}=8\times\texttt{nb\_gamma\_x}$,

• –

$\texttt{fits[][][]}=8\times 4\times\texttt{nb\_gamma\_param}\times\texttt{nb\_gamma\_fits}$,

• –

$\texttt{dof[]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• –

$\texttt{spins[]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• –

$\texttt{masses\_primary[]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• –

$\texttt{times[]}\approx 8\times\langle N\rangle\times\texttt{BH\_number}\times\texttt{param\_number}$,

• –

$\texttt{energies[]}=8\times\texttt{E\_number}$,

• –

$\texttt{tables[][][][]}=8\times\texttt{nb\_fin\_part}\times\texttt{nb\_init\_en}\times\texttt{nb\_fin\_en}\times\texttt{nb\_fin\_part}$,

• –

$\texttt{initial\_energies[]}=8\times\texttt{nb\_init\_en}$,

• –

$\texttt{final\_energies[]}=8\times\texttt{nb\_fin\_en}$,

• –

$\texttt{partial\_hadronized\_spectra[][][]}=8\times\texttt{nb\_fin\_part}\times\texttt{E\_number}\times\texttt{nb\_fin\_en}$,

• –

$\texttt{partial\_primary\_spectra[][]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})\times\texttt{E\_number}$,

• –

$\texttt{partial\_integrated\_hadronized\_spectra[][]}=8\times\texttt{nb\_fin\_part}\times\texttt{nb\_fin\_en}$,

• –

$\texttt{masses\_secondary[]}=8\times\texttt{nb\_fin\_part}$.

Using the parameters of Appendix D.1, the arrays occupy at most $\mathcal{M}\sim 120\,$MB, mostly because of the hadronization tables.

For BlackHawk_inst we have:

• •

step 1 (BH spectrum):

• –

$\texttt{BH\_masses[]}=8\times\texttt{BH\_number}$,

• –

$\texttt{BH\_params[]}=8\times\texttt{param\_number}$,

• –

$\texttt{spec\_table[][]}=8\times\texttt{BH\_number}\times\texttt{param\_number}$,

• •

step 2 (primary spectra):

• –

$\texttt{BH\_masses[]}=8\times\texttt{BH\_number}$,

• –

$\texttt{BH\_params[]}=8\times\texttt{param\_number}$,

• –

$\texttt{spec\_table[][]}=8\times\texttt{BH\_number}\times\texttt{param\_number}$,

• –

$\texttt{gammas[][][]}=8\times 4\times\texttt{nb\_gamma\_param}\times\texttt{nb\_gamma\_x}$,

• –

$\texttt{gamma\_param[]}=8\times\texttt{nb\_gamma\_param}$,

• –

$\texttt{gamma\_x[]}=8\times\texttt{nb\_gamma\_x}$,

• –

$\texttt{fits[][][]}=8\times 4\times\texttt{nb\_gamma\_param}\times\texttt{nb\_gamma\_fits}$,

• –

$\texttt{dof[]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• –

$\texttt{spins[]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• –

$\texttt{masses\_primary[]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• –

$\texttt{instantaneous\_primary\_spectra[][]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})\times\texttt{E\_number}$,

• –

$\texttt{energies[]}=8\times\texttt{E\_number}$,

• •

step 3 (during hadronization):

• –

$\texttt{instantaneous\_primary\_spectra[][]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})\times\texttt{E\_number}$,

• –

$\texttt{energies[]}=8\times\texttt{E\_number}$,

• –

$\texttt{tables[][][][]}=8\times\texttt{nb\_fin\_part}\times\texttt{nb\_init\_en}\times\texttt{nb\_fin\_en}\times\texttt{nb\_fin\_part}$,

• –

$\texttt{initial\_energies[]}=8\times\texttt{nb\_init\_en}$,

• –

$\texttt{final\_energies[]}=8\times\texttt{nb\_fin\_en}$,

• –

$\texttt{masses\_secondary[]}=8\times\texttt{nb\_fin\_part}$,

• –

$\texttt{instantaneous\_hadronized\_spectra[][][]}=8\times\texttt{nb\_fin\_part}\times\texttt{E\_number}\times\texttt{nb\_fin\_en}$,

• •

step 3 bis (during integration):

• –

$\texttt{instantaneous\_primary\_spectra[][]}=8\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})\times\texttt{E\_number}$,

• –

$\texttt{energies[]}=8\times\texttt{E\_number}$,

• –

$\texttt{initial\_energies[]}=8\times\texttt{nb\_init\_en}$,

• –

$\texttt{final\_energies[]}=8\times\texttt{nb\_fin\_en}$,

• –

$\texttt{instantaneous\_hadronized\_spectra[][][]}=8\times\texttt{nb\_fin\_times}\times\texttt{E\_number}\times\texttt{nb\_fin\_en}$,

• –

$\texttt{instantaneous\_integrated\_hadronized\_spectra[][]}=8\times\texttt{nb\_fin\_part}\times\texttt{nb\_fin\_en}$.

Using the parameters of Appendix D.1, the arrays occupy at most $\mathcal{M}\sim 120\,$MB, mostly because of the hadronization tables.

VIII.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 obtain 12 characters per written number, which is 12 bytes. For BlackHawk_tot we have:

• •

file BH_spectrum.txt: $\mathcal{M}=12\times(\texttt{BH\_number}+1)\times(\texttt{param\_number}+1)$,

• •

file life_evolutions.txt: $\mathcal{M}\approx 12\times\texttt{BH\_number}\times(1.+2\times\texttt{param\_number})\times\texttt{BH\_number}\times\texttt{param\_number}\times\langle N\rangle$,

• •

files *_primary_spectrum.txt: $\mathcal{M}\approx 12\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})\times\texttt{E\_number}\times\langle N\rangle\times\texttt{BH\_number}\times\texttt{param\_number}$,

• •

files *_secondary_spectrum.txt: $\mathcal{M}\approx 12\times\texttt{nb\_fin\_part}\times\texttt{nb\_fin\_en}\times\langle N\rangle\times\texttt{BH\_number}\times\texttt{param\_number}$.

Using the parameters of Appendix D.1, the total written disk space is $\mathcal{M}\sim 70\,$MB.

For BlackHawk_inst we have:

• •

file BH_spectrum.txt: $\mathcal{M}=12\times 3\times\texttt{BH\_number}$

• •

file instantaneous_primary_spectra.txt: $\mathcal{M}=12\times\texttt{E\_number}\times(\texttt{particle\_number}+\texttt{grav}+\texttt{add\_DM})$,

• •

file instantaneous_secondary_spectra.txt: $\mathcal{M}=12\times\texttt{nb\_fin\_en}\times\\ \texttt{nb\_fin\_part}$.

Using the parameters of Appendix D.1, the total written disk space is $\mathcal{M}\sim 70\,$kB.

IX Other applications

In this Section we present some hints about how to modify BlackHawk. Most of these modifications will require adding extra parameters in the parameter files and thus a modification of the routines read_params and read_*_infos, and of the structure struct param.

IX.1 Computing new numerical tables

The user may be interested in recomputing the tables described in Appendix C, 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 subfolder:

src/tables/

• •

create/modify the infos.txt information file and the corresponding read_*_infos routine,

• •

add a case 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.

IX.2 Using another BH mass and spin distribution

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 routines M_dist and param_dist,

• •

modify the parameter tmin, or the built-in Eq. (1) 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 $-1$, putting the table file in the subfolder:

src/tables/users_spectra/

and giving its full file name (including the .txt extension for example) to the parameter table. The format has to be:

• •

the first column for BH masses $M$, the first line for BH secondary parameter $x$ and the bulk of the table for the comoving number densities ${\rm d}n$ (with ${\rm d}M$ taken around $M$ and ${\rm d}x$ taken around $x$),

• •

masses and densities in CGS units (g and cm*-3* respectively), secondary parameters in dimensionless form,

• •

numbers in standard scientific notation,

• •

only additional text: a string of characters on the up-left corner (e.g. “mass/spin” for Kerr BHs).

IX.3 Adding primary particles

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

• •

increase the parameter particle_number in the greybody factor table information file infos.txt or add the new particle(s) with a switch similar to the one for the graviton,

• •

recompute the $f(M,x_{j})$ (and $g(M,a^{*})$ for Kerr BHs) 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 and modify the read_gamma_tables and read_gamma_fits routines,

• •

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.

We want to point out that as the graviton emission, as well as DM emission, are included in BlackHawk, new layers of new particles may become difficult to handle.

IX.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 (for example the file names or columns in the writing routines),

• •

add the corresponding contribution(s) to the routine contribution_instantaneous.

IX.5 Other types of black holes

If the user wants to compute the Hawking emission of BHs different from the already included 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,

• •

compute the corresponding $f$, $g$ and eventually new evolution parameters tables and add the corresponding reading routines and associated information in the infos.txt file,

• •

compute the new greybody factor tables and update the corresponding reading and interpolating routines read_gamma_tables, read_gamma_fits and dNdtdE.

Depending on the complexity of the BH model, the user may need to implement some or all of the above modifications. We want to point out that several BH metrics are already present inside BlackHawk v2.0: Reissner–Nordström BHs, higher-dimensional BHs and polymerized BHs, on top of the usual Kerr BHs.

X Conclusion

BlackHawk is the first public code generating both primary and secondary Hawking evaporation spectra for any distribution of black holes, and their evolution in time.555While the latest version of the manual was being written, Refs. [56, 57] appeared, describing the code ULYSSES aimed at leptogenesis, in which Hawking radiation of primordial black holes is computed in a similar way as BlackHawk. Its authors have checked for consistency with the results from BlackHawk. The primary spectra are obtained using greybody factors, and the secondary ones result from the decay and hadronization of the primary particles. The black holes and spectra evolution are obtained by considering the energy (and angular momentum for Kerr black holes) losses 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 primary 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. BlackHawk v2.0 embeds new features compared to the previous version, most of which deal with beyond Standard Model physics: emission of dark matter by Hawking radiation, spin $3/2$ particles or black hole solutions beyond standard general relativity and usual formation models (Reissner–Nordström, higher-dimensional and polymerized black holes). The code has thus been extended to allow for beyond the Standard Model studies (for illustrating examples and a full review of the BlackHawk related literature, see the release note [1]).

Acknowledgments

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. We also thank L. Morrison, S. Profumo and A. M. Coogan for help using the Hazma package. The authors thank the CERN Theory Department for its hospitality during which part of this work was done.

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_{\rm B}=G=4\pi\varepsilon_{0}=1$ to CGS or SI.

A.1 Energy

The energy conversion from GeV to Joules is

[TABLE]

A.2 Mass

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

[TABLE]

A.3 Time

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

[TABLE]

A.4 Distance

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

[TABLE]

A.5 Temperature

The dimensional link between energy and temperature is $[T]=[E/k_{\rm B}]$, and the conversion from GeV to Kelvins is

[TABLE]

A.6 Electric charge

The dimensional link between energy and electric charge is $[Q^{2}]=[4\pi\varepsilon_{0}\hbar c]$, hence electric charge is dimensionless in the natural system of units and we have, through the definition of the fine structure constant $\alpha$,

[TABLE]

where $q^{*}$ is the charge number.

Appendix B Particle information

Appendix C Numerical tables

In this Appendix, we detail the way the numerical tables are computed inside BlackHawk. While increasing the transparency of the code, this also gives hints concerning how to modify the tables or compute them with different methods.

C.1 Greybody factors

C.1.1 General framework

The greybody factors describe the probability that a particle, generated by quantum fluctuations at the horizon of a BH, escapes to spatial infinity. They are the only factor differing from the blackbody radiation in Eq. (34). They are obtained by solving the equations of motion of (beyond the) SM particles in curved spacetime, with special boundary conditions (for a detailed discussion on quantum fields, see e.g. [58, 59, 60]). For massive bosons and spin $1/2$ fermions, the free equations of motion are respectively the Proca and Dirac equations

[TABLE]

while for massive spin $3/2$ fermions it is the Rarita–Schwinger equation

[TABLE]

where $\mu$ is the particle rest mass and $\gamma_{\mu}$ are the Dirac matrices, with $\sigma^{\mu\nu}\equiv i[\gamma^{\mu},\gamma^{\nu}]/2$. Then, we proceed to a separation of these equations into a radial and an angular part; the separation is always possible for spherically symmetric and static BHs, as well as for Kerr BHs. To do that, we decompose the wave function $\Psi$ using the symmetries of the metric (Killing vectors $\partial_{t}$ and $\partial_{\varphi}$) in the Boyer–Lindquist coordinates $t,r,\theta,\varphi$

[TABLE]

where $E$ is the particle energy and $S_{lm}^{s}$ is the spin-weighted spheroidal harmonics for angular momentum $l$, projection $m$ and field spin $s$. The angular part of the wave equation gives some separation constant $\lambda_{s,lm}$ whose analytical expansion in terms of $\gamma\equiv EMa^{*}$ for Kerr BHs is given e.g. in [61] (valid for all spherically symmetric and static BHs as well with $a^{*}=0$), and the radial part satisfies what is denoted as the Teukolsky equation, whose general mathematical form is very complicated (for a complete mathematical overview of this topic, see e.g. [62]). This equation can be numerically solved directly, but the convergence of the solutions at the integration boundaries — namely at BH horizon and at spatial infinity — is not evident. Several authors have thus further transformed this equation into a Schrödinger-like wave equation with short-ranged potentials, different for each field spin, allowing for a more careful numerical treatment.

C.1.2 Kerr black holes

Chandrasekhar and Detweiler have shown that the Teukolsky equation can be reduced to a Schrödinger-like wave equation for Kerr BHs [63, 64, 65, 66]. They found necessary to define a Eddington–Finkelstein radial coordinate $r^{*}$

[TABLE]

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

[TABLE]

where $r_{\rm S}=2M$ is the BH Schwarzschild radius. Unfortunately, the inverse of this equation has to be found numerically and is generally difficult to determine with accurate precision. The Schrödinger-like wave equation is for all spins $s$

[TABLE]

The method to transform the Teukolsky equation into this simple wave equation was proposed in the Chandrasekhar and Detweiler papers [63, 64, 65, 66]. It is indeed difficult to find short-range potentials allowing for precise numerical computations. They give the form of such potentials in [63, 64] for spin 2, [65] for spins 0 and 1 and [66] for spin $1/2$. The potentials are666We found that the spin 0 potential had a missing “$r$” in [65].

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\nu\equiv\lambda_{2\,lm}+4$. Ref. [67] provides ingredients for the spin $3/2$ which we have for now included only in the Schwarzschild case. The potential is in this case

[TABLE]

where the intermediate quantities are this time

[TABLE]

[TABLE]

[TABLE]

In the Schwarzschild limit ($a^{*}=0$), we recover the Regge–Wheeler potentials [68]. As the angular momentum projection $m$ only appears multiplied by $a$, in this case the calculation is simplified since only one common value for all $m$ has to be chosen once $l$ is fixed. The sum over $m$ thus reduces to a factor $2l+1$ and the $r(r^{*})$ relation of Eq. (54) is analytical.

The $r^{*}$ variable change used in these potentials leads to divergences in the potentials, when $r^{2}=r_{\rm div}^{2}\equiv-\alpha^{2}$. This can happen for sufficiently low energies and high (negative) angular momentum projections, and it corresponds to the superradiance regime (for a very good review on the topic, see [69]). As discussed in the Chandrasekhar and Detweiler papers, the technique to avoid this divergence is to integrate Eq. (56) up to slightly before the divergence (i.e. $r_{\rm div}-\epsilon$). At this point, the behaviour of the potential $V_{s}$ is known, Eq. (56) is simplified, and the asymptotic form of the function $\psi_{s}$ can be obtained for $\epsilon\rightarrow 0$. By continuity of the wave function $\psi_{s}$ one can extrapolate this form up to slightly after the divergence (i.e. $r_{\rm 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 or spin $3/2$ potentials because of the $q-\beta_{\pm}\Delta$ and $q-\beta_{\pm}\sqrt{\Delta}$ terms. For this extra divergence, we try to integrate with one of the potentials (e.g. $\kappa_{+}$, $\beta_{+}$), and in case of a problem we try with the other potentials (e.g. $\kappa_{+}$, $\beta_{-}$), as it seems that at least one of the four combinations does not generate any divergence.

C.1.3 Spherically symmetric and static black holes

For spherically symmetric and static BH solutions of general form

[TABLE]

the authors, among others, have recently shown how any of those BH metric leads to a separation between the angular part and the radial part, and finally to a Schrödinger-like wave equation of a form similar to Eq. (56) with the following potentials [47, 48]

[TABLE]

where the spin-dependent parameters $\nu_{i}$ are given by $\nu_{0}=l(l+1)=\nu_{1}$, $\nu_{2}=l(l+1)-2$ and $\nu_{1/2}=l(l+1)+1/4$ and the Eddington–Finkelstein coordinate is defined by

[TABLE]

with the notation $\partial_{*}\equiv\partial_{r^{*}}$. The $r(r^{*})$ relation can be analytical or not, depending on the details of the metric (see Appendix A of [48]), but usually does not show a pathological behaviour. These potentials reduce to the Regge–Wheeler ones in the Schwarzschild limit $F(r)=G(r)=1-r_{\rm S}/r$ and $H(r)=r^{2}$ [68].

C.1.4 Boundary conditions

As boundary conditions to solve Eq. (56), we use a purely ingoing wave. The solution at the horizon has the form

[TABLE]

At infinity, the solution has the form

[TABLE]

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

[TABLE]

Practically, we compute the value of a single helicity/color dof emissivity

[TABLE]

for:

• •

50 values of $0\leq a^{*}\leq 0.9999$ for Kerr BHs,

• •

50 values of $0\leq Q^{*}\leq 0.999$ for Reissner–Nordström BHs,

• •

7 values of $0\leq n\leq 6$ for higher-dimensional BHs,

• •

22 values of $0\leq\varepsilon\leq 100$ for polymerized BHs (with $a_{0}=0$ and $a_{0}=0.11$),

and for a range of 200 energies $0.01\leq x\equiv Er_{\rm S}\leq 5$ (dimensionless). For $x$ out of this range, we have found easier to fit empiric asymptotic forms to the emissivities. At low energies, we have for all spins $s$

[TABLE]

and at high energies for the Kerr BH

[TABLE]

which allows for a reconstruction of the oscillatory behaviour in this limit [61]. For spherically symmetric and static BHs, we chose at high energy to keep only

[TABLE]

as the oscillatory behaviour is damped and the constant cross-section regime is attained rapidly [48]. We checked that the fitting parameters agree with the asymptotic limits of [28] in the Schwarzschild case and of [50] in the Kerr case, and with analytical formulas we derived in [48] in the spherically-symmetric and static case. For Kerr BHs, the Mathematica scripts as well as a C formatting script and a README are provided in the folder:

scripts/greybody_scripts/greybody_factors/spin_*.m

scripts/greybody_scripts/greybody_factors/exploitation.m

scripts/greybody_scripts/greybody_factors/formatting.c

scripts/greybody_scripts/greybody_factors/Makefile

scripts/greybody_scripts/greybody_factors/README.txt

For the other BH solutions, the Mathematica and Python scripts are in the relevant folders:

scripts/greybody_scripts/greybody_factors//spin_.m

scripts/greybody_scripts/greybody_factors//exploitation_.py

If these tables are recomputed, it is advisable to also modify the information in the file:

src/tables/gamma_tables/infos.txt

C.2 Page factors

To compute the integrals of Eqs. (38) and (39), we use the greybody factor tables and the fits computed in the previous Section. The peak of Hawking emission lies around the BH temperature for spherically symmetric and static BHs (see e.g. [48]), and is uncorrelated for Kerr BHs (see e.g. [61]); in any case it is close to the Schwarzschild temperature for reasonable values of the secondary parameters. Thus the integral does not need to be computed over all energies, but a restrained set of energies in $10^{-5}\times T<E<10^{5}\times T$ is sufficient. The domains of integration are segmented over 1000 logarithmically distributed energies, and computed for 1000 masses $M_{\rm P}<M<10^{46}\,$GeV ($\sim 10^{-5}<M<10^{22}\,$g). We compute these tables with and without the graviton emission, and with or without a single dof of DM emission. In the limit where neutrinos are massless, $f(M,x_{j})$ 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,x_{j})$ and $g(M,a^{*})$ are in GeV4 (corresponding to g${}^{3}\cdot$s*-1* and g${}^{2}\cdot$GeV$\cdot$s*-1*, respectively). We have checked that the value of $f(M,0)$ is consistent with that of [29] in the Schwarzschild case and that the values of $f(M,a^{*})$ and $g(M,a^{*})$ are consistent with [70] in the Kerr case. No numerical results were ever given for the polymerized case. The C scripts used to compute the tables, a README and all the necessary greybody factors tables are provided in the subfolder:

scripts/greybody_scripts/fM/fM.c

scripts/greybody_scripts/fM/Makefile

scripts/greybody_scripts/fM/README.txt

scripts/greybody_scripts/fM/spin_*.txt

If these tables are recomputed, it is advisable to also modify the information in the file:

src/tables/fM_tables/infos.txt

C.3 Hadronization tables

Three particle physics codes have been used to compute hadronization tables: PYTHIA [25], HERWIG [26] and Hazma [27]. In the first two cases, the strategy is to generate the output of a collision (for example $e^{+}+e^{-}\rightarrow u+\overline{u}\rightarrow\dots$), and then to count the number of final particles (here denoted as dots) normalized by the number (here 2) of initial particles (here $u$) satisfying the desired stability criterion: Table 4 for early Universe/BBN particles (pythia_tables and herwig_tables) and Table 5 for present epoch particles (pythia_tables_new). This strategy is adapted for those two codes which are Monte-Carlo simulators. Hazma is quite different as it computes analytically the branching ratios once the initial energy is fixed.

To build the PYTHIA and HERWIG tables, we have simulated for each channel listed in Table 6, $10^{5}$ events for initial energies $E^{\prime}$ (half of the center of mass energy) logarithmically distributed between $5\,$GeV and $10^{5}\,$GeV (PYTHIA) 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 quantities ${\rm d}N_{j}^{i}(E^{\prime},E)$ of Eq. (41). For the Hazma tables, we have directly applied the analytical computation of the branching ratios to arrays of initial energies $10^{-6}-5\,$GeV for all the available particles: electrons, muons and pions. The procedure is well described in [27, 71].

The branching ratios $e^{\pm}\rightarrow\gamma\gamma\rightarrow...$ and $e^{\pm}\rightarrow\nu\overline{\nu}\rightarrow...$ have not been computed with PYTHIA and HERWIG. The contribution from the primary photons, neutrinos and electrons are directly added to the secondary spectra with a branching ratio of 1, which is an approximation. Hazma, on the other hand, handles the final state radiation and the decay of pions and muons into photons and electrons, thus in this case only the neutrino final spectrum is out of reach.

For initial energies lower than the cutoff of the PYTHIA and HERWIG 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. This is a mere numerical trick that was shown by [1, 71] to be irrelevant for BH masses close to the QCD scale $M\sim 10^{15}\,$g. To be clear: PYHTIA and HERWIG are particle physics codes that are very efficient in their domain of validity, for energies between some GeV and some TeV (BH mass $M<10^{13}\,$g), and considers SM particles (e.g. quarks) as elementary degrees of freedom. However, Hazma is much more fitted to compute the low energy showering of (charged) particles and considers pions as fundamental degrees of freedom below the QCD scale of $\sim 200\,$MeV (BH mass $M\gtrsim 10^{14}\,$g). As their is no sizeable secondary generated particles for energies below the MeV scale (BH mass $M\gtrsim 10^{17}\,$g), all these codes should provide the same result.

The scripts used to compute the pythia_tables are given in the folder:

scripts/pythia_scripts/formatting.c

scripts/pythia_scripts/Makefile

scripts/pythia_scripts/README.txt

scripts/table_*.cc

The scripts used to compute the pythia_tables_new are given in the folder:

scripts/pythia_new_scripts/formatting.c

scripts/pythia_new_scripts/Makefile

scripts/pythia_new_scripts/README.txt

scripts/pythia_new_scripts/table_*.cc

The scripts used to compute the herwig_tables are given in the folder:

scripts/herwig_scripts/formatting.c

scripts/herwig_scripts/Makefile

scripts/herwig_scripts/README.txt

scripts/herwig_scripts/BH_TABLE_GENERATOR.cc

scripts/herwig_scripts/MC_PARTICLE_COUNTS.cc

scripts/herwig_scripts/RivetMC_PARTICLE_COUNTS.so

scripts/herwig_scripts/*/LEP.in

scripts/herwig_scripts/*/main.cpp

The script used to compute the Hazma tables is given in the folder:

scripts/hazma_scripts/hazma_tables.py

Please contact one of the authors if you have issues using these scripts. If these tables are recomputed, it is advisable to also modify the informations in the files:

src/tables/hadronization_tables/infos.txt

Appendix D Example runs

In this Section, we provide the expected results of a run of each of the BlackHawk programs, when using the built-in parameters.txt file. They should be strictly equivalent to the user’s results.

D.1 Built-in parameters.txt file

The built-in parameters are:

destination_folder = example

full_output = 1

interpolation_method = 0

BH_number = 1

Mmin = 1e+9

Mmax = 1e+16

metric = 0

param_number = 1

amin = 0.

amax = 0.5

Qmin = 0.

Qmax = 0.7

epsilon_LQG = 0.

a0_LQG = 0.

n = 0.

spectrum_choice = 0

spectrum_choice_param = 0

amplitude_lognormal = 1.

amplitude_lognormal2 = 1.

stand_dev_lognormal = 1.

crit_mass_lognormal = 1.

amplitude_powerlaw = 1.

eqstate_powerlaw = 0.3333

amplitude_critical_collapse = 1.

crit_mass_critical_collapse = 1.

amplitude_uniform = 1.

stand_dev_param_gaussian = 1.

mean_param_gaussian = 0.5

table = spin_distribution_BH.txt

tmin_manual = 0

tmin = 1.e-30

limit = 5000

BH_remnant = 0

M_remnant = 1e-4

E_number = 10

Emin = 5

Emax = 1e+5

grav = 0

add_DM = 0

m_DM = 0.

spin_DM = 0.

dof_DM = 0.

primary_only = 0

hadronization_choice = 0

D.2 BlackHawk_tot output files

When launching BlackHawk_tot with the built-in parameters.txt file, one obtains the following output files (truncated to 5 lines and 5 rows when relevant):

• •

BH_spectrum.txt

Initial BH comoving number density as a function of their mass and parameter.

mass/spin 0.00000e+00

1.00000e+09 1.00000e+00

• •

life_evolutions.txt

Evolution of the BH masses and spins as functions of time.

Total number of time iterations: 1044

t M a

1.23833e-29 1.00000e+09 0.00000e+00

2.47665e-29 1.00000e+09 0.00000e+00

3.71498e-29 1.00000e+09 0.00000e+00

4.95331e-29 1.00000e+09 0.00000e+00

6.19164e-29 1.00000e+09 0.00000e+00

• •

dts.txt

Evolution of the integration timestep as a function of time.

t dt

1.23833e-29 1.23833e-29

2.47665e-29 1.23833e-29

3.71498e-29 1.23833e-29

4.95331e-29 1.23833e-29

6.19164e-29 1.23833e-29

• •

photon_primary_spectrum.txt

Hawking primary spectrum as a function of time.

time/energy 5.00000e+00 1.50267e+01 4.51601e+01 1.35721e+02 4.07886e+02

1.23833e-29 4.10141e+09 1.02586e+11 2.56592e+12 6.41798e+13 1.60529e+15

2.47665e-29 4.10141e+09 1.02586e+11 2.56592e+12 6.41798e+13 1.60529e+15

3.71498e-29 4.10141e+09 1.02586e+11 2.56592e+12 6.41798e+13 1.60529e+15

4.95331e-29 4.10141e+09 1.02586e+11 2.56592e+12 6.41798e+13 1.60529e+15

6.19164e-29 4.10141e+09 1.02586e+11 2.56592e+12 6.41798e+13 1.60529e+15

• •

photon_secondary_spectrum.txt

Hawking secondary spectrum as a function of time.

time/energy 1.00000e-06 1.05207e-06 1.10685e-06 1.16448e-06 1.22511e-06

1.23833e-29 9.55768e+30 1.14321e+31 1.87205e+31 1.90536e+31 1.92737e+31

2.47665e-29 9.55768e+30 1.14321e+31 1.87205e+31 1.90536e+31 1.92737e+31

3.71498e-29 9.55768e+30 1.14321e+31 1.87205e+31 1.90536e+31 1.92737e+31

4.95331e-29 9.55768e+30 1.14321e+31 1.87205e+31 1.90536e+31 1.92737e+31

6.19164e-29 9.55768e+30 1.14321e+31 1.87205e+31 1.90536e+31 1.92737e+31

D.3 BlackHawk_inst output files

When launching BlackHawk_inst with the built-in parameters.txt file, one obtains the following output files (truncated to 5 lines and 5 rows when relevant):

• •

BH_spectrum.txt is the same as for BlackHawk_tot

• •

instantaneous_primary_spectra.txt

Hawking primary spectra for each particle types.

energy/particle photon gluons higgs W+- Z0

5.00000e+00 4.10141e+09 3.28113e+10 0.00000e+00 0.00000e+00 0.00000e+00

1.50267e+01 1.02586e+11 8.20688e+11 0.00000e+00 0.00000e+00 0.00000e+00

4.51601e+01 2.56592e+12 2.05274e+13 0.00000e+00 0.00000e+00 0.00000e+00

1.35721e+02 6.41798e+13 5.13438e+14 8.96189e+19 1.92539e+14 9.62696e+13

4.07886e+02 1.60529e+15 1.28423e+16 2.50679e+20 4.81586e+15 2.40793e+15

• •

instantaneous_secondary_spectra.txt

Hawking secondary spectra for each particle types.

energy/particle photon electron muon nu_e nu_mu

1.00000e-06 9.55768e+30 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

1.05207e-06 1.14321e+31 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

1.10685e-06 1.87205e+31 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

1.16448e-06 1.90536e+31 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

1.22511e-06 1.92737e+31 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00