Physics Beyond the Standard Model with BlackHawk v2.0

We present the new version v2.0 of the public code BlackHawk designed to compute the Hawking radiation of black holes, with both primary and hadronized spectra. This new version aims at opening an avenue toward physics beyond the Standard Model (BSM) in Hawking radiation. Several major additions have been made since version v1.0: dark matter/dark radiation emission, spin $3/2$ greybody factors, scripts for cosmological studies, BSM black hole metrics with their associated greybody factors and a careful treatment of the low energy showering of secondary particles; as well as bug corrections. We present, in each case, examples of the new capabilities of BlackHawk.

Since Hawking first demonstrated that black holes (BHs) evaporate by emitting a radiation close to the thermal radiation of a black body [1,2], this phenomenon has been extensively studied. One of the main outcome of Hawking radiation (HR) is the possibility that small BHs (i.e. with mass M BH M ) formed just after the end of the universe expansion, namely primordial BHs (PBHs), may emit or have emitted radiation that could be observable today or that could have left an imprint in cosmology. This leads to a number of constraints on the PBH abundance, depending on their mass M BH , dimensionless spin a * , or the extended distribution of both these parameters. For recent reviews on PBH formation and constraints, see Refs. [3,4]. When considering PBHs with mass M BH M eva where M eva 5 × 10 14 g is the mass of (Schwarzschild) PBHs evaporating just today if formed at the beginning of the universe, these constraints are given as the fraction of dark matter (DM) f PBH that the PBHs represent today. Indeed, light PBHs can represent some fraction of DM as they contribute to cold non baryonic matter in the early universe thermal history [5]. For lighter BHs with mass M BH M eva , these constraints are given as the maximum possible fraction of the universe β collapsed into PBHs at their formation. These BHs cannot represent a significant fraction of DM today since they have evaporated away. However, if evaporation stops at some point, leaving Planck scale remnants, these may contribute to DM [6][7][8][9][10][11]. To determine the PBH evaporation constraints, one compares the Hawking radiation signals to astrophysical or cosmological observations such as e.g. the extragalactic gamma ray background (EGRB) or the ∆N eff limits respectively. PBHs in the mass range just above the evaporation limit, where energy is in the low energy QCD domain E ∼ MeV, are severely constrained as a fraction of DM. This represented a difficulty for BlackHawk v1.2, which relied on PYTHIA and HERWIG to compute the hadronization of primary particles, valid for E GeV. There was a need for more careful treatment of this energy range, a work initiated by the study of Ref. [12]. Noteworthy, HR is a purely gravitational phenomenon. As such, BHs would radiate any dark sector particle on top of the Standard Model (SM) spectrum, such as supersymmetric states [13], DM or dark radiation (DR) [14][15][16][17][18][19][20][21][22][23][24] -and the graviton [8,19,[25][26][27][28][29] (which contributes to DR).
It is thus of utmost importance to be able to predict precisely the HR spectra of SM and beyond SM (BSM) particles. A lot of analytical and numerical results exist in the literature for this purpose, the ones of [30][31][32][33] being among the most cited. However, people use most of the time the high-energy limit for the greybody factors given in these papers, resulting in order-of-magnitude constraints. There was thus a need for an automatic program to compute the precise spectra for any BH mass and spin, with a wide choice of parameters that allow for any kind of HR study. This was the context of creation of the public code BlackHawk v1.0 [34]. Since its first release, the code has received some modifications and has reached version 1.2. The previous version of the manual can be found on the arXiv [34] (v2), while the code is publicly available on HEPForge: https://blackhawk.hepforge.org/ The code has been recently presented to the TOOLS 2020 conference [35]. BlackHawk is used by many groups from very different domains of astrophysics and cosmology to perform striking studies including, to the knowledge of the authors, evolution of BHs spin [36], EGRB constraints with extended mass distributions and spinning BHs [37] or with higher dimensional Schwarzschild BHs [38], electron and positrons signals from the galaxy with the 511 keV line [39,40], current [41] or prospective [42] X-ray limits, neutrino constraints from Super-Kamiokande [40], JUNO [43] or prospective neutrino detectors [44,45], gamma ray constraints from INTEGRAL [46], COMPTEL with improved lowenergy secondary particles treatment [12], prospective AMEGO instrument [12,47], LHASSO [48] or fine modelisation of the Galaxy [49], prediction of signals from Planet 9 within the PBH hypothesis [50], archival galactic center radio observations [51], interstellar medium temperature in dwarf galaxies [52][53][54] or 21 cm measurements by EDGES with Schwarzschild [55] or Kerr PBHs [56], Big Bang nucleosynthesis (BBN) [57], heat flow from a small BH captured in the Earth core [58], warm DM from light Schwarzschild [22] and Kerr [23] PBHs, dark radiation from light spinning PBHs [23,24], (extended) dark sector emission [13,59], axion-like particle emission [60], HR from extended BH metrics [61,62]. While this release note was being written, the authors became aware of Refs. [63,64] which analyze the production of DM by PBH evaporation in the early universe; while these do not explicitly use the code BlackHawk, apart from the spin 2 greybody factors, they provide results compatible with those of BlackHawk with mostly the same spirit in the numerical part. This paper is a release note of a new version of BlackHawk v2.0 that includes some primordial new features linked to the physics described above: dark sector emission, spin 3/2 greybody factors, BH remnants, BSM BH metrics, low energy hadronization; as well as bug corrections. It then completes the manual published in European Physics Journal C [34] and goes together with the updated version of the code available on the BlackHawk website mentioned above, and an updated version of the manual available on the arXiv [34] (v3). All the installation and run procedures, as well as the complete set of parameters and routines are described in the latter, while we focus here on the new features only. It is organized as follows: in Section II we describe the physics of the new features added to the code and illustrate them with examples of BlackHawk results, in Section III we list the new input parameters, in Section IV we mention some improvements on the optimisation of the code and we conclude in Section V. In the following, we use the natural system of units G = = c = k B = 1.

II. NEW FEATURES
In this Section, we describe the new features added to BlackHawk since the manual was published [34] and we give examples of interesting new results. They have been grouped into distinct categories, even if the changes were made separately: addition of DM, spin 3/2 greybody factors, time dependent features for cosmological studies, greybody factors from new BH metrics, low energy hadronization. Most of these add-ons imply a modification of some parameters and routines, that are listed respectively in Section III and Appendix A.
A. Adding a DM particle A persistent question of the BlackHawk users was about the possibility of adding a (BSM) particle to the computed SM spectra. The procedure to implement this modification was described in the BlackHawk manual, but the authors found it very convenient to add this feature to the vanilla version of the code. This addition of a new particle was used in warm DM [22,23] and DR [24] studies where it allowed precision improvement compared to previous analytical approximations.
We need to understand what changes when a particle beyond the SM is Hawking radiated. In fact, as the process is purely gravitational and as the emission of individual particles is independent, any number of additional degrees of freedom can be added without changing the rest of the spectra. The only modification is at the level of the BH lifetime. Indeed, in the integration over the total spectrum to obtain the mass and spin change of the BH, additional dofs provide additional contributions: the mass and spin loss rates are enhanced, thus the BH evaporates faster. When considering only one additional dof, as in [22][23][24], the change in the lifetime is negligible, but it could be much higher with numerous additional dof as considered in [13]. In BlackHawk, this change has been taken into account in alternative f (M, a * ) and g(M, a * ) tables computed by an updated (and simplified) version of the script fM.c. Tables fM add*.txt, fM add* nograv.txt, gM add*.txt and gM add* nograv.txt are read instead of the usual fM.txt, fM nograv.txt, gM.txt and gM nograv.txt, for the addition of a single massless dof of spin *. The tables f (M, a * ) and g(M, a * ) have further been recomputed with increased precision. The greybody factor tables and fits necessary to recompute the f, g tables are also given in the code source. All these files are available in the folder: /scripts/greybody scripts/fM If the additional particle considered has a rest mass µ that has to be taken into account, new tables must be computed using the fM.c script to limit the emission of this particle to the usual condition E > µ. Finally, the primary spectrum output is modified by the addition of one column to the file instantaneous primary spectra.txt (program BlackHawk inst) or by a new file DM primary spectrum.txt (program BlackHawk tot). The additional dof is assumed to have no interaction with the SM. Hence, there is no modification of the SM secondary spectra. If it were the case, the user would have to modify the routines computing the secondary spectra to implement any new branching ratio (e.g. DM decay into SM particles). We want to stress that within BlackHawk, the graviton is not considered as a BSM particle and is already treated separately since BlackHawk v1.0.
Tab. I below shows the lifetime reduction when one adds a single massless (DM) dof to the SM. As expected, the relative change is small, and depends on the new particle spin, with a decreased impact as the spin increases for a Schwarzschild BH (spin a * = 0). On the contrary, when the BH is close to Kerr extremal (spin a * = 0.9999), higher spin new particles have a greater impact compared to the Schwarzschild case. The spin hierarchy of the impact of adding one dof is not inverted, as one could expect knowing that higher spin particles have a rate of emission much higher than lower spin particles for a close to extremal Kerr BH. The reason is that as the BH evolves, it loses its spin quite fast and behaves like a Schwarzschild BH for most of its lifetime [36], where the low spin particles have the higher impact.

B. Spin 3/2 particle
When doing DM and DR studies, it became obvious that particles with spin 3/2 may play a role [22,23]. These particles could be motivated by e.g. supersymmetric extensions of the SM; one example is the gravitino, massive spin 3/2 superpartner of the graviton.
The Newman-Penrose form of the Rarita-Schwinger equation of motion for the massless spin 3/2 field in Kerr-Newman geometry, the corresponding radial Teukolsky equation, as well as the short range potential V 3/2 (r(r * )) obtained following the Chandrasekhar transform of the Teukolsky equation are described in [65][66][67]. We must then integrate the equation of motion of the spin 3/2 field from the BH horizon to spatial infinity The greybody factor for some partial spin-weighted spherical harmonic s, l is given as the probability that the emitted particle reaches infinity, expressed in terms of the wave function amplitudes as where the superscript "out" refers the the outgoing part of the wave. For the moment, in BlackHawk, we have only implemented the Schwarzschild special case of the spin 3/2 Hawking radiation. The massive spin 3/2 particle can be treated as the other spins, by putting an artificial energy cutoff in the Hawking spectrum at the particle rest mass E > µ.
We have computed the greybody factors thanks to a script similar to the ones used for spins 0, 1, 2 and 1/2 particles, in the Schwarzschild case. The Mathematica script for spin 3/2 is spin 1.5.m; it has been added to the BlackHawk folder:

/scripts/greybody scripts/greybody factors
This script performs the integration of Eq. (1) and tabulates the greybody factors summed over the spherical harmonics decomposition of the wave function. For details about this computation see the complete manual [34]. We advise to use as a "spatial infinity" boundary condition at least r * ∞ /r H > 300/x, where x ≡ Er S is the adimensioned energy, for good numerical convergence. This is particularly important for the low energy greybody factors. The script fM.c has also been updated as already developed in Section II A to include spin 3/2 emission in the BH lifetime computation.
In Fig. 1 (left panel) we show the Hawking radiation rate for all massless spins fields with a single dof where s is the particle spin, T = 1/4πr S is the Schwarzschild BH temperature and where Γ l s is the greybody factor for angular momentum l, energy E and BH mass M BH . In the same Figure (right panel) we also show the high energy limit of the cross section for all spins compared to the geometrical optics (GO) approximation [30] σ GO ≡ 27 4 where r S ≡ 2M is the Schwarzschild horizon radius.

C. Time dependent features, cosmological studies
Several studies performed with BlackHawk consider Hawking radiation through cosmological eras, e.g. the early universe before BBN [22][23][24]. It can be relevant to consider not the full time-dependent spectrum of emitted particles, but only the stacked spectrum of those particles at some time after the BH evaporation. This quantity must be computed by taking the redshift into account. Indeed, particles emitted at the beginning of the BH evaporation see their energy diluted compared to those emitted at the end of the BH lifetime. We thus provide into BlackHawk a C script stack.c that performs this integration, with adjustable cosmology parameters such as the matter/radiation domination at the time of evaporation, the time of matter-radiation equality, etc. This script is available in the folder: /scripts/cosmology scripts and is accompanied by a README.txt file containing instructions. This kind of integration requires that we can extract the time steps of the BH evaporation precisely. It was not possible with the previous version of BlackHawk as only the absolute time was available in the output files, with a exponential precision of 5 digits. At the end of the BH evaporation, the time steps become so tiny that this precision is not sufficient to distinguish the time steps, making the integration procedure incomplete. Thus, we have added the output file dts.txt to the result set of a BlackHawk run which contains the explicit time steps on top of the absolute time. This file is read by the script stack.c.
The script computes the integrated spectrumQ i of particle i from BH formation t form to total evaporation t evã The redshift between running time t and the evaporation time t eva is taken into account to shift the energies as well as the BH distribution dilution. This results in the two factors 1 + z(t) where the redshift as a function of time is given by for radiation domination (RD) and matter domination (MD) respectively. The case of a mixed or alternate domination of radiation and matter could be straightforwardly implemented in the routine redshift() inside that script. On top of that, it has become important to treat precisely the boundaries of integral (7). This is why we have implemented in BlackHawk an explicit formula for the value of the BH formation time t form (see e.g. [3,22]) where γ 0.2 is the fraction of a Hubble sphere that collapses into a BH when the density perturbation re-enters the Hubble horizon. This formula is valid for BH formation during a radiation dominated era.
In Fig. 2 we show the (massless) graviton spectrum of a Schwarzschild BH stacked over its lifetime considering RD or MD. The shape can be compared to the instantaneous (massless) spin 2 spectrum of Fig. 1 (left panel).
On a related subject, multiple scenarios of quantum gravity predict that the evaporation of BHs stops at some mass greater than the Planck mass. What is left of the BHs, commonly denoted as BH relics or BH remnants, not to mix up with DM relics, does not evaporate and keeps a constant mass. If the BH was originally spinning or had an electric charge, it is possible that random fluctuations at the end of evaporation attribute a non-zero spin/charge to the remnant [68]. It is also possible that the usual HR paradigm breaks down at some mass scale above the Planck mass, where quantum effects are believed to be important, thus advocating for cautious treatment of the emission below this mass scale. These stable BH remnants, of a mass of the order of magnitude of or superior to the Planck mass, can constitute some fraction of DM [6][7][8][9][10][11] (see however the new discussion started by [69]). As the remnant mass depends on the quantum gravity model, and in particular on the number of extra spatial dimensions [70], we are agnostic about it and leave it as a free parameter of the code. Apart from the particular case of M remnant M BH , the BH lifetime and stacked spectra should not be perturbed much by the incomplete evaporation.

D. New BH metrics and greybody factors
Constraints over the PBH fraction of DM are more and more stringent. Would an experiment detect a signal compatible with PBH Hawking radiation, this would constitute a test of the BH quantum behaviour close to the horizon complementary to the study of BH quantum normal modes during the ringdown phase after a merger (see e.g. [71]). Then, one needs to know precisely what is the effect of exotic BH metrics on the HR signals. In that context, we have implemented greybody factors computed with different BH metrics in BlackHawk, namely higher dimensional BHs, charged (Reissner-Nordström) BHs and polymerized BHs. The first were already studied in [72], and more recently Ref. [38] added the higher dimensional greybody factors inside BlackHawk to re-evaluate the extragalactic gamma ray background constraints. Polymerized BHs are one example of a regular BH solution to the Einstein equations that exhibits a horizon with Hawking radiation. Hence, their HR signals, as well as those of other regular BH solutions, are of particular interest and are the subject of numerous recent studies (see e.g. [73][74][75][76][77]). Some analytical and numerical results about HR from polymerized BHs were presented in [78][79][80][81][82]. The literature on charged BHs is older as they are not astrophysically motivated [68]. We completed those results with a full analytical study of the HR by spherically symmetric BHs described by metrics of the type which are a subset of Petrov type D metrics. In these metrics, the 4-dimensional angular part is dΩ 2 = dθ 2 +sin(θ)dϕ 2 . We derived the short range potentials V s for the Schrödinger-like wave equation equivalent to the radial Teukolsky equation, for massless fields of spins 0, 1, 2 and 1/2 in Ref. [61]. We presented the Hawking radiation spectra in the companion paper [62]. In the literature, there existed a lot of results for tr-symmetric metrics of the common type F (r) = G(r) ≡ h(r) , and H(r) = r 2 .
Examples include higher dimensional BHs (emission on the brane), charged BHs, (A)dS BHs, etc. The general case where F (r) = G(r) received attention only recently when regular BH solutions of that kind were derived; we computed the greybody factors for polymerized BHs which are one physically motivated example within the loop quantum gravity (LQG) paradigm. The greybody factors for these metrics are computed in the same manner as the Kerr ones in the previous version of the code. The equations of motion of the massless fields of spins 0, 1, 2 and 1/2 (i.e. the Newman-Penrose equations) are separated into a radial and an angular part. The radial equation can then be transformed into a Schrödinger-like wave equation with a short range potential V s depending on the particle spin and the metrics details Solving this equation numerically with planar wave boundary conditions gives access to the reflection and transmission coefficients, the latter being nothing but the greybody factor. Analytical results were obtained in the usual low and high energy limits. For more details about these metrics and their Hawking radiation, see Refs. [61,62]. In Fig. 3, we show the full Hawking spectra for the massless (uncharged) spin 0 scalar in the case of a charged BH, a higher-dimensional BH, and a polymerized BH (left panel). In the same figure (right panel) we focus on the high energy oscillations of the cross-section.
Inside BlackHawk, these new metrics and their modified greybody factors and f Page factors are available through the parameters.txt file with the new parameter metric. The script fM.c has been modified to include the possibility of computing the function f (M, , a 0 ) for polymerized BHs, for 11 values of ∈ [0, 0.79] and a 0 = {0, 0.11}, and 11 values of ∈ [1, 100] and a 0 = 0, where the modified Page factor is given by where Γ i (E, M BH , , a 0 ) is the modified greybody factor for the polymerized metrics with parameters and a 0 . Contrary to the Kerr case, where the BH spin a * evolves during evaporation, here and a 0 are purely constant. We want to stress a particular point: if the polymerized BH mass gets too small, namely if r + (M, ) 2 ∼ a 0 , then the validity of the solution is yet to be investigated and we do not guarantee the reliability of the Hawking radiation results. The f tables have been computed with and without the graviton contribution and are stored in: It is of course possible to mix the new metrics with DM contribution to the Hawking spectrum (in that case, the Page factor f must be recomputed for polymerized BHs) following the details of Section II A, at the exception of the spin 3/2 additional particle, for which the greybody factors have not been computed within the new metrics. Using the precise Hawking radiation rates for primordial black holes, Ref. [62] derived the first ever constraints on polymerized primordial BHs as DM.

E. Low energy photon spectrum
It was obvious during the development of BlackHawk that secondary spectra would cause a difficulty due to the limits in the theoretical and numerical tools at hand to determine the evolution of SM particles after having been Hawking radiated. In the first version of the code, we used the public particle physics codes PYTHIA [83] and HERWIG [84] to convolve the primary spectra with hadronization and decay branching ratios. These codes are designed to match the data of particle accelerators and thus their domain of validity corresponds to the energy range of these accelerators, something like ∼ 5 GeV to ∼ 10 TeV. We decided to extend this range to 100 TeV, meaning that there is no unknown physics until there (the fundamental interactions are the same, there is no new (interacting) particle in this range). For higher and lower energies, we decided to extrapolate the branching ratios of PYTHIA and HERWIG, when the secondary particles were kinematically allowed. However this should break down at low energy. Indeed, as stated in [30,31], the secondary spectra hypothesis is that primary SM particles are emitted at all energies if kinematically allowed (E > µ where µ is the particle rest mass), then these fundamental particles are hadronized or decayed. For example, a down quark with a mass of 4.7 MeV can be emitted at an energy of E > 4.7 MeV. But as we know, quarks can not exist outside of bound states, and the lightest bound state is the neutral pion with µ π 0 135 MeV. The questions is now: are quarks of type q Hawking emitted between E = µ q and E = µ π 0 ? A question which has been reformulated elsewhere as: what are the fundamental particles at some low energy E? Indeed, we can argue that below some energy scale (typically the QCD scale), the fundamental degrees of freedom of the theory are the pions and not the individual quarks. Below the pion rest masses, there is no viable QCD degree of freedom, even if light quarks and gluons are kinematically allowed. That is the scheme chosen by [12] when deriving the secondary photon spectrum at low energy (below the QCD scale).
We thus modified BlackHawk such that it emits pions instead of single quarks if the user chooses so, with the new hadronization choice in the hadronization choice parameter. Then, we take back the development of [12] on the behaviour of the particles at low energy and use the public Python code Hazma [85] to evolve the primary particles.  4. Comparison between the old PYTHIA extrapolated (solid red) and new Hazma computed (solid black) BlackHawk instantaneous total photon spectrum for a MBH = 5.3×10 14 g. From left to right: primary photon spectrum (grey solid); neutral pion decay (orange dot-dashed); electron FSR (blue dashed); muon decay+FSR (purple dotted); charged pion decay+FSR (green dotted). To be compared with Fig. 2 upper left panel of [12].
We restrict ourselves to the secondary photon and electron spectra (Hazma does not provide neutrino spectra), which is an improvement over [12] which only considered secondary photons. Depending on the available energy, electrons, muons and pions can contribute to the secondary spectra, in addition to the primary photons and electrons. The contributions are: • final state radiation (FSR): pairs of created opposite charge particles radiate photons, this concerns electrons, muons and charged pions; • decays: single unstable particles decay into photons or electrons, this concerns muons, neutral and charged pions.
We follow [12] to compute the FSR photons thanks to the electroweak splitting functions [86] directly inside BlackHawk from the primary electron, muon and charged pions spectra. We also follow [12] and use the public code Hazma [85] to compute the low-energy pions and muon decay rates, which we have tabulated in: /src/tables/hazma tables/photon.txt which covers center of mass energy of 1 keV to 5 GeV. We improve on [12] by including the secondary low energy electron and positron spectra, from the decay rates of muons and charged pions tabulated in: /src/tables/hazma tables/electron.txt The Python script we have used to compute these tables is available at: /scripts/hazma scripts/hazma tables.py Then, we add up all the contributions to the primary photon and electron spectra to obtain the total secondary spectra. Primary π ± → e ± µ ± → e ± total (old) total (new) There are several drawbacks to this methodology. First, when FSR occurs, there is a loss of center of mass energy into photons, which should alter the charged particles spectra before decays, even if at second order in the fine structure constant. That difficulty arises because the Hazma code convolves initial spectra with analytical transfer functions and does not perform a MCMC statistical analysis of competing processes. Second, and this is the main difficulty with the QCD phenomenology, it is not clear what is precisely the QCD scale Λ QCD (between µ π 0 and 300 MeV) that separates direct pion emission from single quarks and gluons emission that hadronize into jets. This is the reason why we do not fix a Λ QCD at some value, but rather provide both hadronization possibilities. The spectra should interpolate between the two limits when the center of mass energy explores the QCD range. As we can not decide (outside of a chosen model) what are the number of degrees of freedom and their nature in the QCD energy range, this hadronization method is limited to the program BlackHawk inst, because we need this information to compute the Page evolution factors inside BlackHawk tot. The PBH masses at stake are M PBH 10 14 g, thus we can consider that these PBHs are living indefinitely at our time scale -the age of the universe -and that BlackHawk inst is enough for PBH studies linked to photon or electron emission in this scheme.
We show in Fig. 4 (respectively Fig. 5) the total photon (respectively electron) spectrum computed with the old extrapolated version of BlackHawk, compared to the new spectra computed with the method of [12] for a BH of mass 5.3 × 10 14 g.
We point out that this new feature is dedicated to the low energy photon and electron emission, and does not for now provide low energy spectra for neutrinos, which however are only related to sub-dominant (but complementary) constraints in the considered PBH mass range [40,[43][44][45]. We stress that the new spectra may alter the constraints of previous gamma ray and electron-positron studies from Hawking evaporation in the narrow mass range 5 × 10 14 g M PBH 10 16 g, where PBHs can only represent a negligible fraction of DM anyway [3,4].

III. PARAMETERS
In this Section we describe the new parameters.txt file resulting from the modifications and add-ons described in Section II, as well as a simplification of the presentation of the file by a removal of "hard-coded" parameters such as the dimensions of the numerical tables. These parameters are read by new routines in new infos.txt files. On another hand, as we have added new metrics, some parameters changed their name from * a (Kerr case) to * param in order to be generalized. As a consequence, the extended distributions for the BH spin a * have been adapted to the BH charge Q, while for polymerized and higher dimensional BHs there is only one value of the secondary parameter allowed for each run: (a number →) param number = 1. We focus here only in the new or modified parameters which are the interesting entries for BlackHawk users. For a complete description of the new and modified routines see Appendix A in this paper or the complete updated BlackHawk v2.0 manual [34]. Some bug corrections are listed in Appendix B. The parameters file/structure have several interesting new entries: • the new parameter metric switches between the Kerr metric (0), the polymerized metric (1), the Reissner-Nordström metric (2) and the higher-dimensional metric (3); • for the new BH metrics, we have added the parameters Qmin and Qmax (dimensionless -positive -BH charge between 0 and 1), epsilon LQG and a0 LQG (polymerized metric, being the polymerization parameters between 0 and 1 and a 0 the minimal dimensionless area), and finally n and M star (higher dimensional metric, n is the number of extra dimensions and M * the rescaled Planck mass); • tmin manual: as described in Section II C, there is now a possibility to choose between a manual BH formation time t form (tmin manual = 1, tmin = t form ) and the fidutial t form given by Eq. (9) (tmin manual = 0), the latter being only relevant for monochromatic BH distributions and radiation domination at formation; • nb final times has been removed as it is now fixed by the integration procedure; • nb gamma spins describes how many particle spins are available for some BH metric, and nb gamma fits encodes the number of fitting parameters used to extend the greybody factor tables to low and high energy; • add DM, m DM, spin DM and dof DM are new parameters linked to the DM emission described in Section II A. To add DM emission one has to set add DM = 1 (otherwise let it to 0), then one can choose the DM mass m DM in GeV, the DM spin s DM ∈ {0, 1, 2, 0.5, 1.5} and the DM number of dof (be carefull that only 1 dof has been taken into account when computing the new f, g tables with DM emission); • BH remnant switches between the usual total evaporation at Planck mass (BH remnant = 0) and an interrupted evaporation at M remnant ≥ M P (BH remnant = 1); • the parameter hadronization choice can now take the value 3, corresponding to the hadronization procedure described in Sec. II E. This choice forces the DM parameters to be used to compute the pion spectrum, thus the low energy hadronization procedure is at present incompatible with DM emission by BHs; • the Page factors tables parameters (Mmin fM, Mmax fM, nb fM masses and nb fM a → nb fM param), the greybody factor tables parameters (particle number, nb gamma a → nb gamma param, nb gamma x, a min → param min and a max → param max) and the hadronization parameters (Emin hadro, Emax hadro, nb init en, nb fin en, nb init part and nb fin part) do not appear anymore in the parameters.txt file as they are a priori fixed, they are now read in info files stored in respectively: and we recommend that they be not modified unless the tables are recomputed; • the following parameters changed their name: BHnumber → BH number, Enumber → E number, anumber → param number.

IV. OPTIMISATION
In this Section we briefly present a simple optimisation procedure that has been implemented inside BlackHawk in order to fasten the computations and diminish the disc memory occupied by the output files when some specific studies do not necessitate the full Hawking radiation spectra. Following the addition of the write * parameters which decide whether some particle spectrum is written down to the output, we added the tables compute primary[] and compute secondary[] which allow the code to skip some particle spectrum computations. We must warn the user that if some primary particle "computation" parameter is set to 0, then this particle will not participate in the secondary spectra as well. However, for studies focused on a single particle of some type, e.g. the graviton, DM or some secondary particle, then the other secondary particles can be "turned off". The BH time evolution will not be affected as it is determined by the numerical tables of the Page factors f, g.

V. CONCLUSION
We have described the new features available in the public code BlackHawk v2.0, together with interesting illustrations: the addition of dark matter (and dark radiation) emission, the massless spin 3/2 field greybody factors for the Schwarzschild metric, the possibility to keep a black hole remnant at the end of evaporation, the greybody factors associated with spherically symmetric and static black hole solutions more general than the Schwarzschild case and the careful computation of the low energy photon and electron spectra. We hope that this will open a new era of Hawking radiation studies with already promising results. An updated version of the complete manual is available on the arXiv [34] (v3) and the BlackHawk v2.0 code is available on HEPForge at https://blackhawk.hepforge.org/. double *initial energies, double *final energies, struct param *parameters) has been extended and used to generate the new header file: /src/hadro hazma.h compiled within the BlackHawk programs if HARDTABLES is defined.
We have modified the evolution routines in order to account for the new BH metrics, associated with so-called BH "secondary parameters" (the BH charge, the value of for polymerized BHs or the number of extra dimensions). In fact, all these parameters behave in the same manner as the spin. Some can even be treated as constants, like the number of extra dimensions.
Most routines also take into account the fact that the number of primary particles is now given by: particle number + grav + add DM

New routines
The new routines are: • int read fM infos(struct param *parameters) reads the f, g tables info in the new file: /src/tables/fM tables/infos.txt and fills the parameters Mmin fM, Mmax fM, nb fM masses and nb fM param; • int read gamma infos(struct param *parameters) reads the greybody factors tables info in the new files: /src/tables/gamma tables/infos.txt and fills the parameters particle number, nb gamma param, nb gamma x, param min and param max; • int read hadronization infos(struct param *parameters) reads the hadronization tables info in one of the new files: • void add FSR instantaneous(double **instantaneous primary spectra, double **instantaneous integrated hadronized spectra, double *energies, double *final energies, double *masses primary, struct param *parameters) computes the FSR of electrons, muons and charged pions as described in Sec. II E; we have been careful about the number of degrees of freedom, as for example the "electron" spectrum in BlackHawk in fact takes into account the e ± multiplicity, we have to divide the formulas of [85] by 2 because they also account for it a second time; • new void free *() functions have been defined to free arrays of different formats; • double exp adapt(double x) computes the quantity exp(x) − 1 in the case of small x, using the Taylor development of the exponential function up to 5th order.