Limits on primordial black holes detectability with Isatis: a BlackHawk tool

Primordial black holes (PBHs) are convenient candidates to explain the elusive dark matter (DM). However, years of constraints from various astronomical observations have constrained their abundance over a wide range of masses, leaving only a narrow window open at 1017g≲M≲1022\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{17}\,\mathrm{g} \lesssim M \lesssim 10^{22}\,$$\end{document}g for all DM in the form of PBHs. We reexamine this disputed window with a critical eye, interrogating the general hypotheses underlying the direct photon constraints. We review 4 levels of assumptions: (i) instrument characteristics, (ii) prediction of the (extra)galactic photon flux, (iii) statistical method of signal-to-data comparison and (iv) computation of the Hawking radiation rate. Thanks to Isatis, a new tool designed for the public Hawking radiation code BlackHawk, we first revisit the existing and prospective constraints on the PBH abundance and investigate the impact of assumptions (i)–(iv). We show that the constraints can vary by several orders of magnitude, advocating the necessity of a reduction of the theoretical sources of uncertainties. Second, we consider an “ideal” instrument and we demonstrate that the PBH DM scenario can only be constrained by the direct photon Hawking radiation phenomenon below Mmax∼1020\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{\mathrm{max}} \sim 10^{20}\,$$\end{document}g. The upper part of the mass window should therefore be closed by other means.


Introduction
Primordial black holes (PBHs) are early universe objects predicted by lots of inflation models; they could result from large overdensities collapse [1] as well as more exotic events such as phase transitions or topological defect collapse (see the reviews [2,3] and references therein or the lecture notes [4]). They are not the outcome of star collapse, hence their mass can span a wide range from the Planck mass M Pl 10 −5 g up to "stupendously large" values M ∼ 10 50 g [5]. One most interesting aspect of PBHs is that they can explain all or part of the missing dark matter (DM) density in the universe [2][3][4]). As such, the constraint on their abundance translates into a constraint on the fraction of DM, f PBH ≡ PBH / DM , they can represent.
Because all BHs evaporate due to Hawking radiation [6,7], losing their mass in a time related to their initial mass by τ ∼ M 3 init , their contribution to DM today implies that their initial mass was more than a few times 10 14 g. Reference [8] underlines that accretion may overcome Hawking radiation during the radiation domination era, causing PBH to grow during this period, so that the mass of PBHs evaporating just now is reduced from M 10 15 g to M ∼ 10 14 g. This justifies that we safely show the constraints down to M ≥ M min = 10 14 g. The emission of Hawking radiation by PBHs, in the form of all kinds of particles and particularly photons, gives rise to observational constraints. PBHs with mass M ∼ 10 20 g emit in the keV band; as the mean energy of Hawking radiation is inversely proportional to the PBH mass E ∝ 1/M, heavier PBHs M ≥ M max = 10 20 g are very hard to see through their photon Hawking radiation, as will be shown in the following.
During the last years, numerous constraints have been established on the PBH abundance in this disputed mass range [M min , M max ], with at stake the answer to the question of whether light PBHs are indeed the elusive DM. What we will henceforth call direct photon constraints are set by comparing the PBH emission spectrum to the galactic center (GC), or any other target, photon flux [9][10][11], to the stacked isotropic X-ray or gamma-ray backgrounds (denoted as EXGB) [12][13][14][15][16] or to a combination of the two [17][18][19]. There also exist direct photon (and neutrino) constraints on the local rate of PBH final burst as they reach the end of their evaporation, but this is specific to a fine-tuned PBH mass and thus we do not consider them here (see the very complete article [20]). Indirect photon constraints are related e.g. to the change in the re-ionization history at the cosmic microwave background [21,22]. Other studies are based on the direct [23] or indirect detection of electrons-positrons [22,[24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]; and finally there are some papers that focus on neutrino constraints [25,[39][40][41]. Most of the more recent studies now take into account the effect of having a nonzero PBH spin, which enhances the photon emission and thus results in more stringent constraints for Kerr (rotating) PBHs [9,13,18,19,25,28,32,37,38,41], or an extended PBH mass distribution, that should be more realistic regarding the channels of PBH formation [13,15,16,18,24,25,29,32,41]. Exotic studies have also derived constraints on non-standard (other than Kerr) PBHs [42][43][44] or PBH-DM mixed scenarios [45][46][47][48]. To do so, recent constraints make use of the public code BlackHawk [49,50] 1 to compute the precise Hawking radiation spectra, and we are using the latest version 2.1 of this code in our analysis. Photon constraints are handy thanks to several simplifications: photons are massless particles, thus any PBH mass is compatible with the emission of photons with the corresponding mean energy; 1 BlackHawk is available at https://blackhawk.hepforge.org. photons travel in straight lines, so their detection does not rely on complex propagation models; and the detection of photons by all kinds of instruments is a mastered technique that does not suffer from large measure uncertainties and has already been performed accurately in a very wide range of energies. This is why we focus on photons in this study, leaving electrons (all simplifications do not hold) and neutrinos (only the last simplification does not hold) for future work. We further restrict the scope of our analysis to most minimalist Schwarzschild PBHs with a monochromatic mass distribution and unclustered spatial distribution. A discussion on more elaborate PBH models is postponed to the last Section of this paper.
As stated above, the fraction f PBH is already partially constrained in the M min − 10 16 g range, leaving small doubt on the possibility that PBHs can constitute all of DM in the low mass part of this window, with constraints robustly set to f PBH < 10 −5 . However, the high mass range 10 16 -10 22 g has been the subject of multiple studies that claimed exclusion of 100% PBH-DM up to 10 17 and even 10 18 g (e.g. EXGB constraints in the case of an extended mass distribution of spinning PBHs [13]). The aim of this article is to make quantitatively explicit the underlying assumptions and approximations, to contextualize the claimed limit validity and to underline the way they can (or cannot) be compared to each other. This is of utmost importance because not the whole mass window that is currently scrutinized for 100% PBH-DM is accessible to Hawking radiation constraints, as will be demonstrated; extended mass functions that are fitted to the available remaining window constrain back the fine-tuned PBH formation models (and thus the early universe conditions) used to derive them; and the robustness estimation of the limits can help design optical instruments that would be sensitive to the most extreme parameter space accessible to Hawking radiation measures. We regroup the assumptions into 4 categories used throughout this paper, all of which discussed in details in the following: (i) instrument characteristics, (ii) computation of the (extra)galactic photon fluxes, (iii) statistical treatment, (iv) computation of the Hawking radiation. This paper is organized as follows: Sect. 2 gives the required basic formulas for Hawking radiation; Sect. 3 proposes a nomenclature of type (iii) assumptions; Sect. 4 discusses the robustness of the current existing and prospective Hawking radiation constraints on PBHs and Sect. 5 determines the characteristics and limits of an instrument for high PBH mass Hawking radiation studies; we conclude in Sect. 6. Throughout the paper, we use natural units c = G = k B =h = 1, making the constants appear uniquely when dimensionality is unclear.

Basic formulas
Since it has been shown that BHs emit quasi-thermally all particles in the "spectrum of Nature" [6,7], people have tried to extract a PBH signal from the observational data, without success so far. This absence of signal has been interpreted as a constraint on the PBH abundance over a wide range of masses. For an initial mass from the Planck scale M Pl ≈ 10 −5 g to some M min = 10 14 g, PBHs would have already evaporated away by now, forbidding them to be a sizeable component of the DM today, but their Hawking radiation could have left imprints in the big bang nucleosynthesis element abundances or in the cosmic microwave background [2]. Above an initial mass of some M min , PBHs would still be around, filling the universe with all kind of radiation, from photons to hadronized jets.
The emission rate of a massless particle i with spin s i by a Schwarzschild BH is given by [7] where i is the "greybody factor" (GF) that encodes the probability that the emitted particle reaches spatial infinity away from the BH horizon. The GFs deviate from the pure blackbody spectrum, with a power-law suppression at lower energies that depends on the particle's spin. Those GFs were numerically computed in the public code BlackHawk that we use in this study. 2 The BH horizon temperature is related to the BH mass through so that when the mass decreases because of evaporation, the temperature increases. The emission rate of Eq. (1) is further cutoff at the particle rest mass E > μ i , which is trivial for massless photons. It should be noted that both the BH temperature T and the GFs i depend strongly on the BH metric (see e.g. [43,51] for the general formulas for spherically symmetric BH solutions).
Once the instantaneous emission rates Q i of all particles i are known, they can be integrated over to obtain the BH mass loss rate [52] which determines the BH mass evolution equation Estimates of the function {(M) can be found e.g. in [53]. Inside BlackHawk we use a full numerical result. Considering that {(M) is a constant, we obtain a relationship between the BH initial mass M init and its lifetime τ τ (M init ) ∼ 13.9 × 10 9 yr M init 5 × 10 14 g 3 . (5) A more detailed study of the evolution of BHs can be found in e.g. [54].
We would like to emphasize that the above calculations assume the Hawking formula (1) [6,7] and the GF computation inside BlackHawk, based upon semi-classical general relativity, hold throughout (most) of the BH history. They may break down only at the very end of evaporation t ∼ τ , where τ is the usual BH lifetime calculated above, when M ∼ M Pl and quantum gravity effects become relevant. This has no effect on the constraints under discussion. However, it has been claimed that the so-called "memory burden" of BHs [55] -their capacity to store information -can slow down the evaporation process to extremely low rates or even stop it, within timescales t < τ. If that were true, the evaporation constraints set by e.g. photon background measurement would be dramatically alleviated and PBHs of mass M 10 14 g could represent a fraction or all of DM [56]. Hereafter, we assume that the usual paradigm holds until This quasi-thermal emission, that we call "primary", is not the final output of Hawking radiation. Most particles in the Standard Model are not stable (weak gauge bosons, charged leptons) or cannot exist outside confined hadrons (gluons and quarks). The primary emission must be convolved with analytical or numerical branching ratios Br i→ j to obtain the "secondary" spectra of stable particles that can be detected in instruments (photons, electrons, neutrinos, protons and their antiparticles). The secondary spectrum for particle j is theñ In particular,Q γ is the rate of emission of direct photons from PBHs, meaning that no interaction with the interstellar medium (ISM) or with other astronomical objects has been considered (scattering, absorption...). Inside BlackHawk v2.1 we use the Python package Hazma [57], that is relevant for low energy hadronization and decays (E some GeV), to compute the branching ratios Br i→ j . 3 This package considers that pions are emitted as fundamental particles instead of single quarks and gluons, and subsequently decayed into photons and leptons. Hazma, which relies on analytical decay and final state radiation formulas, undoubtedly suffers from inherent approximations. We should also implement Bremsstrahlung effects for emitted charged particles [58] that may dominate at very low energy, around the keV scale. Using Hazma, Ref. [10] showed that as the emission of all the massive particles (the lightest being the electron) is exponentially suppressed for PBHs with mass M 10 17 g only primary photons, neutrinos and maybe gravitons are radiated. Secondary photons however dominate the spectrum for PBHs with lower masses, with a mean energy Ẽ well below the PBH temperature. The implementation of Hazma inside BlackHawk v2.0 was a necessary improvement compared to the extrapolation of the PYTHIA [59] and HERWIG [60] results used in the previous version BlackHawk v1.2. The PBH nature (beyond Schwarzschild) and mass distribution, as well as the computation of the GFs and branching ratios, belong to the type (iv) assumptions.
In the following two subsections, we give the basic formulas for (extra)galactic PBH photon flux computation, together with a discussion of related type (ii) assumptions. Both contributions are important because coincidentally, they are of the same order of magnitude at MeV energies for PBHs right in the relevant Hawking radiation window M ∼ 10 15 −10 18 g (a fact that was first noted by Ref. [61]). Low energy photons principally come from the redshifted extragalactic spectrum while high energy photons originate from the galaxy [17][18][19].

Galactic flux
If PBHs represent some fraction of DM, then they are present in the Milky way halo with some spatial distribution ρ gal (r ) (in galactocentric coordinates) which we take as spherical for simplicity. The photons they emit through Hawking radiation propagate in straight lines from their origin PBH to Earth. Thus, the flux of GC PBH photons received by an instrument on Earth per unit energy, time, surface and solid angle is 4 where with the field of view considered and A gal ≡ M is the normalisation constant for a monochromatic mass M distribution of Schwarzschild PBHs. Note that this formula could be adapted for any compact source other than the GC with the relevant surface factor J integrated over the volume of that source. We see that the flux depends on the mass distribution of DM in the Milky Way halo and on the precise location R 0 4 These formulas are given in some form in e.g. [9][10][11][17][18][19]. of the Solar System in this halo, contained in the definition of the line of sight (LOS) [19].

Extragalactic flux
PBHs form in the early universe, with a formation time t form related to their initial mass by [2] t form ∼ 10 −23 s M init 10 15 g , assuming radiation domination. Hence, they emit particles through Hawking radiation during all cosmological eras until today, and this continuous flux piles up with a dilution factor a(t) = 1 + z(t). The flux on an instrument today per unit energy, time, surface and solid angle is given by 5 with t min = t CMB for photons, because the universe is opaque to light before that, and t max = min(τ, t today ). The normalization constant for a monochromatic mass M distribution of Schwarzschild PBHs is A egal ≡ M/ρ DM where ρ DM ≡ DM ρ crit is the global density of DM today. We see that the flux depends on the value of the cosmological DM density fraction, but also on the redshift history of the universe; we have further neglected the optical depth of the ISM (line absorption of light) apart from the cut-off at the CMB time in integral (10).

A nomenclature of constraints
In the previous Section, we have computed the flux of photons emitted directly by PBHs, with a galactic and an extragalactic contribution. Now, we detail how the fluxes (7) and (10) are compared to different sets of data to derive PBH abundance constraints. This discussion is quite trivial and cumbersome but we think that it is of central importance to recall the basics. Altogether, this Section thus deals with type (iii) assumptions.

Existing data
The first and very classical method to obtain a constraint on the PBH abundance from Hawking radiation consists of directly comparing the predicted flux for some PBH mass and abundance -we recall that we focus only on a monochromatic distribution of Schwarzschild PBHs -and some set of data measured by an instrument. Photon data are either presented in a differential (energy/events per unit time) or integrated (total energy/events) form, as a function of photon energy. Both should be equivalent as we do not expect the PBH signal to vary during the time of observation. Then, to constrain a signal, one can compare the spectrum to each energy bin of the data set, asking that the PBH photon flux does not overrun the measured flux. One can also consider the whole instrument energy band and compare the integrated quantities. This is a model independent method, but there are already two ways to compare the signal to the data, and there remains a choice to make as for the confidence level (CL) we chose to exclude PBHs (data, data +σ , data +2σ , ...). One can also fit the data with some function (from a simple segmented function to a complex analytical formula), motivated by a physical interpretation or not. This is already a model dependent approximation, as it infers data in unmeasured energy bins from data in measured ones. One can finally go one step further, by assuming that some fraction of the measured signal comes from astronomical sources. Thus, there is even less available parameter space for the PBH abundance. This is highly model dependent, and contains a hidden feedback loop: most astronomical backgrounds are estimated (calibrated) thanks to the data. Furthermore, the data are often "cleaned" using catalogs of identified point-sources to obtain the diffuse components. This introduces a bias that we have used in the present analysis: if PBHs are highly clustered in the Galaxy, they would resemble point-like sources and be cleaned away in the diffuse component search procedure. Constraints for clustered PBHs require a dedicated treatment (see e.g. [62,63]).

Prospective instruments
When data is not yet available in some energy range, one can put a conservative constraint on the PBH signal by saying that if the prospective instrument designed to explore this particular energy range is built and measures nothing, then this means that the signal is below the sensitivity of the instrument, with some confidence level. When data is already available, it is very complicated to decide what to do. In fact, the more precise instrument to come can totally revolutionize the measures due to an error of appreciation of the functioning of the previous ones. An independent conservative constraint would be to predict what sensitivity to a PBH signal some instrument can in principle reach, assuming that all the signal comes from PBHs or that some (model-dependent) background is to subtract beforehand. This is quite radical and nobody expects that all the measures taken by long-time working instruments are to be thrown away. A more reasonable method is to build a model of background, and to use it to estimate what would be the prominence of the PBH signal "above" the background (signal-to-noise ratio -SNRmethod or χ 2 method). This also contains a feedback loop: Fig. 1 The different constraint methods. We have restricted the tree to methods we have encountered in the literature: for existing instruments, "direct data" means that the instrument data points are used, and "fitted data" that a fit running through the data points is chosen; for prospective instruments, "direct bckg" means that the data points used as agnostic background, "fitted bckg" that some fit is used and "full sens" that preexisting data are ignored and the new detector maximal sensitivity is used. The designation we use in the text are: "type 1" methods M 1 with sub-methods M 1 1 and M 2 1 for existing instruments; "type 2" methods M 2 with sub-methods M 1 2 , M 2 2 and M 3 2 for prospective instruments background models are often calibrated thanks to the older instrument data, thus a deviation from the expected background compensated by a PBH signal would appear as no signal at all. Once again, the PBH signal and the background can be compared over the whole energy range or inside each energy bin separately.

Summary
Hence, we see that even the simplest data-signal comparison method contains non-trivial features that complicate the comparison between constraints set with different choices. This is even worse for the prospective instrument methods as the dependency to the background modelling is complex. We do not intend to build a hierarchy of the methods chosen by different authors but just to quantitatively highlight their differences. Reference [10] is the first, to our knowledge, to compare the PBH constraints from different prospective instruments with rigorously the same statistical method. One must also check that instrument characteristics and theoretical choices for the PBH flux calculation are coherent from one study to another. We have built a tree in Fig. 1 to summarize the methods described above, restricted to the ones used in the literature cited in the Introduction. Please refer to this tree for all the subsequent abbreviated method nomenclature.

A numerical tool: the Isatis program
For the sake of this study, we have designed Isatis, a numerical public tool that relies on the BlackHawk PBH spectra to compute the constraints with a controlled set of assumptions. The code is presented in details in the Appendix A. The idea is to use BlackHawk to obtain the direct PBH photon spectrum, and then to derive the constraint on the PBH abundance for a list of optical instruments, existing or prospective. All constraint assumption types (i)-(iv) listed in the Introduction, with different statistical methods from the nomenclature given in Fig. 1 available, can be modified. Hence, the quantitative individual impact of the assumptions on the PBH "direct" photon constraints can be investigated. This allows to compare consistently the constraints from the literature. Identifying the dominant parameters further gives an insight on where to look for improvements in the constraints. In this Section, we revisit the literature constraints from GC and EXGB PBH photon fluxes and examine type (iii) assumption impact.
To obtain the PBH constraints, we must use method M 1 . We show the results for the submethod M 1 1 in Fig. 3: direct comparison of the signal to the data requiring that in each energy bin, where the fluxes d gal/egal /dE are given by Eqs. (7) and (10), E X low/up are the bounds of the energy bins of instrument X and d X gal/egal /dE (resp. X ) is the photon flux (resp. error bar) measured by X and shown in Fig. 2. CL is the confidence level, translated in the number of error bars considered above or below the central value of the data points. We have chosen a NFW DM profile [78] with the "convenient" set of parameters from [79] for the galactic flux, and the standard radiation domination from CMB (t ≈ 1.2 × 10 13 s) to today (t ≈ 4.4 × 10 17 s) for the extragalactic flux, with DM given by Planck [80]. In Fig. 3, we have performed an empirical power-law fit to the most stringent constraints, with parameters for the GC, which overestimates the results at ∼ 5 × 10 15 g and for the EXGB, which can be compared to [3].
As an example, let us now examine the impact of the CL parameter, which takes into account the error bars of the photon fluxes. This is a type (iii) assumption. To do so, we repeat the analysis of Fig. 3 with CL = {−1, 1, 2} (lower error bar, upper error bar and twice the upper error bar), and present in Table 1 the maximum relative discrepancy in f PBH , as compared to the case CL = 0 (central values). Obviously, instruments with large error bars are the most affected: looking at Eq. (11), we easily deduce that an increase of the instrument flux by a factor α results in a constraint relieved by the same factor α. This is directly shown by the linearity between CL and the relative increase in f PBH in Table 1 (slightly broken for asymmetric error bars). We observe that the CL parameter can have an impact of up to a factor of a few on f PBH for instruments with very large error bars.
The impact would be exactly of the same kind if the flux data are replaced by some fitting function: the relative discrepancy in the constraint would be directly proportional to the relative discrepancy between the data and the fit. On the other hand, an increase of the PBH signal by a factor α results in a constraint more stringent by the same factor α.

Prospective instruments
A large set of prospective instruments, among which AdEPT, AMEGO, ASTROGRAM (results for the AS-ASTROGRAM design are new to this work), GECCO, GRAMS, MAST, PANGU and XGIS-THESEUS, are designed to explore the MeV energy range more accurately than previously done with COMPTEL and EGRET, which will improve both the description of the EXGB and the diffuse GC emission. This energy scale is crucial to set constraints on the fraction of DM f PBH in the disputed mass range 10 15 -10 18 g. Data is already available in this energy band, so new instruments will only reduce the error bars. The total number of photons -background plus PBH signal -observed by a prospective instrument X is given by [10,11,19]   . For completeness, we show the following background models: [69] (dashed), [70] (dot-dashed) and [17] (dotted). Right: isotropic extragalactic photon flux measured by HEAO+balloon [71], COMPTEL [72], EGRET [73] and Fermi-LAT [74] (model A). For completeness, we show the following background models: [71] (dashed, black), [75] (dashed, grey), [76] (dot-dashed), [14] (dotted), [17] Fig. 2. For comparison, we show the limits set by Ref. [9] for INTEGRAL (dashed) and Ref. [10] for COMPTEL (dot-dashed), stressing the fact that the statistical method is not the same as here. Right: PBH constraints from EXGB observations of the 4 instruments of Fig. 2.
For comparison, we show the limits set by Ref. [12] (dashed) and Ref.
where T obs is the duration of observation, is the solid angle field of view (fov), A eff is the effective area (which defines the energy band probed, see Fig. 4, left panel) and W (E, E ) is a window function accounting for the finite energy resolution of the instruments. Following the literature, we take this to be a Gaussian [83]. 8 All the useful references concerning these instruments are listed in Table 3 in Appendix A. Following method M 2 , the PBH constraint is obtained by requiring that the SNR stays below some detection threshold where we have separated the PBH and the background contributions to the total photon count. An alternative method would consist in demanding that this SNR is not attained in any energy bin of the background. We have checked with Isatis that this results in a relative change of the constraint f PBH up to a factor of a few.
Most of the difficulty resides in the choice of background. As discussed above, most backgrounds are calibrated to the data and thus automatically contain feedback loops. In Fig. 4 (right panel), we show the PBH constraints obtained by using the same fitted background as Ref. [10] (see [70] for details), i.e. method M 2 2 ; with SNR = 5, T obs = 10 yr and = 5 o × 5 o . The other parameters follow Sect. 4.1. While not reproduced here, we checked that we obtain results very close to that of Ref. [10] for AdEPT, GECCO and GRAMS (see also [84]), Ref. [17] for AMEGO and Ref. [19] for XGIS-THESEUS. There are discrepancies with Ref. [11] for GECCO and Ref. [43] for AMEGO, probably linked to the different statistical treatment. Stranger are the discrepancies we find with respect to Ref. [10] for AMEGO, eASTRO-GRAM, MAST and PANGU: the constraints obtained by [10]  Relative extremal variation in f PBH when using the same background as [19] or the central values of the data points compared to the background of [70] for the instruments of Fig. 3

Instrument
Ref. [ As a second example, we explore the effect of the background choice, that is also a type (iii) assumption. In Table 2, we show the extremal relative change in f PBH when choosing the same background as Ref. [19] and the agnostic background consisting of the central values of the data points (i.e. method M 1 2 ), as compared to the background of [10]. We observe that the choice of background, for the 3 examples shown, can have an impact up to a factor of a few on f PBH in the case of the instrument XGIS-THESEUS.

Reverse engineering the Hawking radiation constraints
In the previous Section, we have revisited existing constraints and exposed their robustness relative to some type (iii) assumptions, namely the background choice and the statistical method of comparison between the PBH signal and data. In this Section, we examine thoroughly all the other assumptions. Indeed, changing the instrument characteristics at will inside Isatis gives an unprecedented access to a reverse procedure: what would be the capabilities of an "ideal" instrument if we could fix arbitrarily all the parameters? Answering this first question carefully could certainly give hints towards design choices for future instruments. But we can go deeper: whatever be the capabilities of the instruments, is there a limit to the PBH mass range we can constrain with direct photon PBH constraints? In other words, can we close the window for all DM into PBHs still open between ∼ 10 17 -10 22 g? If the answer to this second question is yes,  [10,19,81,82]. We have truncated the effective area of MAST as we are not concerned in > GeV photons for the PBH mass range probed here. Right: PBH constraints from GC+EXGB observations using the method M 2 2 with the background of [70] this will certainly weigh in the favor of dedicated instruments designs. If it is no, theoretical efforts will need to be pursued to use other kinds of constraints in the (to be determined) remaining window. As already stated, for a given PBH mass M, the constraint f PBH (M) for an "ideal" prospective instrument is impacted by several assumptions regrouped in 4: (i) instrument characteristics, (ii) uncertainties on the (extra)galactic fluxes, (iii) different statistical treatments and (iv) uncertainties on Hawking radiation. Each of those can then be decomposed into several contributions, that we review in detail below. We define a test instrument with a given set of technical characteristics, along with fidutial parameters for the (extra)galactic fluxes and standard assumptions for the Hawking radiation, and we fix the statistical method: (i) observation time T obs = 10 yr, fov = 5 o × 5 o around the GC, relative energy resolution (E) = 1%, constant effective area A eff (E) = 10 3 cm 2 for E = 1 keV − 1 GeV, (ii) standard radiation domination from the CMB to today for the extragalactic flux and DM from Planck [80], NFW profile for the Milky Way halo with the "convenient" parameters of [79] resulting in J gal (cf. Sect. 4.1), (iii) statistical method M 2 2 (central values of the data points as background) which considered with point (ii) gives d bckg /dE, and SNR = 5, (iv) Schwarzschild PBHs with a monochromatic distribution, Hawking primary spectra Q i computed by BlackHawk with the low-energy Hazma branching ratios Br i→ j (E, E ) for the secondary spectra. 9 9 In particular, we do assume that the semi-classical Hawking radiation holds on throughout most of the PBH lifetime, see discussion in Sect. 2.
These characteristics result in a fidutial constraint that we denote by f PBH (M). Any modification of some assumption results in a different constraint f PBH (M) that we parametrize through where the α's are the detailed contributions from types (i)-(iv) assumptions. We detail the mathematical form of these in the following, and give quantitative results for the modified f PBH thanks to Isatis.

Instruments parameters
A direct look at Eqs. (14) and (15) shows that and as the number of photons captured for both the background and the signal is directly proportional to these quantities. In the limit in which the fov covers only the desired target, we also have Energy resolutions are overall small and might get smaller in the future as detection techniques improve, so that the window function in Eq. (14) tends towards a Dirac distribution, reducing the pollution from neighbouring energy bins. Thus, changing the energy resolution should have a negligible impact α( ) ∼ 1. All these estimations are confronted to

"Test" instruments were implemented inside
Isatis to obtain these constraints quantitative results from Isatis in Fig. 5. First, we observe that the fidutial constraint extends up to ∼ 10 19 g, because the energy coverage of that "ideal" instrument goes down to the keV scale. Second, we note that the constraint is close to a power-law of the PBH mass f PBH (M) ∼ (M (g)/8× 10 18 ) 2 , a feature that comes from four facts: the primary emission peaks at an energy that is proportional to the PBH mass; the number density of PBH is inversely proportional to their mass; the background fit of [70] is approximately a powerlaw; and the fidutial effective area is a constant. The effect of the secondary spectra is somewhat visible in the slope breaking below 10 15 g. The global variations of f PBH compared to f PBH from the modification of the various parameters has precisely the behaviour expected from Eqs. (17), (18) and (19). We verify that the energy resolution has no visible impact.

Flux parameters
As the column density of the target impacts only N PBH , we can estimate The precise impact is obviously more complex depending on the predominance of the target flux over the diffuse extragalactic flux, which also depends on the energy considered.
As seen before, we expect the target flux to dominate at high energies, constraining low PBH masses, and the redshifted isotropic flux to dominate at low energies, constraining high PBH masses. In Fig. 6 we show the extremal change in f PBH for different galactic DM profiles: • the generalized Navarro-Frenk-White (NFW) profile defined by [78,85] ρ DM (r ) = ρ c r c r where the classical NFW is obtained for γ = 1 • the Einasto profile defined by [85,86] Both of these profiles depend on 3 independent parameters: a characteristic density ρ c , a characteristic radius r c and a power γ . Even with modern measures, these parameters suffer from large uncertainties. We use the 68% CL parameters of [85] (Tables II and III corresponding to two baryonic models B1 and B2) to be as general as possible in Fig. 6. 10 Spanning the whole 68% CL range of the parameters results in constraints f PBH varying by 2 (resp. 1) orders of magnitude inside the same profile and baryonic model for NFW (resp. Einasto) profile; by ∼ 50% (resp. ∼ 1%) from one baryonic model to the other (B1 or B2) for NFW (resp. Einasto) profile; and by 20% (resp. 30%) from one profile to the other (NFW or Einasto) for model B1 (resp. B2). The largest variations are obtained with the generalized NFW profile. We are far from having explored the complete diversity of the DM halo profiles existing in the literature, but we can already predict that α(gal) ∼ 10 −1 -10 1 . We do not consider targets other than the GC. Even if the density factor is determined more precisely in e.g. M31 or DM concentrated dSphs like Draco because we observe them as a whole and are not sensitive to their halo profile, it is much smaller than that of the GC and results in less stringent constraints (10 −2 factor reduction for M31 and Draco [10]). A modified redshift history of the universe or DM density may have a small impact on the time stacked isotropic spectrum. We have checked that this is negligible regarding the small uncertainties on the cosmological parameters for the recent universe (after the CMB), resulting in α(egal) ∼ 1. 10 References [10,85] use the incredibly precise radius of the Sun orbit R = 8.122 ± 0.031 kpc recently obtained by the measurement of the orbit of S2 [87], which differs from the older "convenient" value of [79]. Reference [10] further uses the central values of the NFW profile of Table III

Statistical treatment
As the photon background only affects N bckg , we have already concluded in Sect. 4.2 that within method M 2 . The different choices of backgrounds for prospective instruments in Table 2 result in α(bckg) ∼0.5-3 at most. This impact is completely equivalent to that of changing the CL parameter inside method M 1 1 , as shown by Table 1. On the other hand, the SNR is directly proportional to f PBH , thus As this paper was finalized, the author became aware of Ref. [88], where the constraints from the EXGB are revisited with data from instruments different from those presented in Fig. 2. The contribution from the GC included, as well as Hawking radiated electron-positron annihilation. Reference [88] takes into account astrophysical contribution to the diffuse photon flux in order to obtain more stringent constraints on the PBH abundance, thus using what we denoted as method M 2 1 . They obtain constraints more stringent by 2 orders of magnitude in the mass range 10 15 -10 17 g.

PBH Hawking radiation
Last but not least, there are uncertainties related to the very computation of the PBH Hawking radiation spectra. These only affect N PBH in Eq. (15), thus we predict and The uncertainties on the primary emission rates for Schwarzschild PBHs are directly related to the tables contained inside BlackHawk, which have been confronted to the literature with an accuracy of less than a percent, such that α(Q γ ) ∼ 1.
There could be large uncertainties in the branching ratios computed by the particle physics codes. At very high energy E TeV BlackHawk relies on HDMSpectra [89], in the documentation of which is explained that the precise treatment of the electroweak cascades can alter the branching ratios into photons by several orders of magnitude compared to PYTHIA. This has not yet been confronted to accelerator data. At LHC energies E ∼ GeV − TeV, BlackHawk relies on PYTHIA [59] or HERWIG [60] with good correspondence to the data, even with the QCD uncertainties [90]. 11 At low energies E GeV, BlackHawk relies on Hazma [57], which is also based on accelerator data, but does not take into account e.g. the Bremsstrahlung radiation of charged particles specific to PBH radiation [58] which should dominate at the keV scale. The choice of the QCD scale QCD at which pions are emitted as primary particles and of the dynamical rest masses of the quarks and gluon introduces further uncertainties (see [92] for a recent discussion). Ref. [89] has shown that the branching ratios for "low final energies" E/E 10 −4 in Eq. (6) as computed with PYTHIA (and thus HDMSpectra) suffer from order-of-magnitude uncertainties linked to the difficult tracing of electroweak cascades on such stretched scales. Hence, we can conclude that even with very precise primary spectra, the secondary spectra (integrated over primary energies from all scales) are associated with an order-of-magnitude possible variation from all these codes α(Br i→γ ) ∼ 10 −1 − 10 1 .

Summary
Summarizing the results from this Section, we see that compared to our fidutial case f PBH , the value of the PBH constraint for all masses can vary within each set of assumptions: 11 Updated estimates of the QCD uncertainties will be released soon [91].
Overall, we obtain with all the sources of uncertainties and if we restrict ourselves to the fidutial instrumental parameters, which is still, as we say in french, une sacrée fourchette -a very large range, as it runs over 11 orders of magnitude for α tot and 5 orders of magnitude for α ii)− iv) . The extremal constraints are schematically shown on Fig. 7. The 100% PBH-DM scenario is excluded up to M ∼ 10 20 g in the most favorable case but the open window goes down to M ∼ 10 17 g in the less favorable one. Of course, the extremal values of the (independent) uncertainties (ii)-(iv) are presumably not attained at the same time, but nevertheless we see that PBH constraints from direct Hawking radiation of photons are still highly model dependent due to several sources of uncertainties. We see that modifying the prospective instrument characteristics allows to close the M 10 20 g window, but we also note that the slope of f PBH is nearly vertical in this region; we discuss this last feature below.

Very low energies
We observe on Fig. 7 that the constraint f PBH saturates at high PBH masses. This is due to the exponential cutoff of Q γ and thus of d PBH /dE at energy E E max ∼ T PBH ≈ 1 keV(M/10 19 g) −1 for Schwarzschild PBHs. Hence, as our description of the background and of the secondary photon spectrum is limited to this keV lower energy bound in this study, we cannot obtain f PBH at M 10 20 g. Increasing the capabilities of our prospective "ideal" instrument would only make the constraint steeper at M ∼ 10 20 g. Thus, pushing the PBH constraint on direct photon emission up to lower energies and thus higher PBH masses, with the aim of constraining the whole (yet) open window M ∼ 10 18 −10 22 g for all DM into PBHs, would require both precise background description at E ∼ 10−1000 eV and precise particle physics codes to compute electroweak showers down to the eV energy range. We can suppose that the primary spectrum dominates below the keV scale as the emission of all particles except for photons, neutrinos and gravitons is exponentially suppressed at M 10 17 g. However, this energy range falls down right in the (very) far UV band (λ ∼ 10 − 1000 Å). To our knowledge, this band is not yet covered within a precise all-sky (resp. GC) survey that would give an access to the isotropic (resp. GC) component of the background. The GALEX [93] (and future UVEX [94]) instruments are sensitive only starting above λ ∼ 1000 Å. Suppose that we extrapolate the "ideal" instrument capabilities, the PBH Hawking radiation rates and the GC+EXRB background down to the eV energy scale: d bckg /dE ∝ E −2 and A eff (E) = A eff (E) at E = 10−1000 eV. For the fidutial values of the parameters (i)-(iv) we obtain f PBH ∼ 10 6 at M = 10 22 g. Closing the window for all DM into PBHs with direct photons would then require to increase the capabilities of our "ideal" instrument by a factor 10 6 (or reduce f PBH by a factor 10 −6 ), which is unrealistic as α min i) = 8 × 10 −4 is already technically challenging. Plus, we expect the spectrum of light at these wavelengths to be overcrowded by astronomical sources (see e.g. Fig. 2 of [94]), further increasing the amelioration factor needed to obtain f PBH ∼ 1 at M = 10 22 g. Thus, we conclude that the window for all DM into PBHs cannot be closed by direct photon detection from PBH evaporation. 12 Complementary constraints, relying on complex modelling, are thus needed: dynamical capture by or microlensing of stars (see the review by [96]), rare collisions with solid objects [97] or GRB lensing [98].

Conclusion
In this paper, we have computed the photon emission spectrum by Hawking radiation of monochromatic Schwarzschild primordial black holes. We restricted ourselves to the direct photons resulting from the primary emission of photons and the secondary emission resulting from electroweak interactions or decay of other primary particles. We obtained the flux on Earth of photons coming from the galactic center plus the isotropic redshifted emission from past ages. We identified all the assumptions underlying the computation of this flux and ended up with a clear nomenclature of the constraining methods M j i , classifying existing and prospective constraints. Then, we used Isatis, a new public tool available inside the Hawking radiation code BlackHawk and described in the Appendix, to examine quantitatively how the constraints presented in the literature depend on those assumptions, regrouped in 4 categories: (i) instrument characteristics, (ii) (extra)galactic flux description, (iii) statistical method used and (iv) Hawking radiation uncertainties. We have shown that, for a given instrument design (i), the constraint on the DM fraction into PBHs f PBH (M) can span up to 5 orders of magnitude due to uncertainties (ii)-(iv). This has a direct impact on the size of the current and prospective available window for the 100% PBH-DM scenario and questions the robustness of the current constraints. We underline the necessity of more precise galactic photon background and Hawking radiation determination, especially at low energies where high mass PBH constraints could close the window. We emphasize that the very validity of the semi-classical computation of Hawking radiation assumed throughout the paper, that could be included as one of the (iv) uncertainties, is a fundamental basis of the whole work. Finally, we have examined to what extent a prospective "ideal" instrument could constrain f PBH < 1 from (currently) M ∼ 10 18 g up to M max ∼ 10 20 g by "reverse engineering" Isatis. We conclude that even if (i) instrument characteristics allow for another 5 orders of magnitude variation of f PBH (M), above M max , direct photon constraints are not effective anymore so that complementary constraints should be developed instead.
There are numerous ways to build up on this study and improve Isatis. First, we have only considered Schwarzschild PBHs, but some early matter domination models predict PBHs born with high spin, which would increase their photon yield and result in more stringent constraints. Non-standard BH solutions, e.g. derived from effective loop quantum polymerization, also predict very different photon spectra. Deviations from the semi-classical Hawking radiation computation at mass scales greater than the Planck scale could also be explored. Second, we have considered a monochromatic mass distribution, which is obviously unrealistic and could be refined within a specific PBH formation model. The main effect is to "spread" the constraints towards higher and lower PBH masses due to the high and low mass tails of the distribution under consideration. Third, we have computed the direct photon spectrum from instantaneous Hawking radiation, but we could go further and obtain the indirect photon spectrum after model dependent interactions of direct photons or electrons/positrons with the interstellar medium. One last possibility is to constrain PBHs through their (in)direct electron/positron, neutrino, graviton or (putative) DM particle yield. This is left for future work, but a thorough list of references can be found on the BlackHawk website (for those who rely on the code) or in the latest review [4]. Finally, we have assumed that PBHs are not clustered when examining the Galactic constraints, but are rather smoothly distributed following the halo density function. Clustered PBHs would appear as point-like sources and would require a suited search.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All data generated or analysed during this study are included in this published article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Isatis: a public tool to compute consistent PBH constraints
This Appendix is dedicated to the presentation of Isatis, a new public numerical tool to compute PBH Hawking radiation constraints with a controlled set of assumptions.

Technical aspects
Isatis is a public C tool that relies on the BlackHawk [49,50] code to compute PBH constraints. Be careful to download the last version of BlackHawk (v2.1 as we write this article) before using it. Isatis is available as one of the "add-ons" (/scripts folder) of BlackHawk and as such can be downloaded together with the main code on the same website: https://blackhawk.hepforge.org The Isatis main folder contains: • the main file Isatis_photons.c that defines all the Isatis routines, • a parameter file parameters.txt that defines the model parameters for the constraints (galactic and extragalactic fluxes assumptions, statistical comparison method), • a file README.txt that summarizes the basic steps in the use of Isatis, • a file plotting.py that can be used to plot the constraints obtained with Isatis, • a folder /constraints that contains all kinds of numerical tables used to compute the constraints (e.g. tabulated background fluxes or cross-sections), • a folder /BH_launcher that we explain below.

The BH_launcher program
The folder /BH_launcher contains another BlackHawk add-on: an automatic comprehensive launcher adapted to Isatis. With this program, you can launch several Black-Hawk runs in parallel in order to obtain a set of spectra for e.g. different PBH initial masses. The script is interactive so its use should be transparent, but a README.txt file is provided anyway. First, please ensure that BH_launcher.c contains the correct path to your BlackHawk version and that the BlackHawk programs BlackHawk_inst and BlackHawk_tot are themselves compiled. We advise to compile it with #define HARDTABLES to save time at execution. To obtain the executable BH_launcher.x type the command make in the terminal from the /BH_launcher folder. BH_launcher.x dos not require any argument so just run it with the command ./BH_launcher.x. The script launches BlackHawk the desired amount of times with the given parameters and generate several outputs: • the BlackHawk output is stored as usual inside the BlackHawk/results/ folder, • BlackHawk/nohup_*.txt files are generated for each run of BlackHawk, • a file BH_launcher/* is generated that contains the number of runs and the name of each run, that may be given as an argument to Isatis.
Launching many (tens of) parallel executions of BlackHawk can be memory consuming for both RAM and disk. We advise that you tweak the files write_*.txt in BlackHawk /src/tables in order to keep only the photon output that is needed by Isatis. With BH_launcher, you can choose the type (iv) assumptions, namely you can play with the capabilities of BlackHawk concerning non-standard PBHs, extended mass distributions, etc.

The Isatis program
The Isatis main program takes two arguments: a parameter file that contains all the assumptions of types (ii)-(iii) to use in the constraint determination and a file which contains the number of BlackHawk runs nb_runs = * that you want to compute PBH constraints from, and the names of those runs as a list. The latter can either be generated automatically by the side program BH_launcher or by hand, listing the results folders of interest that the user may already have generated with BlackHawk. The constraints computation requires that the secondary spectra of particles are computed, and that both programs BlackHawk_inst and BlackHawk_tot have been used on the same parameters file. Otherwise, errors are cast when trying to read nonexisting output files. Please ensure that parameters.txt contains the correct path to your version of BlackHawk.
To compile Isatis.c into Isatis.x just go into the /Isatis main folder and type make Isatis. Then, launch it with the parameter file * and the run list file ** as an argument with the command ./Isatis.x * **. This should generate a result file results_photons_*.txt containing the list of PBH constraints with each line being one of your BlackHawk runs (e.g. with some specific PBH initial mass), and each column the corresponding PBH constraint set by a specific (existing or prospective) optical instrument.
To summarize, a basic use of the Isatis program consists in (supposing everything has been compiled correctly): 1. launching BH_launcher.x to generate a set of Black Hawk spectra with some varying PBH parameters, 2. launching Isatis.x to compute the corresponding PBH photon constraints, 3. using plotting.py to visualize the constraints.

Instrument implementation
For existing instruments, we use the most recent data they have produced. Those data are presented in Fig. 2 together with the corresponding literature. The measured photon fluxes in the direction of the desired target (galactic center, isotropic flux) are tabulated in files constraints /photons/flux_*_**.txt. These are used to obtain the constraints labelled M 1 in the nomenclature of Fig. 1. For prospective instruments, instrument characteristics (type (i) assumptions) are directly hardcoded inside Isatis, with relevant literature listed in Table 3: • the time of observation T obs , • the fov , • the effective area as a function of photon energy A eff (E) (tabulated in files constraints/photons /Aeff_*_**.txt where * is the instrument name and ** the arXiv identifier), • the relative energy resolution as a function of photon energy (E).
We propose several (extra)galactic backgrounds inside Isatis, labeled as constraints/photons/back-ground_**.txt. These are used to obtain the constraints labelled M 2 in Fig. 1. All files contain * the name of the instrument and/or ** the arXiv identifier to help traceability.

For information: what's in the literature?
For transparency, we classify the PBH direct photons constraints from existing (Table 4) and prospective (Table 5) instruments that are present in the literature and listed in the Introduction. This is not an exhaustive list as it is restricted to the bounds usually shown in the recent review efforts. This classification could obviously be extended to older or to be published constraints. In this study, we have used method M 1 1 for prospective instruments.

Perspectives
Modification of Isatis should be straightforward. Adding a new instrument or a new background consists in providing the relevant *.txt flux or effective area file. New galactic profiles can be easily implemented, as well as targets different from the GC, thanks to the corresponding density factor J . The code is already adapted to extended mass functions with general expressions of the normalization constants A of Eqs. (7) and (10). BH_launcher is also handy to modify, as it works following a decision tree. Variation of parameters different from the mass from one run to the next consists in adding branches to that tree.
In the near future, we intend to produce neutrino and electron-positron tools as well, in the form of additional programs Isatis_neutrinos.c and Isatis_ electrons.c. The scheme will be the same as for Isatis_photons.c but adapted to constraints set for direct neutrinos and electrons-positrons. Indirect detection constraints are also a horizon of development.