1. Introduction
Several astrophysical and cosmological evidences point towards the existence of dark matter (DM), something that constitutes ∼26% of the total energy density of the Universe [1]. Despite its abundance, the nature of DM remains mysterious, as it has evaded all nongravitational direct and indirect detections thus far [2–4]. One of the most widely discussed DM candidates is primordial black hole (PBH) [5–9]. PBHs may have formed in the early universe from the collapse of overdensities resulting from quantum fluctuations [10–14] and by other mechanisms such as phase transitions [15–22]. The mass of PBHs can vary in a large range depending on the formation time. In general, the initial mass MPBH of a PBH should be close to the mass enclosed by the Hubble horizon at the formation time t, i.e. MPBH ∼ c3t/G ≃(t/10−23 s) 1015 g, where c is the speed of light and G is the Newton constant. For two typical formation times: the Planck time (t ∼ 10−43 s) and the time just before the big-bang nucleosynthesis (t ∼ 1 s), the initial PBH masses are approximately 10−5 g and 105 M⊙, respectively, where M⊙ is the solar mass. Moreover, realistic production mechanisms predict not just a unique mass for all PBHs but rather an extended mass function.
PBHs may constitute all or a fraction of the DM. The fraction of DM in the form of PBHs is defined as fPBH ≡ ΩPBH/ΩDM, where ΩPBH and ΩDM are the energy density parameters of PBHs and DM relative to the critical density of the present Universe, respectively. There exist numerous observational constraints on fPBH (for recent reviews, see, e.g. [23–26]). For heavy PBHs with masses MPBH ≫ 1017 g, constraints on fPBH arise from gravitational effects of PBHs, including lensing [27–33], gravitational waves [34–41], dynamical constraints from globular clusters, Galaxy disruption, and other observables [42–46]. Light PBHs are expected to emit Standard Model particles with a quasi-thermal energy spectrum through Hawking radiation [47, 48], where the temperature of the Hawking spectrum is inversely proportional to the PBH mass. Due to Hawking radiation, PBHs lose their mass at a rate ${\rm{d}}{M}_{{\rm{PBH}}}/{\rm{d}}t\propto {M}_{{\rm{PBH}}}^{-2}$ [49, 50], implying that lighter PBHs evaporate more quickly. It has been shown that PBHs with masses MPBH < 5 × 1014 g have lifetimes shorter than the age of the Universe [51], thus cannot contribute to present DM. In this work, we are particularly interested in the mass range above the evaporation limit (5 × 1014 g) and below the lowest lensing limit (5 × 1017 g). In this mass range, PBHs are expected to emit particles with typical energies from ${ \mathcal O }(10\,{\rm{MeV}})$ down to ${ \mathcal O }(10\,{\rm{keV}})$. A number of constraints on fPBH have been set by considering Hawking radiation in various astrophysical observations, such as extragalactic and galactic γ-rays [52–56], the 511 keV line from the Galactic Center [57–60], cosmic microwave background (CMB) [26, 61, 62], 21 cm radio signals [63, 64], Lyman-α forest measurements [65, 66], neutrinos [57, 67–69], cosmic-ray (CR) electrons [70], CR positrons [71] and CR antiprotons [72, 73].
Recently, it has been noted that strong constraints on fPBH can be derived from direct and indirect CR electron and positron observables when adopting diffusive re-acceleration models of CR propagation [60, 71]. Such models, characterized by a significant Alfvén velocity ${V}_{a}\sim { \mathcal O }(10)\,{\rm{km}}\,{{\rm{s}}}^{-1}$, have received strong support from a number of independent analyses [74–84]. In these models, the energies of the evaporated all-electrons (electrons plus positrons), which typically have initial energies of ${ \mathcal O }(10\,{\rm{MeV}})$, can be enhanced by roughly two orders of magnitude through scattering with the Galactic random magnetic fields during propagation. This effect makes it possible to constrain fPBH using direct and indirect CR electron and positron observables at energies around the GeV scale. For instance, our previous work [71] showed that, in well-constrained diffusive re-acceleration models, a significant portion of the flux of evaporated positrons could be constrained by current AMS-02 data, deriving constraints approximately an order of magnitude stronger than those derived from Voyager-1 all-electron data for MPBH ∼ 1016 g. Similarly, in [60], it was shown that the diffuse x-ray emissions from the up-scattering of Galactic ambient photons due to the inverse Compton effect of the evaporated all-electrons could be constrained by the XMM-Newton data, and the obtained constraints could be improved by two orders of magnitude for MPBH ∼ 1016 g as Va increasing from 13.4 to 40 km s−1 in their benchmark diffusive re-acceleration models. Motivated by these works, we investigate another possible indirect detection channel of PBHs—Galactic diffuse synchrotron emission. The Galactic diffuse synchrotron emission arises from CR all-electrons during their propagation in the Galactic magnetic field (GMF), and constitutes an indirect observable of the interstellar CR all-electrons. The observation of Galactic synchrotron emissions at frequencies above ∼20 MHz can provide indirect measurements of interstellar CR all-electrons with energies above ∼100 MeV, and therefore can set meaningful constraints on fPBH within diffusive re-acceleration models.
To discuss the impact of diffusive re-acceleration on the constraints derived from Galactic synchrotron observations, we fit parameters of a set of benchmark CR propagation models to the AMS-02 [85] and Voyager-1 [86] data on the boron-to-carbon nuclei (B/C) flux ratio. Our model set includes four diffusive re-acceleration models with different diffusion halo half-heights zh (fixed to 4, 6, 8, 10 kpc, respectively), and one diffusion break model with zh fixed to 4 kpc, which does not include the diffusive re-acceleration (i.e. with Va fixed to zero). For diffusive re-acceleration models, our fits obtain significant Alfvén velocities of Va ∼ 20 km s−1. We also employ an additional diffusive re-acceleration model with the parameters fitted in [79], in which solar modulation is modeled using the Helmod code [87–92], to complement our simple treatment of force-field approximation. We show that, for Va ∼ 20 km s−1, a significant fraction of the evaporated all-electrons can be boosted to the energies of ∼100 MeV during their Galactic propagation. Consequently, synchrotron emissions generated by the evaporated all-electrons, hereafter referred to as synchrotron signals of PBHs, can be constrained by the low-frequency radio continuum surveys (from 22 MHz to 1.4 GHz). With diffusive re-acceleration models, we obtain stringent constraints on fPBH. The most conservative constraints are stronger than those derived from the Voyager-1 all-electron data [70] by more than one order of magnitude for MPBH ≳ 1 × 1016 g, and also stronger than our previous constraints derived from the AMS-02 positron data [71] for MPBH ≳ 2 × 1016 g.
The remainder of this paper is organized as follows. In section 2 , we provide a brief overview of the all-electron energy spectrum from PBH evaporation. In section 3 , we describe the CR propagation in the Galaxy and the models for solar modulation, and present the parameter fit results for our benchmark CR propagation models. In section 4 , we describe the Galactic synchrotron emission and the adopted GMF models. In section 5 , we discuss the constraints on fPBH derived from Galactic synchrotron emission observations. Finally, we summarize our work in section 6 .
2. Evaporation of primordial black holes
In this section, we adopt the natural system of units with ℏ = kB = c = 1, where ℏ is the reduced Planck constant, kB is the Boltzmann constant, and c is the speed of light. We consider a simple scenario that the spin of PBHs is negligible, which can be obtained from a series of possible formation mechanisms [93–95]. The emission rate of particle species i per unit total energy E from a PBH of mass MPBH is given by [48]1 ) describes the probability that the particle i created at the PBH horizon finally escapes to spatial infinity, and is determined by the equation of motion of the particle in curved spacetime near the horizon. In the geometric optics limit (i.e. the high-energy limit), the graybody factor for electrons can be approximated as ${{\rm{\Gamma }}}_{e}\simeq 27{G}^{2}{M}_{{\rm{PBH}}}^{2}{E}^{2}$. Note that equation (1 ) only describes the primary particles directly emitted from the PBH. The production of secondary particles from decays of unstable primary particles should also be considered. We use the numerical code BlackHawk [97, 98] to compute the energy spectra of evaporated all-electrons, in which both the primary and secondary production processes are calculated. Following BlackHawk 's manual, we calculate the secondary production process using the results of Hazma code [99], since the energies of primary particles cannot reach 5 GeV for PBH masses we interested in.
$\begin{eqnarray}\frac{{{\rm{d}}}^{2}{N}_{i}}{{\rm{d}}t{\rm{d}}E}=\frac{{g}_{i}{{\rm{\Gamma }}}_{i}}{2\pi }{\left[\exp \left(\frac{E}{{T}_{{\rm{PBH}}}}\right)-{(-1)}^{2{s}_{i}}\right]}^{-1},\end{eqnarray}$
where si and gi are the spin and the total degree of freedom of the particle i, respectively, Γi is the energy-dependent graybody factor of the particle i, and TPBH is the temperature of the PBH which is given by [96] $\begin{eqnarray}{T}_{{\rm{PBH}}}\approx 10.6\times \left(\frac{1{0}^{15}\,{\rm{g}}}{{M}_{{\rm{PBH}}}}\right)\,{\rm{MeV}}.\end{eqnarray}$
The graybody factor Γi in equation (The mass function of PBHs, dn/dMPBH (number density per unit mass), depends on the formation mechanisms of PBHs. In this work, we consider two widely used mass functions: monochromatic [13] and log-normal [12] mass functions. A nearly monochromatic mass function is expected if all PBHs are formed at a same epoch, and a log-normal mass function can arise from inflationary fluctuations [100, 101]. These two mass functions are given by
$\begin{eqnarray}\frac{{\rm{d}}n}{{\rm{d}}{M}_{{\rm{PBH}}}}=\left\{\begin{array}{rlr} & {A}_{1}\delta ({M}_{{\rm{PBH}}}-{M}_{{\rm{c}}}) & \,\rm{(monochromatic)}\,\\ & \frac{{A}_{2}}{\sqrt{2\pi }\sigma {M}_{{\rm{PBH}}}^{2}}{\rm{\exp }}\left[-\frac{{{\rm{ln}}}^{2}({M}_{{\rm{PBH}}}/{M}_{{\rm{c}}})}{2{\sigma }^{2}}\right] & \,\rm{(log-normal)}\,\end{array}\right.\end{eqnarray}$
where Mc is the characteristic mass, σ is the width of log-normal mass distribution, and A1 and A2 are normalization factors determined by the energy density of PBHs, ρPBH, as follows: $\begin{eqnarray}{\int }_{{M}_{{\rm{\min }}}}^{\infty }{M}_{{\rm{PBH}}}\frac{{\rm{d}}n}{{\rm{d}}{M}_{{\rm{PBH}}}}{\rm{d}}{M}_{{\rm{PBH}}}={\rho }_{{\rm{PBH}}},\end{eqnarray}$
where ${M}_{{\rm{\min }}}=5\times 1{0}^{14}\,{\rm{g}}$ corresponds to the mass of which PBHs have completely evaporated today [51], thus the PBHs lighter than ${M}_{{\rm{\min }}}$ cannot contribute to present DM. For PBHs heavier than ${M}_{{\rm{\min }}}$, the evaporation timescales drastically exceed the timescales of CR propagation in the Galaxy, so we neglect the variation of MPBH during evaporation.The energy spectrum of all-electrons evaporated by PBHs with a mass function dn/dMPBH is given by BlackHawk .
$\begin{eqnarray}\frac{{{\rm{d}}}^{2}{n}_{{e}^{\pm }}}{{\rm{d}}t{\rm{d}}E}={\int }_{{M}_{{\rm{\min }}}}^{\infty }\frac{{\rm{d}}n}{{\rm{d}}{M}_{{\rm{PBH}}}}\frac{{{\rm{d}}}^{2}{N}_{{e}^{\pm }}}{{\rm{d}}t{\rm{d}}E}\,{\rm{d}}{M}_{{\rm{PBH}}},\end{eqnarray}$
where ${{\rm{d}}}^{2}{N}_{{e}^{\pm }}/{\rm{d}}t{\rm{d}}E$ is the emission rate of all-electrons from a single PBH of mass MPBH with both the primary and secondary emission components considered. For the monochromatic mass function, above MPBH integral simplifies easily. For the log-normal mass function, the integral is calculated by 3. Cosmic ray propagation
3.1. Cosmic ray propagation in the galaxy
Once injected into the Galaxy, the all-electrons evaporated by PBHs propagate as CRs. When discussing CR propagation, we assume cylindrical symmetry of the Galaxy, and use cylindrical coordinates (R, z). We assume that CR propagation occurs within a cylindrical diffusion halo of radius Rh = 20 kpc and half-height zh (typically a few kpc). The number density of CR particles per unit momentum p at position $\overrightarrow{r}$ at time t, $\psi (\overrightarrow{r},p,t)$, is related to the phase space distribution function $f(\overrightarrow{r},\overrightarrow{p},t)$ as $\psi (\overrightarrow{r},p,t)=4\pi {p}^{2}f(\overrightarrow{r},\overrightarrow{p},t)$, assuming an isotropic momentum distribution of CR particles. A free-escape boundary condition (ψ = 0) is imposed at the boundary of diffusion halo. The Galactic propagation of CRs is described by the diffusion equation: Galprop code [102–106] to solve the diffusion equation numerically. Following the equation (6 ), Galprop evolves ψ long enough time to make ∂ψ/∂t = 0 satisfied, and the result is the steady-state solution of the diffusion equation.
$\begin{eqnarray}\begin{array}{rcl}\frac{\partial \psi }{\partial t} & = & q(\overrightarrow{r},p)+{\rm{\nabla }}\cdot \left({D}_{xx}{\rm{\nabla }}\psi -{\overrightarrow{V}}_{c}\psi \right)\\ & & +\frac{\partial }{\partial p}{p}^{2}{D}_{pp}\frac{\partial }{\partial p}\frac{1}{{p}^{2}}\psi \\ & & -\frac{\partial }{\partial p}\left[\dot{p}\psi -\frac{p}{3}({\rm{\nabla }}\cdot {\overrightarrow{V}}_{c})\psi \right]-\frac{1}{{\tau }_{f}}\psi -\frac{1}{{\tau }_{r}}\psi ,\end{array}\end{eqnarray}$
where q is the source term, which represents the number of CR particles injected into the Galaxy per unit volume, per unit momentum, and per unit time, Dxx is the spatial diffusion coefficient, ${\overrightarrow{V}}_{c}={\rm{sign}}(z)\left({V}_{0}+\frac{{\rm{d}}V}{{\rm{d}}z}| z| \right)\hat{z}$ is the convection velocity driven by the galactic wind with V0 generally assumed zero, Dpp is the diffusion coefficient in momentum space, describing the effect of diffusive re-acceleration, $\dot{p}\equiv {\rm{d}}p/{\rm{d}}t$ is the momentum loss rate, which accounts for physical processes through which CR particles lose energy during propagation, including ionization, Coulomb interactions, bremsstrahlung, inverse Compton and synchrotron processes, τf and τr are the timescales of particle fragmentation and radioactive decay, which account for destruction of CR particles by interaction with interstellar gas and decay of unstable CR particles, respectively. We employ the state-of-the-art In the diffusion equation (6 ), the spatial diffusion coefficient Dxx is parameterized as follows:
$\begin{eqnarray}{D}_{xx}={D}_{0}{\beta }^{\eta }{\left(\frac{\rho }{{\rho }_{0}}\right)}^{\delta },\end{eqnarray}$
where ρ is the rigidity of CR particles, β = v/c is the velocity of CR particles relative to the speed of light c, D0 is the normalization of diffusion coefficient at reference rigidity ρ0, and δ is the spectral power index, for Kolmogorov type of turbulence δ = 1/3 [107] and for an Iroshnikov–Kraichnan cascade δ = 1/2 [108, 109]. The exponent η is introduced to accommodate the low-energy behavior of measured CR spectra as in recent works [79, 110–113], while traditionally set to η = 1. An η < 1 increases the diffusion coefficient at low energies. A negative η phenomenologically accounts for physical nonlinear phenomena at low energies, such as the turbulent dissipation and wave damping that produce a very sharp rise of diffusion coefficient at the rigidity less than around 1.5 GV [114, 115]. Additionally, two breaks of δ could be introduced in certain types of models (see below), in which δ changes to δl and δh when ρ < ρl and ρ > ρh for each break respectively, with ρl ∼ 4 GV and ρh ∼ 300 GV typically.The momentum-space diffusion coefficient Dpp, describing the effect of diffusive re-acceleration, is parametrized as follows:
$\begin{eqnarray}{D}_{pp}=\frac{4{V}_{a}^{2}{p}^{2}}{3{D}_{xx}\delta (4-{\delta }^{2})(4-\delta )},\end{eqnarray}$
where Va is the Alfvén velocity, which characterizes the propagation of disturbances in GMFs. The scattering of charged particles by the random motion of the magnetic fields characterized by the Alfvén velocity leads to a certain amount of second-order Fermi acceleration during propagation, which can significantly modify the low-energy CR spectra. We consider Va to be constant throughout the diffusion halo, which should be understood as an effective parameter. In some semi-analytic frameworks, the re-acceleration is assumed to be confined in the Galactic disc and normalized to a half-width of h ∼ 0.1 kpc, which allows for fast calculations for CR propagation [116]. The value of Va obtained in this semi-analytical approach should be roughly rescaled by a factor of $\sqrt{h/{z}_{{\rm{h}}}}$ when compared with the one adopted in this work.The CR source term q in equation (6 ) includes both the primary and secondary components. The primary CRs are believed to be accelerated by supernova remnants (SNRs) and pulsar wind nebulae, whose source term is modeled as the product of a broken power-law rigidity spectrum and a spatial distribution of the primary source $n(\overrightarrow{r})$. Assuming the spectrum has m breaks at ρi (i = 0, 1, 2, …, m − 1), spectral indices γi and γi+1 below and above ρi, the primary source term is written as 6 ) in ψ (excluding the source term q), the normalization of primary source term is absorbed into the proton flux normalization or primary isotopic abundances, which are free parameters of the model.
$\begin{eqnarray}q(\overrightarrow{r},\,p)=n(\overrightarrow{r})\,{\left(\frac{\rho }{{\rho }_{0}}\right)}^{-{\gamma }_{0}}\,\displaystyle \prod _{i=0}^{m-1}{\left[\frac{{\rm{\max }}(\rho ,{\rho }_{i})}{{\rho }_{i}}\right]}^{{\gamma }_{i}-{\gamma }_{i+1}},\end{eqnarray}$
where the spatial distribution $n(\overrightarrow{r})=n(R,z)$ follows the distribution of SNRs or pulsars [117], given by $\begin{eqnarray}\begin{array}{rcl}n(R,z) & \propto & {\left(\frac{R}{{R}_{\odot }}\right)}^{a}\exp \left(-b\frac{R}{{R}_{\odot }}\right)\\ & & \times \exp \left(-\frac{| z| }{{z}_{0}}\right),\end{array}\end{eqnarray}$
where a = 1.9, b = 5.0 [118], z0 = 0.2 kpc, and R⊙ = 8.5 kpc. We enforce n(R, z) remained constant from R = 10 kpc and fixed to 0 for R > 15 kpc. Due to linearity of equation (Secondary CRs arise from the interaction between primary CRs and interstellar gas during propagation. For a secondary species s, the source term is given by Galprop calculations, j only runs over H I, H2, and H II, while the contributions from interactions with He are accounted for by rescaling the corresponding H interaction cross sections with effective factors.
$\begin{eqnarray}\begin{array}{rcl}{q}_{s}(\overrightarrow{r},p) & = & \displaystyle \sum _{j}{n}_{j}(\overrightarrow{r})\\ & & \displaystyle \times \displaystyle \sum _{i}\int {\rm{d}}{p}_{i}\,c{\beta }_{i}{\psi }_{i}(\overrightarrow{r},{p}_{i})\frac{{\rm{d}}{\sigma }_{ij\to s}}{{\rm{d}}p}(p,{p}_{i}),\end{array}\end{eqnarray}$
where i indexes primary CR species producing the species s, j indexes interstellar gas components, nj is number density of the gas, and dσij→s/dp is the differential production cross section. The interstellar gas consists mostly of hydrogen (H) and helium (He) with a number ratio of 9:1. H has different states: atomic (H I), molecular (H2) and ionized (H II), and He is mostly neutral. In practical PBHs in the Galaxy can be regarded as stable sources of CR all-electrons due to Hawking radiation. Assuming PBHs constitute a fraction fPBH of the DM, their density profile follows that of the DM as ${\rho }_{{\rm{PBH}}}(\overrightarrow{r})={f}_{{\rm{PBH}}}{\rho }_{{\rm{DM}}}(\overrightarrow{r})$. We adopt NFW DM profile [119] with scale radius rs = 24.42 kpc, and Burkert DM profile [120] with scale radius rs = 12.67 kpc. Both DM profiles are normalized such that ρDM(r⊙) = 0.4 GeV cm−3 at the Solar position [121]. The NFW profile is cuspy towards the Galactic center, while the Burkert profile exhibits a constant density core, and predicts a smaller DM density in the innermost part of the Galaxy. The source term for the evaporated all-electrons, with a mass function dn/dMPBH of PBHs, is given by 2 ), and dn/dMPBH is normalized to ${\rho }_{{\rm{PBH}}}(\overrightarrow{r})$ = ${f}_{{\rm{PBH}}}\,{\rho }_{{\rm{DM}}}(\overrightarrow{r})$ through equation (4 ), thus depends on $\overrightarrow{r}$.
$\begin{eqnarray}\begin{array}{rcl}{q}_{{e}^{\pm }}(\overrightarrow{r},p) & = & {\displaystyle \int }_{{M}_{{\rm{\min }}}}^{\infty }\frac{{\rm{d}}n}{{\rm{d}}{M}_{{\rm{PBH}}}}\\ & & \times \beta \frac{{{\rm{d}}}^{2}{N}_{{e}^{\pm }}}{{\rm{d}}t{\rm{d}}E}\,{\rm{d}}{M}_{{\rm{PBH}}},\end{array}\end{eqnarray}$
where ${M}_{{\rm{\min }}}=5\times 1{0}^{14}\,{\rm{g}}$, ${{\rm{d}}}^{2}{N}_{{e}^{\pm }}/{\rm{d}}t{\rm{d}}E$ is the emission rate of a single PBH of mass MPBH (see section 3.2. Cosmic ray propagation in the heliosphere
As CRs propagate from the heliopause (HP), i.e. heliosphere boundary, to the Earth, they undergo affects by heliospheric magnetic fields generated by the solar wind. These physical processes suppress CR spectra below ∼50 GV measured at the Earth compared to the interstellar spectra and introduce time-dependent variations related to solar activities [122]. This effect is known as solar modulation. Since all CR measurements—except Voyager-1 (since August 2012) and Voyager-2 (since December 2018)—are conducted within the heliosphere, CR propagation studies require the modeling of both galactic propagation and solar modulation. The propagation of CRs in the heliosphere is described by the Parker equation [123]: Helmod code [87–92] offers modulation tables based on numerical solutions of equation (13 ) with parameters calibrated using the up-to-date CR data. These tables relate local interstellar (LIS) and top-of-atmosphere (TOA) CR spectra, denoted by ΦLIS and ΦTOA, via Helmod code is able to quantitatively reproduce the time variation of CR fluxes, such as that of protons, and the predicted LIS proton flux is in remarkable agreement with the Voyager-1 data [91]. Other numerical codes solving the Parker equation include SOLARPROP [124] and HelioProp [125], etc.
$\begin{eqnarray}\begin{array}{rcl}\frac{\partial U}{\partial t} & = & {\rm{\nabla }}\cdot \left({K}^{S}{\rm{\nabla }}U-({\overrightarrow{V}}_{{\rm{sw}}}+{\overrightarrow{v}}_{{\rm{d}}})U\right)\\ & & +\frac{1}{3}{\rm{\nabla }}\cdot {\overrightarrow{V}}_{{\rm{sw}}}\frac{\partial }{\partial T}\left({\alpha }_{{\rm{rel}}}TU\right),\end{array}\end{eqnarray}$
where $U(\overrightarrow{r},t,T)$ is the number density of CR particles per unit kinetic energy T, KS is the symmetric part of diffusion tensor, ${\overrightarrow{V}}_{{\rm{sw}}}$ is the solar wind velocity, ${\overrightarrow{v}}_{{\rm{d}}}$ is the magnetic drift velocity, and αrel = (T + 2T0)/(T + T0) with T0 denoting the rest energy of the CR particle. The public $\begin{eqnarray}{{\rm{\Phi }}}^{{\rm{TOA}}}(\rho )=\int {{\rm{\Phi }}}^{{\rm{LIS}}}({\rho }^{{\prime} })\,G(\rho ,{\rho }^{{\prime} };t)\,{\rm{d}}{\rho }^{{\prime} },\end{eqnarray}$
where $G(\rho ,{\rho }^{{\prime} };t)$ describes the probability for a CR particle with initial rigidity ${\rho }^{{\prime} }$ at the HP to be observed at TOA with rigidity ρ, and t denotes the time interval of the CR measurement. The Numerically solving the Parker equation is computationally demanding for CR propagation studies. Consequently, the much simpler force-field approximation [126] is widely adopted. This model provides an analytical solution of the Parker equation under a series of assumptions: steady-state system, spherical symmetry, and zero streaming. Assuming the diffusion coefficient is separable k(r, ρ) = βk1(r)k2(ρ), where β = v/c and r is the distance from the Sun, the solution relates LIS and TOA CR spectra via
$\begin{eqnarray}\begin{array}{rcl}\frac{{{\rm{\Phi }}}^{{\rm{TOA}}}}{{({\rho }^{{\rm{TOA}}})}^{2}} & = & \frac{{{\rm{\Phi }}}^{{\rm{LIS}}}}{{({\rho }^{{\rm{LIS}}})}^{2}}\\ & \Rightarrow & {{\rm{\Phi }}}^{{\rm{TOA}}}({T}^{{\rm{TOA}}})\\ & = & \frac{{T}^{{\rm{TOA}}}({T}^{{\rm{TOA}}}+2M)}{{T}^{{\rm{LIS}}}({T}^{{\rm{LIS}}}+2M)}{{\rm{\Phi }}}^{{\rm{LIS}}}({T}^{{\rm{LIS}}}),\end{array}\end{eqnarray}$
where ρ is the rigidity, T is the kinetic energy, and M is the particle mass. The rigidities ρLIS and ρTOA are connected through $\begin{eqnarray}{\int }_{{\rho }^{{\rm{TOA}}}}^{{\rho }^{{\rm{LIS}}}}\frac{\beta {k}_{2}(\rho )}{\rho }{\rm{d}}\rho ={\int }_{{r}_{{\rm{TOA}}}}^{{r}_{{\rm{HP}}}}\frac{{V}_{{\rm{sw}}}(r)}{3{k}_{1}(r)}{\rm{d}}r\equiv \phi ,\end{eqnarray}$
where φ is the modulation potential. For k2(ρ) ∝ ρ, this gives $\begin{eqnarray}{T}^{{\rm{LIS}}}={T}^{{\rm{TOA}}}+| Ze| \phi ,\end{eqnarray}$
where ∣Ze∣ is the absolute charge of the CR particle. In practice, φ is considered as a free nuisance parameter in the fitting procedure of CR propagation models. There exists a strong degeneracy between φ and other model parameters. The validity of force-field approximation is challenged by variation of CR monthly fluxes measured by PAMELA and AMS-02 [127]. There are proposed modifications to the force-field approximation by introducing rigidity-dependent [128–130] or charge-asymmetric [131, 132] modulation potential. However, for CR propagation studies concerning CR data measured over a long time period, the force-field approximation is widely adopted.3.3. Benchmark propagation models
In sections 3.1 and 3.2 , we have introduced the basic physical concepts of CR propagation in the Galaxy and heliosphere. However, the specific parameter values, including those in the diffusion equation (6 ), primary CR source term equation (9 ), and those for solar modulation (in the parker equation (13 ) or the modulation potential φ), need to be determined through fitting to CR measurement data [106]. Secondary-to-primary nuclei flux ratios, such as the most widely used B/C flux ratio, can probe the grammage (representing the integrated gas density along the CR trajectory before escape from the diffusion halo) of CRs and are primarily sensitive to propagation parameters, as the source dependence largely cancel out. It is known that re-acceleration can provide a natural mechanism to reproduce the low-energy B/C flux ratio with Kolmogorov type of turbulence and a spectral break in the power-law spectra of primary CRs at a rigidity of a few GV [133]. Diffusive re-acceleration models with a significant Alfvén velocity ${V}_{a}\sim { \mathcal O }(10)\,{\rm{km}}\,{{\rm{s}}}^{-1}$ have received strong support from a number of independent analyses [74–81], in which the values of Va are within 30–45 km s−1. Several recent analyses obtained slightly smaller Va ∼ 20 km s−1 [81–84], when including additional model free parameters, such as introducing a low-energy δ break or a free parameter η, and considering the high-energy break of δ around 300 GV. Notably, analyses fitting to proton, antiproton, and He data preferred a substantially lower Va [76, 134], suggesting different CR species may sample distinct interstellar environments due to spatially dependent propagation effects.
Since the effects of propagation parameters are partially degenerated, it is possible to construct alternative CR propagation models without the re-acceleration (i.e. with Va fixed to zero) but including a low-energy break of δ at ρl ∼ 4 GV. These types of diffusion break models can reproduce the similar structure in B/C flux ratio. So far, both scenarios can well explain the data of secondary-to-primary nuclei flux ratio [83, 86, 113, 135]. Notably, it is possible to distinguish the two scenarios by considering additional observables. For example, the existence of significant re-acceleration can be probed by synchrotron and soft γ-ray emissions [136], and a break in the injection primary source can be examined by γ-ray emissions from gas clouds [137].
For constraining fPBH, the diffusion halo half-height zh is a critical parameter. Since CRs escape freely at the boundary of diffusion halo, only the PBHs within the diffusion halo can contribute to observable signals. Thus, the uncertainty of zh can directly affect the constraints on fPBH. Because of the famous degeneracy between zh and D0, it is hard to constrain the value of zh through fitting to the data of secondary-to-primary nuclei flux ratios. The flux ratio of unstable-to-stable secondary isotopes can provide insights into the residence time of CRs within the diffusion halo. The life time of 10Be is about 2 million years, comparable to the typical CR residence time. Hence the 10Be/9Be flux ratio is an ideal observable for constraining zh. However, due to the large uncertainty of Be isotopic data (such as the data of ACE-CRIS, ACE-SIS, ISOMAX, and PAMELA) and their relatively low energy upper limit ≲2 GeV/n, it is still hard to obtain stringent constraints on zh. As a substitution, the data of Be/B flux ratio were used to constrain zh in several recent works. For example, in [138] a lower bound zh ≥ 5 kpc was set, and in [139]${z}_{{\rm{h}}}={5}_{-2}^{+3}\,{\rm{kpc}}$ was found. A relatively smaller ${z}_{{\rm{h}}}=3.{8}_{-1.6}^{+2.8}\,{\rm{kpc}}$ was found in the AMS-02 Be/B analysis in [140], and zh = 4.7 ± 2 kpc was obtained using the 10Be/9Be data in the same work. In [79], zh = 4.0 ± 0.6 kpc was obtained using the CR data of light nuclei (with charge number Z ≤ 28), but the uncertainties of production cross sections of secondary CRs were not considered. In 2022 ICHEP, AMS-02 collaboration preliminarily reported the Be isotopic measurements,5 which have reached an unprecedented energy of 12 GeV/n. However, the uncertainties of the production cross sections of secondary CRs still prevent the community to obtain a robust constraint on zh. Using the preliminary data of 7Be and 10Be nuclei, zh = 5.67 ±0.76 kpc was obtained in [141], where the cross section uncertainties were considered by adopting a cross section parametrization that fully utilizes the available experimental data. Besides the uncertainties of production cross sections, it is also reported that the Galactic gas model [142] and inhomogeneous diffusion [143] can affect the constraints on zh derived from the 10Be/9Be data. Given the huge difficulties in constraining zh, instead of choosing one typical zh from the earlier works, we fit parameters of 4 diffusive re-acceleration models with zh fixed to 4, 6, 8, 10 kpc respectively, covering zh values reported in the literature. With these models, we can illustrate the affect of zh on the obtained constraints on fPBH.
To discuss the impact of diffusive re-acceleration on the constraints derived from Galactic synchrotron observations, we fit parameters of a set of benchmark CR propagation models to the AMS-02 [85] and Voyager-1 [86] data on the B/C flux ratio. Besides, in calculations of constraints on fPBH, we also employ an additional diffusive re-acceleration model with the parameters fitted in [79], in which solar modulation is modeled using the Helmod code [87–92], to complement our simple treatment of force-field approximation. Our benchmark models are listed as follows:
DRz4-10 models. Four diffusive re-acceleration models with different zh (fixed to 4, 6, 8, 10 kpc, respectively). As discussed above, the values of zh are selected to illustrate the affect of zh on the obtained constraints on fPBH. We have four free parameters for CR propagation, which are D0 (normalized at 4 GV), η, δ, and Va.
DBz4 model. A diffusion break model with zh fixed to 4 kpc. We have four free parameters for CR propagation, which are D0 (normalized at 4 GV), ρl, δl and δ. In this model, Va is fixed to zero.
GH model. The model was built with analysis framework of
Other parameters of the DRz4-10 and DBz4 models are described below. The high-energy spectral hardening around 300 GV has been discovered by several CR measurements [144–147]. We model this feature by including the high-energy δ break in the diffusion coefficient Dxx [148, 149], which is described by 18 ) for a clear expression. For primary source parameters, we assume one spectral break in the primary CR source term equation (9 ), having four free parameters: the break rigidity ρ0, the spectral indices γ0 and γ1 below and above ρ0, and the primary abundance of 12C isotope AC12 (with the proton abundance fixed at 1.06 × 106). For solar modulation, we use the force-field approximation with a single free parameter: modulation potential φ for the data measured by AMS-02. Altogether we have nine free parameters for both types of models.
$\begin{eqnarray}{D}_{xx}=\left\{\begin{array}{rlr} & {D}_{0}{\beta }^{\eta }{(\rho /{\rho }_{0})}^{\delta } & \,\rm{for}\,\rho \lt {\rho }_{h},\\ & {D}_{0}{\beta }^{\eta }{({\rho }_{h}/{\rho }_{0})}^{\delta }{(\rho /{\rho }_{h})}^{{\delta }_{h}} & \,\rm{for}\,\rho \gt {\rho }_{h},\end{array}\right.\end{eqnarray}$
where ρh = 308 GV and δh = δ − 0.175 are fixed to the best-fit parameters obtained by [84]. Because we mostly focus on energies below 10 GeV, and to reduce parameter space dimensionality, we adopt these values from [84] rather than fit them independently. The low-energy δ break for the diffusion break model is omitted in equation (Below, we describe our model fitting procedure, display the estimated model parameters, and show the B/C flux radio predictions for the benchmark models. The numerical details and other fit results are put in appendix A . In the last part of this section, we show energy spectra of evaporated all-electrons for CR propagation models adopted in this work.
We adopt Bayesian inference framework to estimate parameter values, with uniform prior probability distribution functions (PDFs) over prior ranges, and a likelihood function in Gaussian form: MultiNest package [150–152].
$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L }({\rm{\Theta }}) & = & \displaystyle \prod _{d,i}\frac{1}{\sqrt{2\pi {\sigma }_{d,i}^{2}}}\\ & & \times \exp \left(-\frac{{\left({{\rm{\Phi }}}_{d}({E}_{i},{\rm{\Theta }})-{{\rm{\Phi }}}_{d,i}\right)}^{2}}{2{\sigma }_{d,i}^{2}}\right),\end{array}\end{eqnarray}$
where Θ denotes model parameters, d indexes CR measurement datasets, i indexes energy bins, σd,i is the experimental uncertainty, Φd(Ei, Θ) is the model-predicted flux for dataset d at ith energy bin with parameters Θ, and Φd,i is the corresponding measured value. Assuming uncorrelated errors, we compute the experimental uncertainty σd,i by adding the reported statistical and systematic errors in quadrature. The posterior PDF, which is the normalized product of prior PDFs and the likelihood function, is sampled with the In table 1, we list the model parameters, and their prior ranges, posterior means and standard deviations of the benchmark models. In figure 1, we show model predictions and residuals of B/C flux ratio for the benchmark models, and the results for C flux are put in appendix A . It can be seen that all DRz4-10 models provide equally good fits to Voyager-1 and AMS-02 data, yielding nearly identical predictions. While the DBz4 model exhibits slightly insufficient spectral hardening around 300 GV, because we fix the ρh and δ − δh same as in diffusive re-acceleration models. Since the value of δ in diffusion break models is generally larger than that in diffusive re-acceleration models, the fixed parameters for the high-energy δ break certainly induce some biases in our diffusion break model. A small defect of our models is that they all underpredict the B/C flux ratio data point of Voyager-1 at the lowest energy, a feature also noted in [84, 86]. Nevertheless, our fits show that both types of models can well reproduce the B/C flux ratio data. Notably, our fits obtain significant Alfvén velocities Va ∼ 20 km s−1 for all diffusive re-acceleration models.
Table 1. Model parameters, and their prior ranges, posterior means and standard deviations (show in brackets) of the benchmark models. |
| Parameter | Prior range | DBz4 | DRz4 | DRz6 | DRz8 | DRz10 |
|---|---|---|---|---|---|---|
| D0 (1028 cm2 s−1) | [1, 10] | 4.23 (0.18) | 4.20 (0.17) | 5.86 (0.23) | 7.04 (0.28) | 7.84 (0.31) |
| ρl (GV) | [0.5, 10] | 5.12 (0.36) | — | — | — | — |
| δl | [−1, 0] | −0.374 (0.090) | — | — | — | — |
| δ | [0.2, 0.6] | 0.490 (0.008) | 0.440 (0.011) | 0.437 (0.011) | 0.436 (0.011) | 0.435 (0.011) |
| η | [−2, 2] | — | −0.628 (0.137) | −0.607 (0.135) | −0.595 (0.134) | −0.583 (0.138) |
| Va (km s−1) | [0, 60] | — | 21.1 (1.7) | 20.7 (1.7) | 20.1 (1.7) | 19.5 (1.6) |
| γ0 | [0.2, 3.2] | 0.885 (0.122) | 0.848 (0.116) | 0.869 (0.111) | 0.879 (0.113) | 0.885 (0.111) |
| γ1 | [1.8, 3.2] | 2.346 (0.008) | 2.369 (0.008) | 2.372 (0.008) | 2.373 (0.008) | 2.374 (0.008) |
| ρ0 (GV) | [0.1, 20] | 1.52 (0.11) | 1.54 (0.10) | 1.54 (0.10) | 1.54 (0.10) | 1.54 (0.09) |
| AC12 | [1000, 6000] | 3487 (28) | 3421 (28) | 3410 (28) | 3407 (28) | 3404 (28) |
| φ (MV) | [0, 1500] | 470 (27) | 606 (27) | 606 (26) | 605 (26) | 605 (26) |
Figure 1. Model predictions and residuals of B/C flux ratio for the benchmark models. Left: for Voyager-1. Right: for AMS-02. |
The DAMPE experiment has reported high-precision measurements of the B/C flux ratio in the energy range from 10 GeV/n to 5.6 TeV/n [147], providing stringent constraints on propagation parameters such as ρh and δh. Because neither our fitting procedure nor that of [84] incorporated the DAMPE B/C data, we have checked the compatibility of our benchmark models with the DAMPE measurements. We find that all models exhibit a mild tendency to underpredict the B/C flux ratio at energies above ∼100 GeV/n. For the DRz4-10 models, the deviation is at the level of approximately one experimental uncertainty, while for the DBz4 model the maximal deviation reaches roughly twice the experimental uncertainty. Overall, the predicted B/C flux ratios remain broadly consistent with the DAMPE data.
In table 2, we summarize the model parameters relevant to the calculation of synchrotron signals of PBHs, including those obtained from our fits and those of the GH model [79]. In figure 2, we show the LIS energy spectra of evaporated all-electrons at the Solar position for CR propagation models adopted in this work, assuming a monochromatic mass function with Mc = 9 × 1015 g, fPBH = 1.0, and NFW DM profile. It can be seen that, for energies above 10 MeV, the all-electron flux of the diffusion break model (DBz4) drops steeply, while those of diffusive re-acceleration models (DRz4-10 and GH) decrease much slower as the energy increasing. We show that, for Va ∼ 20 km s−1, a significant fraction of evaporated all-electrons can be boosted to energies of ∼100 MeV. As zh increasing from 4 kpc to 10 kpc in DRz4-10 models, the all-electron flux increases by about a factor of 3 at energies of ∼100 MeV. Notably, compared with the DRz4 model, the GH model with the same zh = 4 kpc predicts a smaller all-electron flux at energies of ∼100 MeV, despite having a larger Va = 30 km s−1. While for larger energies of ∼1 GeV, the all-electron fluxes of these two models are in agreement. This feature may arise from the different values of δ in these two models, which can affect spectral slope of the LIS all-electron spectra predicted by the models. The very distinct values of η, however, have no effect for these highly relativistic energies of all-electrons. The impact of this feature also shows up in the resulting constraints on fPBH, as will be discussed in section 5.3 .
Table 2. Model parameters relevant to the calculation of synchrotron signals of PBHs for CR propagation models adopted in this work. |
| Parameter | DBz4 | GH [79] | DRz4 | DRz6 | DRz8 | DRz10 |
|---|---|---|---|---|---|---|
| zh (kpc) | 4 | 4 | 4 | 6 | 8 | 10 |
| D0 (1028 cm2 s−1) | 4.23 | 4.3 | 4.20 | 5.86 | 7.04 | 7.84 |
| ρl (GV) | 5.12 | — | — | — | — | — |
| δl | −0.374 | — | — | — | — | — |
| δ | 0.490 | 0.415 | 0.440 | 0.437 | 0.436 | 0.435 |
| η | — | 0.7 | −0.628 | −0.607 | −0.595 | −0.583 |
| Va (km s−1) | — | 30 | 21.1 | 20.7 | 20.1 | 19.5 |
| dV/dz (km s−1 kpc−1) | — | 9.8 | — | — | — | — |
Figure 2. The LIS energy spectra of evaporated all-electrons at the Solar position for CR propagation models adopted in this work, assuming a monochromatic mass function with Mc = 9 × 1015 g, fPBH = 1.0, and NFW DM profile. |
4. Synchrotron emission and GMF models
4.1. Basic physical processes
The Galactic synchrotron emission is generated by CR all-electrons during propagation in the GMF, which is the major component of the Galactic radio emission from ${ \mathcal O }(10\,{\rm{MHz}})$ to ${ \mathcal O }(10\,{\rm{GHz}})$. With a typical GMF strength of ${ \mathcal O }(\mu {\rm{G}})$, the spectral observation of Galactic synchrotron emissions in this frequency range can provide indirect measurements of interstellar CR all-electrons in the energy range from ${ \mathcal O }(100\,{\rm{MeV}})$ to ${ \mathcal O }(10\,{\rm{GeV}})$. We use the Galprop code to calculate synchrotron emissions [153, 154]. Below, we introduce several physical processes related to the calculation of synchrotron emissions and the GMF modeling.
Let ${\overrightarrow{B}}_{\perp }$ denote the projection of magnetic field vector onto the plane perpendicular to the line of sight. The synchrotron emissivity (i.e. power per unit volume per unit frequency) of an isotropic distribution of monoenergetic relativistic all-electrons is partially linearly polarized, which has polarized components parallel and perpendicular to the ${\overrightarrow{B}}_{\perp }$, denoted by j∥ and j⊥, as follows [155]:
$\begin{eqnarray}{j}_{\parallel }(\nu )=\frac{\sqrt{3}}{2}\frac{{q}_{e}^{3}}{{m}_{e}{c}^{2}}{B}_{\perp }[F(x)-G(x)],\end{eqnarray}$
$\begin{eqnarray}{j}_{\perp }(\nu )=\frac{\sqrt{3}}{2}\frac{{q}_{e}^{3}}{{m}_{e}{c}^{2}}{B}_{\perp }[F(x)+G(x)],\end{eqnarray}$
where qe is the absolute charge of electron, me is the electron mass, c is the speed of light, x = ν/νc with νc = 3qeB⊥γ2/(4πmec), where γ is the Lorentz factor. The functions F(x) and G(x) are defined as $\begin{eqnarray}\begin{array}{rc}F(x) & =\,x{\displaystyle \int }_{x}^{\infty }{K}_{5/3}({x}^{{\prime} }){\rm{d}}{x}^{{\prime} },\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rc}G(x) & =\,x{K}_{2/3}(x),\end{array}\end{eqnarray}$
where K5/3 and K2/3 are the modified Bessel functions of orders 5/3 and 2/3, respectively. The total emissivity jtot and polarized emissivity jpol are given by $\begin{eqnarray}{j}_{{\rm{tot}}}(\nu )={j}_{\perp }(\nu )+{j}_{\parallel }(\nu ),\end{eqnarray}$
$\begin{eqnarray}{j}_{{\rm{pol}}}(\nu )={j}_{\perp }(\nu )-{j}_{\parallel }(\nu ).\end{eqnarray}$
We adopt the IAU coordinate convention for Stokes parameters [156]. The intrinsic polarization angle χ0 is measured from the Galactic north direction to the polarization direction (perpendicular to the ${\overrightarrow{B}}_{\perp }$), increasing towards the Galactic east, with χ0 varying in [0, π). The emissivities for Stokes parameters Q and U, denoted by jQ and jU, are given by $\begin{eqnarray}{j}_{Q}={j}_{{\rm{pol}}}\cos (2{\chi }_{0}),\quad {j}_{U}={j}_{{\rm{pol}}}\sin (2{\chi }_{0}).\end{eqnarray}$
Radiation due to the acceleration of an electron/positron moving in the Coulomb field of an ion is known as bremsstrahlung or free–free emission. The inverse process, in which radiation is absorbed by an electron/positron moving in the Coulomb field of an ion, is known as free–free absorption. For interstellar ionized gas, the free–free emission becomes significant at frequencies above a few GHz, and the free–free absorption affects the radiation below about 100 MHz. As modeled in Galprop [154, 157], the free–free opacity kff and emissivity eff at frequency ν are given by
$\begin{eqnarray}\begin{array}{r}{k}_{ff}(\nu ,{n}_{{\rm{TE}}},{T}_{{\rm{TE}}})=0.0178\,{g}_{ff}(\nu ,{T}_{{\rm{TE}}})\frac{{n}_{{\rm{TE}}}^{2}}{{\nu }^{2}{T}_{{\rm{TE}}}^{3/2}},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{r}{e}_{ff}(\nu ,{n}_{{\rm{TE}}},{T}_{{\rm{TE}}})=5.444\times 1{0}^{-39}{g}_{ff}(\nu ,{T}_{{\rm{TE}}})\frac{{n}_{{\rm{TE}}}^{2}}{{T}_{{\rm{TE}}}^{1/2}},\end{array}\end{eqnarray}$
where ${g}_{ff}(\nu ,{T}_{{\rm{TE}}})=10.6+1.9\,{\mathrm{log}}_{10}({T}_{{\rm{TE}}})-1.26\,{\mathrm{log}}_{10}(\nu )$, nTE is the number density of thermal electrons, which equals the number density of ionized hydrogen (H II), and TTE is the temperature of thermal electrons. Since kff and eff both scale with the local ${n}_{{\rm{TE}}}^{2}$, the opacity and emissivity must be multiplied by a clumping factor C to capture the effect of gas inhomogeneities on the squared-density. Following [154], we adopt TTE = 7000 K and C = 100. The optical depth τν along the line of sight is given by $\begin{eqnarray}{\tau }_{\nu }(s)={\int }_{0}^{s}C\,{k}_{ff}(\nu ,{n}_{{\rm{TE}}},{T}_{{\rm{TE}}})\,{\rm{d}}r,\end{eqnarray}$
where s is the line-of-sight distance, and nTE depends on position r along the line of sight.The polarization angle of an electromagnetic wave rotates when passing through a magnetized plasma. This effect is known as Faraday rotation. The rotated angle is proportional to the square of the wavelength λ. For an initial polarization angle χ0, the polarization angle χ after passing through the plasma becomes
$\begin{eqnarray}\chi ={\rm{RM}}\cdot {\lambda }^{2}+{\chi }_{0},\end{eqnarray}$
where the rotation measure (RM) quantifies the linear change rate of χ − χ0 with respect to λ2, and is determined by the plasma property as $\begin{eqnarray}{\rm{RM}}=\frac{{q}_{e}^{3}}{2\pi {m}_{e}^{2}{c}^{4}}{\int }_{{\rm{los}}}{\rm{d}}r\,{n}_{{\rm{TE}}}\,{B}_{\parallel },\end{eqnarray}$
where nTE is the number density of thermal electrons, B∥ is the magnetic field component parallel to the line of sight, and the subscript LOS stands for line-of-sight integration.Let ${n}_{{\rm{CR}}}(\overrightarrow{r},\gamma )$ denote the number density of CR all-electrons with Lorentz factor γ at position $\overrightarrow{r}$. For a detector of area dA oriented towards a small solid angle dΩ, the received radiation power per unit frequency is given by 33 ), given by
$\begin{eqnarray}\begin{array}{rcl}\frac{{\rm{d}}W}{{\rm{d}}\nu } & = & {\displaystyle \int }_{{\rm{los}}}\frac{{\rm{d}}A}{4\pi {r}^{2}}\left(\displaystyle \int {j}_{{\rm{tot}}}(\nu ,\gamma )\,{n}_{{\rm{CR}}}(\overrightarrow{r},\gamma ){\rm{d}}\gamma \right)\\ & & \times {{\rm{e}}}^{-{\tau }_{\nu }(r)}\,{r}^{2}{\rm{d}}{\rm{\Omega }}{\rm{d}}r,\end{array}\end{eqnarray}$
where the subscript LOS stands for line-of-sight integration, and r is the distance to the detector. The measured intensity I is given by $\begin{eqnarray}\begin{array}{rcl}I(\nu ) & = & \frac{{\rm{d}}W}{{\rm{d}}A{\rm{d}}{\rm{\Omega }}{\rm{d}}\nu }=\frac{1}{4\pi }{\displaystyle \int }_{{\rm{los}}}\\ & & \times \left(\displaystyle \int {j}_{{\rm{tot}}}(\nu ,\gamma ){n}_{{\rm{CR}}}(\overrightarrow{r},\gamma ){\rm{d}}\gamma \right){{\rm{e}}}^{-{\tau }_{\nu }(r)}{\rm{d}}r.\end{array}\end{eqnarray}$
The polarized intensity can be calculated analogously by replacing jtot with jpol. For Stokes parameters Q and U, the effect of Faraday rotation must be considered. It is helpful to introduce the complex polarized intensity: $\begin{eqnarray}\tilde{P}\equiv P{{\rm{e}}}^{2{\rm{i}}\chi }=Q+{\rm{i}}U,\end{eqnarray}$
where $P=| \tilde{P}| $ is the polarized intensity, and χ is the polarization angle. The equation follows directly from the fact that the synchrotron emissivity is partially linearly polarized. The measured $\tilde{P}$ is then obtained by replacing jtot with jpole2iχ in equation ( $\begin{eqnarray}\tilde{P}(\nu )=\frac{1}{4\pi }{\int }_{{\rm{los}}}\left(\int {j}_{{\rm{pol}}}(\nu ,\gamma ){{\rm{e}}}^{2{\rm{i}}\chi (r)}{n}_{{\rm{CR}}}(\overrightarrow{r},\gamma ){\rm{d}}\gamma \right){{\rm{e}}}^{-{\tau }_{\nu }(r)}{\rm{d}}r,\end{eqnarray}$
where χ(r) = RM · λ2 + χ0 is the observed polarization angle of the emission from a volume element at line-of-sight distance r, which depends on position due to variations of both RM and χ0 along the line of sight.It is useful to provide results under the assumption of a constant spectral index α of CR all-electrons, implying nCR(γ) = N0γ−α. In this case, the integrals over γ in equations (33 ) and (35 ) can be evaluated analytically:
$\begin{eqnarray}\begin{array}{rc}I(\nu ) & ={\displaystyle \int }_{{\rm{los}}}{J}_{I}\,{{\rm{e}}}^{-{\tau }_{\nu }(r)}\,{\rm{d}}r,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rc}\tilde{P}(\nu ) & ={\displaystyle \int }_{{\rm{los}}}{J}_{P}\,{{\rm{e}}}^{2{\rm{i}}\chi (r)}\,{{\rm{e}}}^{-{\tau }_{\nu }(r)}\,{\rm{d}}r,\end{array}\end{eqnarray}$
with $\begin{eqnarray}\begin{array}{l}{J}_{I}=\frac{\sqrt{3}{q}_{e}^{3}{N}_{0}}{4\pi {m}_{e}{c}^{2}(1+\alpha )}{B}_{\perp }^{\frac{1+\alpha }{2}}{\left(\frac{2\pi \nu {m}_{e}c}{3{q}_{e}}\right)}^{\frac{1-\alpha }{2}}\\ \,\times {\rm{\Gamma }}\left(\frac{\alpha }{4}+\frac{19}{12}\right){\rm{\Gamma }}\left(\frac{\alpha }{4}-\frac{1}{12}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{J}_{P}=\frac{\sqrt{3}{q}_{e}^{3}{N}_{0}}{16\pi {m}_{e}{c}^{2}}{B}_{\perp }^{\frac{1+\alpha }{2}}{\left(\frac{2\pi \nu {m}_{e}c}{3{q}_{e}}\right)}^{\frac{1-\alpha }{2}}\\ \,\times {\rm{\Gamma }}\left(\frac{\alpha }{4}+\frac{7}{12}\right){\rm{\Gamma }}\left(\frac{\alpha }{4}-\frac{1}{12}\right),\end{array}\end{eqnarray}$
where Γ is the Gamma function. For a spectral index α = 3, which is approximately true for the CR all-electron spectrum at energies of ${ \mathcal O }(10\,{\rm{GeV}})$, both JI and JP scale with ${N}_{0}{B}_{\perp }^{2}{\nu }^{-1}$, providing intuitions to the frequency spectrum and the magnetic field dependence of the Galactic synchrotron emissions.4.2. GMF models
To calculate synchrotron signals of PBHs, a realistic model of the GMF is required. However, our understanding of the GMF is limited by both our position within the Galactic disc and the challenges in interpreting observables: RM, synchrotron intensity (I), and polarized synchrotron intensity (P). As described above, these observables are line-of-sight integrated quantities which probe different GMF components: the RM probes the weighted average of GMF component parallel to the line of sight, with the density of thermal electrons as the weight; while the synchrotron emission (I and P) probes the weighted average of the squared perpendicular GMF component, with the density of CR electrons as the weight. The contribution of CR positrons is generally negligible in construction of GMF models, which mainly concerns synchrotron emission observations at 408 MHz and 23 GHz. However, CR positrons can have significant affects on synchrotron emissions at lower frequencies of ≲100 MHz as discussed in [136]. There are no direct measurements of the GMF that do not depend on other components of the interstellar medium. Furthermore, these observables can be contaminated by other physical processes, such as free–free emission and dust emission.
The interstellar medium is known to be turbulent, leading to fluctuations in the GMF on scales larger than ∼100 pc [158, 159]. The modeling of turbulent interstellar medium is another challenge in construction of GMF models [160]. To simultaneously explain the observations of RM, I, and P, in many works, the GMF is modeled with three effective component fields that contribute differently to these observables (for a recent review, see [161]). These are regular field, striated random field,6 and isotropic random field. The regular field maintains coherence on kpc scale and contributes to all observables: RM via equation (31 ), I via equation (33 ), and P via equation (35 ). The striated random field exhibits random strength variations but maintains alignment with a specific direction (typically but not necessarily, the direction of the regular field). Only its direction and variance of the strength are relevant. The striated random field contributes to I and P but not to RM, since its line-of-sight projections cancel out. The isotropic random field varies randomly in both strength and direction, and only the strength variance is relevant. It contributes solely to I due to the isotopic random direction.
The synchrotron emissivity of the random fields is modeled in several ways. In the Galprop [154] framework, the striated random field with a field strength root-mean-square ${B}^{{\rm{str}}}$ is modeled as two opposing regular fields with same strength of ${B}^{{\rm{str}}}$. Its total (polarized) synchrotron emissivity is ${j}_{{\rm{tot}}({\rm{pol}})}^{{\rm{str}}}(\nu ,{B}_{\perp }^{{\rm{str}}})=2{j}_{{\rm{tot}}({\rm{pol}})}^{{\rm{reg}}}(\nu ,{B}_{\perp }^{{\rm{str}}})$, where ${j}_{{\rm{tot}}({\rm{pol}})}^{{\rm{reg}}}$ is the emissivity of regular field (see equations (24 ) and (25 )) with B⊥ replaced by ${B}_{\perp }^{{\rm{str}}}$, the projection of the magnetic field onto the plane perpendicular to the line of sight. For the isotropic random field with a field strength root-mean-square Bran, the total synchrotron emissivity is modeled as the average emissivity of isotropically distributed regular fields with the same strength of Bran, approximated as in [153]: Galprop scenario under the assumption of α = 3, and compare synchrotron spectra predicted by these two model scenarios. The two scenarios obtain nearly identical results for CR electron spectrum with α ≈ 3. While for several spectra of evaporated all-electrons after propagation, the maximal deviation reaches ∼20% within the frequency range from 1 MHz to 10 GHz. The predictions of the two model scenarios are in agreement, so that we can use the GMF model constructed in [162, 163] to calculate the synchrotron signals of PBHs with the Galprop code.
$\begin{eqnarray}\begin{array}{l}{j}_{{\rm{tot}}}^{{\rm{ran}}}(\nu ,{B}^{{\rm{ran}}})=\frac{2\sqrt{3}{q}_{e}^{3}}{{m}_{e}{c}^{2}}{B}^{{\rm{ran}}}\,{x}^{2}\\ \,\times \,\left[{K}_{4/3}{K}_{1/3}-\frac{3}{5}x({K}_{4/3}^{2}-{K}_{1/3}^{2})\right],\end{array}\end{eqnarray}$
where x = ν/νc with νc = 3qeBran γ2/(2πmec), and K4/3, K1/3 are modified Bessel functions of variable x. On the other hand, in the framework of [162, 163], a constant spectral index α = 3 is assumed for CR electrons, simplifying the synchrotron emissivity to $\begin{eqnarray}\begin{array}{l}J\propto \langle {({B}_{\perp }^{{\rm{reg}}}+{B}_{\perp }^{{\rm{str}}}+{B}_{\perp }^{{\rm{ran}}})}^{2}\rangle \\ \,=\,{B}^{2}{\sin }^{2}\theta +\beta {B}^{2}{\sin }^{2}\theta +\frac{2}{3}{({B}^{{\rm{ran}}})}^{2},\end{array}\end{eqnarray}$
where ⟨ ⟩ means average over the random fields, θ is the angle between the line of sight and the regular field, and β is the striated-to-regular emission ratio defined in the articles. However, for the all-electrons evaporated by PBHs, their energy spectra deviate from α = 3 power laws. We match the model parameters of this scenario with those of the The turbulent magnetic fields remain nearly coherent over the small scale of ∼100 pc, implying that stochastic fluctuations may not fully average out on kpc scale. The model scenarios described above, which capture only the mean effect of the turbulent magnetic fields, neglect their stochastic features. A more elaborate modeling method involves stochastic field realizations. For example, as implemented in the Hammurabi code [164, 165], the turbulent magnetic fields are modeled with realizations of Gaussian stochastic fields with the power spectrum informed by magnetohydrodynamic theory. In such framework, observational data correspond to the prediction of a possible stochastic field realization. Consequently, a precise match of the model prediction to the observational data is not expected. However, such models can predict the level of variation induced by the turbulent magnetic fields, referred to as the the galactic variance. The galactic variance itself is an observable and can be used to quantify the statistical significance of the model residuals [166]. Owing to the high computational cost of stochastic field realizations, few GMF models are constructed with this method, notable examples like [167–169] are restricted to the Galactic plane. For our purpose of computing the synchrotron signals of PBHs, the model capturing the mean feature of turbulent magnetic fields suffices. Nevertheless, the stochastic nature of turbulent magnetic fields can introduce additional uncertainties. For instance, in the sky map with a pixel resolution of 3. 7∘ (with Nside = 16 in HEALPix scheme), fluctuations on the scale of 100 pc cannot average out within ∼1.5 kpc far from the observer. These remained local fluctuations act as foregrounds of observational data, resulting in an unconsidered uncertainty in our analysis.
Despite the modeling difficulties, several GMF models have been constructed under specific assumptions about interstellar gas, CR electron density, and GMF geometrical structure [161]. In this work, we adopt two benchmark models which include all three GMF component fields (regular, striated random, and isotropic random) and achieve agreements with RM and synchrotron (both I and P) observational data, as follows:
SUNE model. The model of the same name constructed in [154]. It utilizes the best-fit regular field from [170, 171], including the ASS + RING disc field and an updated toroidal halo field from [171]. Its striated random field aligns with the regular field and maintains a constant strength ratio relative to the strength of regular field through out the Galaxy. The isotropic random field decreases exponentially in both radial and vertical directions. The random field parameters are fitted using the radio continuum surveys from 22 MHz to 2.3 GHz, and the total and polarized intensity data from 7 year Wilkinson Microwave Anisotropy Probe (WMAP) at frequencies between 23 GHz and 94 GHz. The WMAP MCMC templates are adopted to account for dust and spinning-dust emissions. Free–free absorption and emission are computed with specific interstellar ionized gas model, and the CR electron density follows the z04LMPDS model in [172].
JF12 model. A more complex model constructed in [162, 163]. Its regular field includes: 1. a disc field with eight spiral arms, 2. a toroidal halo field, 3. an X-shaped halo field extending perpendicular to the Galactic plane. The striated random field is parameterized by a parameter β in the whole galaxy, defined as the striated-to-regular emission ratio. We use the updated value of β in [163]. The isotropic random field includes: 1. a disc field mirroring the structure of the regular disc field, 2. a halo field combining radial exponential and vertical Gaussian profiles. The random field parameters are fitted using the total and polarized intensity data from the 7 year WMAP at 23 GHz, and with the assumption of a constant spectral index α = 3 for CR electron density.
5. Constraints on the PBH abundance
5.1. Observational data and error estimation
The diffuse Galactic radio emissions arise from several astrophysical processes. At frequencies below a few GHz, synchrotron radiation produced by CR electrons propagating in the GMF dominates the sky brightness. At higher frequencies, free–free emission from ionized gas and thermal emission from interstellar dust become increasingly important, while anomalous microwave emission from spinning dust grains may also contribute in the GHz range. In this work, we focus on the low-frequency regime, where the observed radio emissions are primarily of synchrotron origin. Moreover, low-frequency radio waves can be affected by other astrophysical processes when passing through the interstellar plasma. In particular, free–free absorption can attenuate the emission near the Galactic plane at frequencies below ∼100 MHz, while Faraday rotation can lead to significant depolarization below ∼1 GHz.
After propagation, the energies of evaporated all-electrons are mostly below a few GeV (see figure 2). Thus, the synchrotron signals of PBHs are primarily below a few GHz. We use the data from radio continuum surveys at frequencies between 22 MHz and 1.4 GHz to constrain fPBH. The datasets include:
22 MHz. Dominion Radio Astrophysical Observatory Northern Hemisphere survey [173].
45 MHz. Northern sky: Japanese Middle and Upper Atmosphere radar array [174], Southern sky: Maipú Radio Astronomy Observatory 45 MHz array [175], and the combined all-sky map from [176].
85 MHz. Parkes radio telescope [177].
150 MHz. Parkes–Jodrell Bank all-sky survey [177].
408 MHz. Bonn–Jodrell Bank–Parkes all-sky survey [178], and reprocessed by [179].
850 MHz. Dwingeloo telescope [180].
1420 MHz. Northern sky: 25 m Stockert telescope [181, 182], Southern sky: 30 m Villa-Elisa telescope [183, 184].
All datasets are publicly accessible on the LAMBDA website.7 The observed intensity is recorded in the sky map of brightness temperature T(ν), which is converted from the intensity I(ν) as $\begin{eqnarray}T(\nu )=\frac{{c}^{2}\,I(\nu )}{2{k}_{{\rm{B}}}\,{\nu }^{2}},\end{eqnarray}$
where c is the speed of light, kB is the Boltzmann constant, and ν is the frequency. It is the temperature that a body with a Rayleight-Jeans spectrum would need to emit the same intensity at the frequency ν. Thus, the brightness temperature is defined in unit of Kelvin (Kbr). To obtain conservative constraints on fPBH, we do not subtract any offset from the observational sky maps.All observational sky maps are pixelated using the HEALPix scheme [156], of which the minimum Nside = 64. To suppress the effect of fluctuations of the turbulent magnetic fields and to estimate the observational errors, we degrade the sky maps to Nside = 16, corresponding to a resolution of approximately 3.7∘. The observational error for each degraded pixel is estimated as the standard deviation of values of all high-resolution sub-pixels contained within the degraded pixel. As discussed in [162, 163], this error estimation method not only accounts for experimental errors, but also includes the astrophysical variance caused by turbulent magnetic fields and inhomogeneities in the interstellar medium. For sky maps with the lowest original resolution (Nside = 64), the standard deviation is computed with 16 sub-pixels contained within each degraded pixel. For robustness of our constraints, we reject pixels (with Nside = 16) which contain sub-pixels with bad-values at the original resolution.
5.2. Synchrotron signals of PBHs
In calculation of synchrotron signals of PBHs, we can safely neglect the extragalactic contributions to the signals and consider only the synchrotron emissions produced within the Galaxy. This is due to the very small averaged magnetic field strength B ∼ nG on cosmological scales [185], which strongly suppresses the emissivity since it scales approximately as ∝B2. In figures 3 and 4, we show sky maps of synchrotron signals of PBHs with NFW and Burkert DM profiles (left and right) at 22 MHz, 150 MHz, and 408 MHz (top to bottom), assuming a monochromatic mass function with Mc = 5 × 1014 g, GH CR propagation model, and SUNE and JF12 GMF models, respectively. All maps are computed with fPBH = 10−8. The 22 MHz and 150 MHz maps correspond to the frequencies at which the strongest single-frequency constraints on fPBH are obtained for the SUNE and JF12 models, respectively (see below). The 408 MHz map is shown because it corresponds to the frequency at which the observational sky map is minimally contaminated by non-synchrotron processes. Dark regions on the Galactic disc towards the Galactic center in the maps at 22 MHz and 150 MHz arise from free–free absorption, which diminishes at higher frequencies and becomes negligible at 408 MHz. It can be seen that Burkert profile predicts weaker synchrotron signals than those predicted by NFW profile towards the region near the Galactic center, where the Burkert profile predicts a smaller DM density.
Figure 3. Sky maps of synchrotron signals of PBHs with NFW and Burkert DM profiles (left and right) at 22 MHz, 150 MHz, and 408 MHz (top to bottom), assuming a monochromatic mass function with Mc = 5 × 1014 g, GH CR propagation model, and SUNE GMF model. All maps are computed with fPBH = 10−8, and are plotted with Nside = 64 and linear color mapping. |
Figure 4. Same as figure 3 but for JF12 GMF model. |
In figure 5, we show the spectra of synchrotron signals of PBHs for GH, DRz4, DRz6 and DBz4 CR propagation models, assuming SUNE GMF model, NFW DM profile, and monochromatic mass functions of selected Mc = 5 × 1014 g, 1 × 1015 g, and 5.6 × 1016 g. The spectra are averaged within the intermediate latitudes 10∘ < ∣b∣ < 40∘. Observed intensities averaged within the same region are also shown. The spectrum of DBz4 model with Mc = 5.6 × 1016 g is much lower than the observed intensities, thus is not shown. For a tight plot, the spectra of GH model are scaled by fPBH = 2 × 10−8, 2 × 10−8 and 2 × 10−4 for each Mc, and the spectra of DBz4 model are scaled by fPBH = 2 × 10−4 and 1.0 for the smaller two Mc (as labeled in the figure). It can be seen that the synchrotron signal intensity of GH model is larger than that of DBz4 model significantly, by a factor of 104 for Mc = 5 × 1014 g at the lowest observational frequency 22 MHz, and the factor becomes much larger for higher frequencies and larger Mc. By requiring the synchrotron signals not exceeding the observed intensities, constraints on fPBH can be roughly estimated. Since the spectra of synchrotron signals are softer than the observed spectrum, we expect that the observational data at lower frequencies can set stronger constraints. For the GH model, strong constraints can be estimated, such as fPBH < 2 × 10−8 for Mc = 5 ×1014 g, and one order of magnitude weaker constraint for Mc = 1 × 1015 g. Even for heavy PBHs with Mc = 5.6 × 1016 g, fPBH < 2 × 10−2 can be estimated. In contrast, for the DBz4 model, only the lightest PBHs with Mc < 1 × 1015 g can be constrained. The estimated constraints are much weaker, for example, fPBH < 2 × 10−4 for Mc = 5 × 1014 g.
Figure 5. Spectra of synchrotron signals of PBHs for GH, DRz4, DRz6 and DBz4 CR propagation models, assuming SUNE GMF model, NFW DM profile, and monochromatic mass functions of selected Mc = 5 × 1014 g, 1 × 1015 g, and 5.6 × 1016 g. The spectra are averaged within the intermediate latitudes 10∘ < ∣b∣ < 40∘. Observed intensities averaged within the same region are also shown. For a tight plot, the spectra of GH model are scaled by fPBH = 2 × 10−8, 2 × 10−8 and 2 × 10−4 for each Mc, and the spectra of DBz4 model are scaled by fPBH = 2 × 10−4 and 1.0 for the smaller two Mc (as labeled). Only the spectra with Mc = 5.6 × 1016 g for DRz4 and DRz6 models are shown, which are also scaled by fPBH = 2 × 10−4 as the spectrum of GH model with the same Mc (the fPBH is not labeled). |
In figure 5, we also show the spectra of synchrotron signals of PBHs for DRz4 and DRz6 models with Mc = 5.6 × 1016 g, which are scaled by fPBH = 2 × 10−4 (not labeled in the figure). Compared to the GH and DRz4 models whose diffusion halo half-height zh = 4 kpc, the DRz6 model predicts larger synchrotron signal intensities for Mc = 5.6 × 1016 g due to its larger zh = 6 kpc. The synchrotron signal intensity of the DRz4 model is larger that of the GH model for the same Mc at frequencies of ${ \mathcal O }(10\,{\rm{MHz}})$, which arise from the larger flux of evaporated all-electrons at energies of ${ \mathcal O }(100\,{\rm{MeV}})$ predicted by the DRz4 model (see figure 2). For higher frequencies of ≳100 MHz, the DRz4 and GH models predict nearly equal synchrotron signal intensities, because the fluxes of evaporated all-electrons predicted by these two models become in agreement at energies of ≳1 GeV.
At this point, we would like to clarify that the inclusion of PBHs as additional CR sources does not affect the self-consistency of the adopted GMF models, which are generally derived under the assumption of a specific CR electron spectrum. The parameters of the GMF model are tightly constrained by RMs and high-frequency radio/microwave observations (e.g. ≳20 GHz for WMAP). Among these observables, the RMs are completely independent of CR models, and the high-frequency radio emissions are dominated by CR electrons with energies above ∼10 GeV. As we show in figure 2, for the all-electrons evaporated by PBHs, their energy cannot reach 10 GeV even though the re-acceleration effect is considered in propagation. Therefore, the inclusion of PBHs as additional CR sources affects the high-energy CR electron spectrum marginally in the high-energy regime relevant for the calibration of GMF model parameters. We also have checked its correctness by extending the synchrotron signal spectra to higher frequencies of ∼20 GHz, and found that the signals drop rapidly above a few GHz and become much smaller than the observed radio emissions.
5.3. Results
We constrain fPBH by requiring that, for each pixel in the observational sky maps (with Nside = 16), the synchrotron signal of PBHs does not exceed the observed intensity plus twice the observational error. Since the synchrotron signals are proportional to fPBH, we rescale fPBH until the limit condition is satisfied. To obtain conservative constraints, we do not attempt to model or subtract any astrophysical background components. Instead, we conservatively allow the entire observed radio intensity to be accounted for by synchrotron signals of PBHs.
In figure 6, we show the constraints on fPBH for different CR propagation models, assuming monochromatic mass functions, SUNE GMF model, and NFW DM profile. For the DBz4 model which does not include the re-acceleration, we obtain modest constraints only for Mc < 1 × 1015 g, as can be expected from figure 5. In contrast, for all diffusive re-acceleration models, we obtain stringent constraints. For instance, the most conservative constraints (obtained with the GH model) require fPBH < 1.25 × 10−4 for Mc = 1.9 × 1016 g. It can be seen that, as the diffusion halo half-height zh increases from 4 kpc to 10 kpc (in the DRz4-10 models), the constraints strengthen by approximately an order of magnitude. For the GH and DRz4 models whose zh = 4 kpc, the obtained constraints differ slightly for Mc ≳ 5 × 1015 g.
Figure 6. Constraints on fPBH for different CR propagation models, assuming monochromatic mass functions, SUNE GMF model, and NFW DM profile. |
In figure 7, we show the constraints on fPBH for SUNE and JF12 GMF models, assuming monochromatic mass functions, GH CR propagation model, and NFW DM profile. The shaded region illustrates the impact of varying the overall normalization of the SUNE magnetic field by a factor of two on the resulting constraints. It can be seen that, although these two GMF models have distinct geometrical structures and spatial distributions of magnetic field strength, the derived constraints are in remarkable agreement. Adopting alternative CR propagation models leads to similar results compared to the GH case. For completeness, we present these results in appendix B . To quantify the impact of magnetic field strength uncertainty, we vary the overall normalization of the SUNE magnetic field by a factor of two. This results in the range of constraints shown by the shaded region in figure 7. In particular, increasing the magnetic field strength by a factor of two strengthens the constraints by approximately a factor of 3.6. Nevertheless, the overall qualitative conclusions remain unchanged.
Figure 7. Constraints on fPBH for SUNE and JF12 GMF models, assuming monochromatic mass functions, GH CR propagation model, and NFW DM profile. The shaded region illustrates the impact of varying the overall normalization of the SUNE magnetic field by a factor of two on the resulting constraints. |
Although the sky maps of the synchrotron signals of PBHs (see figures 3 and 4) for NFW and Burkert DM profiles look different, the constraints obtained with these two DM profiles differ slightly. We show the comparison of the constraints obtained with different DM profiles in appendix B . This agreement can be understood by that the main difference between the sky maps of synchrotron signals comes from the sky region near the Galactic center, where the Galactic synchrotron emission is too bright to set strong constraints. As will be shown later, the constraints are derived from the pixel located in dark areas in the observational sky map, and none of them is close to the Galactic center.
We regard the constraints obtained with the GH model (the blue line in figure 6) as our main conclusion. Because, firstly, the constraints are the most conservative among the results obtained with diffusive re-acceleration models. Secondly, the best-fit zh = 4.0 kpc of the GH model matches the CR electron densities used in the construction of both the SUNE and JF12 GMF models, improving consistency of our calculation.
In figure 8, we show constraints on fPBH derived from individual observational sky maps, assuming monochromatic mass functions, GH CR propagation model, NFW DM profile, and SUNE GMF model in the left, JF12 GMF model in the right, respectively. For the SUNE model, the constraints weaken as the observational frequency increasing, and the strongest constraints are derived from the 22 MHz observational sky map. For the JF12 model, the observational sky maps at frequencies from 45 MHz to 150 MHz set equal and the strongest constraints, while the constraints derived from the 22 MHz observational sky map are weaker slightly. Besides, the moderately weaker constraints derived from the 408 MHz observational sky map are particularly robust. Because the 408 MHz observational sky map is minimally contaminated by non-synchrotron processes, and the observational sky maps at frequencies from 22 MHz to 85 MHz are not full-sky.
Figure 8. Constraints on fPBH derived from individual observational sky maps, assuming monochromatic mass functions, GH CR propagation model, and NFW DM profile. Left: for SUNE GMF model. Right: for JF12 GMF model. |
In figure 9, we show the ratio of synchrotron signals to observational limits (intensity + 2 × error) at 22, 150 and 408 MHz (top to bottom) for SUNE and JF12 GMF models, respectively, assuming a monochromatic mass function with Mc = 5 × 1014 g, GH CR propagation model, and NFW DM profile. Synchrotron signals are scaled by specific values of the constraints on fPBH derived from individual observational sky maps such that ratio = 1 corresponds to the pixel where the constraint is obtained. Observational sky maps at corresponding frequencies are also shown for comparison. Results for the Burkert profile are put in the appendix B , which do not largely differ from those for the NFW profile. From the figures, it can be seen that the Galactic center is too bright to set strong constraints. Instead, the most constraining sky regions are:
For the SUNE model. In the 22 MHz ratio map, the region ∼30∘ above the Galactic center within ∼10∘ radius, corresponding to the dark area enclosed by the North Polar Spur in the observational sky map.
For the JF12 model. In the 150 MHz ratio map, the region ∼15∘ below and above the Galactic disc at longitudes 240∘ < l < 300∘, corresponding to the dark area on the Galactic disc at the same longitudes in the observational sky map.
Note that in the right panel of figure 8, the constraints derived from the 22 MHz observational sky map for the JF12 model are slightly weaker than the strongest constraints, because the outer-disc region (at longitudes 240∘ < l < 300∘) absents in the 22 MHz observational sky map. The difference between the most constraining sky regions for these two GMF models likely origins from JF12's prominent disc random field (especially in 6th and 7th spiral arms), which yields larger synchrotron signals on the Galactic disc.
Figure 9. The first two columns: ratio of synchrotron signals to observational limits (intensity + 2 × error) at 22, 150 and 408 MHz (top to bottom) for SUNE and JF12 GMF models, respectively, assuming a monochromatic mass function with Mc = 5 × 1014 g, GH CR propagation model, and NFW DM profile. Synchrotron signals are scaled by specific values of the constraints on fPBH derived from individual observational sky maps such that ratio = 1 corresponds to the pixel where the constraint is obtained. All maps are plotted with Nside = 16 and linear color mapping. The last column: the 22, 150 and 408 MHz (top to bottom) observational sky maps for comparison, plotted with original resolution and logarithmic color mapping. |
We now extend our analysis to log-normal PBH mass functions, which are often considered more physically motivated in formation scenarios. In figure 10, we show constraints on fPBH with log-normal mass functions. The left panel of figure 10 shows results for different CR propagation models, assuming the distribution width σ = 1, SUNE GMF model, and NFW DM profile. Similar to the monochromatic case, we obtain stringent constraints with diffusive re-acceleration models, while the constraints obtained with the diffusion break model are much weaker. Adopting a log-normal mass function modifies the mass dependence of the constraints compared to the monochromatic case. While the constraints at the lowest masses (e.g. ∼1015 g) remain nearly unchanged (compared with the lines in figure 6), the limits derived with the log-normal mass function weaken more slowly with increasing Mc. As a result, around 1016 g the constraints become stronger than the limits of monochromatic case by approximately one order of magnitude for σ = 1. This behavior arises because the extended mass function retains a significant contribution from lighter PBHs, which emit more and higher-energy all-electrons via Hawking radiation, resulting in stronger synchrotron emission. While for the lowest masses, this effect vanishes as we impose a lower bound of the PBH mass (${M}_{{\rm{\min }}}=5\times 1{0}^{14}\,{\rm{g}}$). The right panel of figure 10 compares the constraints across distribution widths σ = 0.5, 1.0, and 2.0, assuming GH CR propagation model and otherwise identical conditions. Constraints strengthen with increasing σ for Mc > 5 × 1015 g, because broader mass functions (at fixed Mc) retain lighter PBHs, whose larger contribution to synchrotron signals effectively strengthen the constraints. Since the PBH mass function affects the injection spectrum of evaporated all-electrons, rather than their spatial dependence, the sky maps of synchrotron signals remain similar across σ values. The effects of changing different GMF models and DM profiles are similar to those in the monochromatic case, which can be regard as a log-normal limit as σ → 0. So we omit redundant exhibitions.
Figure 10. Constraints on fPBH with log-normal mass functions. Left: for different CR propagation models, assuming the distribution width σ = 1, SUNE GMF model, and NFW DM profile. Right: for different distribution widths σ = 0.5, 1.0, and 2.0, assuming GH CR propagation model, SUNE GMF model, and NFW DM profile. |
In figure 11, we compare our constraints on fPBH derived from radio continuum surveys with monochromatic mass functions (assuming GH CR propagation model, SUNE GMF model, and NFW DM profile) with several existing limits. These include constraints derived from the data of extragalactic γ-rays [52], CMB [26], 511 keV γ-rays [59], Voyager-1 all-electrons [70], and AMS-02 positrons [71]. The result from [59] was obtained assuming an NFW inner slope of γ = 1.6, making it rather conservative. While [60] obtained much stronger constraints from the 511 keV γ-ray data. The results from [70] and [71] were obtained without including the astrophysical backgrounds, and using the same method to calculate the constraints on fPBH as in this work, i.e. requiring that the predicted signal does not exceed the data by more than twice the experimental uncertainty. Thus, we can compare our constraints with them directly. Moreover, our previous constraints derived from AMS-02 positron data [71] were obtained using the same GH CR propagation model.
Figure 11. Several constraints on fPBH with monochromatic mass functions. Shadow regions show constraints derived from different observables, including extragalactic γ-rays [52], CMB [26], and 511 keV γ-rays [59]. Lines show constraints derived from Voyager-1 all-electrons data [70], AMS-02 positron data [71], and radio continuum surveys (assuming GH CR propagation model, SUNE GMF model, and NFW DM profile). |
From figure 11, it can be seen that the constraints derived from radio continuum surveys are stronger than those derived from the Voyager-1 all-electron data [70] by more than one order of magnitude for MPBH ≳ 1 × 1016 g, and also stronger than our previous constraints derived from AMS-02 positron data [71] for MPBH ≳ 2 × 1016 g. The CR positrons are believed to be of secondary origin and relatively rare, making them sensitive to exotic contributions. It is therefore noteworthy that radio continuum surveys can set constraints that are comparable to, or even stronger than, those derived from AMS-02 positron data. The CR particles lose energy when traveling through the heliosphere. Based on the force-field approximation, interstellar positrons at least need a kinetic energy ${E}_{{\rm{k}}}\gtrsim {E}_{{\rm{AMS}},{\rm{\min }}}$ + eφ ≈ 0.57 GeV + 0.5 GeV ≈ 1 GeV to be detected by AMS-02, where ${E}_{{\rm{AMS}},{\rm{\min }}}=0.57\,{\rm{GeV}}$ is the kinetic energy of the lowest-energy AMS-02 positron data point, and φ = 500 MV is a conservatively estimated modulation potential. In contrast, the observation of Galactic synchrotron emissions at frequencies of ∼20 MHz can provide indirect measurements of interstellar CR all-electrons with energies of ∼100 MeV, and therefore can set stronger constraints on heavier PBHs, for which the energy spectra of evaporated all-electrons are shifted to lower energies. For the heaviest PBHs with masses MPBH ≳ 1 × 1017 g, their temperatures become too low to emit significant all-electrons, causing all constraints to weaken. Furthermore, the synchrotron emissions can travel though the Galaxy, thus our constraints can complement those derived from local CR measurements.
It is worth noting what fundamentally limits our current conservative constraints and the implications for future radio probes. At low radio frequencies, the instrumental thermal noise is negligible compared to the extremely bright Galactic synchrotron emission. Therefore, our bounds are not primarily limited by instrumental thermal noise, but rather by the intensity of the astrophysical radio background and its spatial variance (through the observational error). The next-generation radio facilities, such as the SKA-low and improved LOFAR surveys, are expected to meaningfully tighten these constraints. They will achieve this through their unprecedented angular resolution, which can resolve and remove extragalactic point sources, reducing the small-scale spatial variance in the observational maps. Furthermore, their advanced broad-band and polarimetric capabilities will make it possible to construct a robust GMF model and subtract the astrophysical radio background. By constraining fPBH using only the much smaller residual intensity rather than the total observed intensity, future limits could be improved by an order of magnitude. On the other hand, since the predicted spectra of synchrotron signals of PBHs are softer than the observed spectrum, future radio surveys extending below the atmospheric window threshold (10 MHz) could provide stronger constraints. It is also possible to obtain meaningful constraints on fPBH with diffusion break models of CR propagation, if the radio survey data at lower frequencies are available.
6. Conclusion
In this work, we have investigated the possibility of constraining PBHs with masses MPBH ≳ 1015 g through Galactic diffuse synchrotron emissions. Using the AMS-02 [85] and Voyager-1 [86] data on the B/C flux ratio, we fit parameters of a set of benchmark CR propagation models, and confirm that a significant Alfvén velocity Va ∼ 20 km s−1 is favored in several diffusive re-acceleration models. We show that, for Va ∼ 20 km s−1, a significant fraction of the evaporated all-electrons, which typically have initial energies of ∼10 MeV, can be boosted to the energies of ∼100 MeV during their Galactic propagation. Consequently, synchrotron emissions generated by the evaporated all-electrons can be constrained by the low-frequency radio continuum surveys (from 22 MHz to 1.4 GHz). With diffusive re-acceleration models, we obtain stringent constraints on fPBH. The most conservative constraints are stronger than those derived from the Voyager-1 all-electron data [70] by more than one order of magnitude for MPBH ≳ 1 × 1016 g, and also stronger than our previous constraints derived from the AMS-02 positron data [71] for MPBH ≳ 2 × 1016 g. We have checked the robustness of our constraints against different DM profiles and uncertainties in the GMF model, and we have also extended our analysis to log-normal PBH mass functions.
Admittedly, the main results of this work critically depend on the diffusive re-acceleration models of CR propagation, which are supported by a number of independent analyses [74–84], but not yet conclusive. Future x-ray observation missions (e.g. e-ASTROGAM, AMEGO [186]) can probe interstellar all-electrons with energies of ${ \mathcal O }(100{\rm{MeV}})$ via their inverse Compton effects, potentially examining the existence of diffusive re-acceleration [136]. Finally, we emphasize that our present bounds are primarily limited by the intensity and spatial variance of the astrophysical radio background. Future low-frequency radio facilities such as SKA-low and improved LOFAR surveys with improved angular resolution and broad-band polarimetric capabilities are expected to enable more robust modeling of the GMF and subtraction of the astrophysical radio background in the analysis. Such improvements could tighten the constraints on fPBH by an order of magnitude. In addition, since the spectra of synchrotron signals of PBHs are softer than the observed spectrum, future radio surveys extending to frequencies below the atmospheric window threshold can also provide stronger constraints on PBHs.
Appendix A Fit results for benchmark models
We fit parameters of a set of benchmark CR propagation models using the AMS-02 [85] and Voyager-1 [86] data on the B/C flux ratio. Our model set includes four diffusive re-acceleration models (DRz4-10 models) with different diffusion halo half-heights zh (fixed to 4, 6, 8, 10 kpc, respectively), and one diffusion break model (DBz4 model) with zh fixed to 4 kpc. In section 3.3 of the main text, we have introduced the set of model free parameters and the fitting procedure, and shown the the estimated model parameters (see table 1) and the B/C flux ratio predictions (see figure 1) for the benchmark models. In this section, we describe the numerical details of our fits, show the C flux predictions for the benchmark models, and display corner plots (showing the joint and marginal posterior distributions) of the posterior PDFs for the benchmark diffusive re-acceleration models.
A.1. Numerical details
The Galprop configurations are fixed as follows: We employ 2D spatial grids assuming cylindrical symmetry of the Galaxy. The diffusion halo extends radially to R = 20 kpc and vertically to z = ±zh, with zh values of 4, 6, 8, and 10 kpc for the respective models. We set linear spatial grids with 41 vertical grid points and 24 radial grid points, the latter one is adjusted automatically by Galprop because we adopt the vectorized Crank-Nicolson solution method (with GALDEF option solution_method = 4 ). The spatial dependence of the CR primary source follows equation (10 ), with a = 1.9, b = 5.0 [118], and z0 = 0.2 kpc. Primary isotopic abundances, except that of 12C, are fixed to standard Galprop values [187]. The proton flux normalization at 100 GeV is 4.49 × 10−2 GeV−1 m−2 s−1 sr−1, which is used to determine the normalization of primary CR source term. The nuclear reaction network includes elements up to Si (with charge number Z = 14), utilizing the cross section parameterization corresponding to the GALDEF option kopt=012 . For interstellar gas components, we adopt the standard Galprop released CO and H I distributions, and the H II model as in [154]. We use the MultiNest package [150–152] to sample posterior PDFs, of which the dimension is 9 for both types of models. For the MultiNest configurations, we set 400 live points, evidence tolerance of 0.5 and sampling efficiency of 0.8.
A.2. Extended results
In figure 12, we show model predictions and residuals of C flux for the benchmark models. It can be seen that all DRz4-10 models provide equally good fits to Voyager-1 and AMS-02 data, yielding nearly identical predictions. While the DBz4 model exhibits slightly insufficient spectral hardening around 300 GV. In figure 13, we show corner plots of posterior PDFs for DRz4, DRz6, and DRz10 models. The result of DRz8 model is not shown for clarity. The corner plots for the DBz4 model are omitted, since our main conclusion does not rely on this model. For a tight plot, D0 is replaced with D0/zh (with unit: 1028 cm2 s−1 kpc−1). It can be seen that all parameters, except for D0/zh, remain approximately constant across zh. As for D0/zh, it decreases slightly as the zh increasing. This feature is also obtained by other works, for example, see [76].
Figure 12. Model predictions and residuals of C flux for the benchmark models. Left: for Voyager-1. Right: for AMS-02. |
Figure 13. Corner plots of posterior PDFs for DRz4, DRz6, and DRz10 models. The result of DRz8 model is excluded for clarity. Contours show 1 and 2-σ credible regions of 2D distributions, containing 39.3% and 86.4% of the posterior probability mass, respectively. Points mark the posterior means, and crosses mark the maximum posterior values. In 1D marginalized posterior plots, solid vertical lines show the posterior means, and dotted vertical lines show positions of the posterior mean ± standard deviation. For a tight plot, D0 is replaced with D0/zh (with unit: 1028 cm2 s−1 kpc−1). |
Appendix B Extended results for constraints on the PBH abundance
This appendix serves as supplementary materials of section 5.3 . Below, we show several results of constraints on fPBH that are of less importance or redundant to be shown in the main text.
In figure 14, we show constraints on fPBH for different DM profiles and for different CR propagation models (as labeled in the figure), assuming monochromatic mass functions and SUNE GMF model. The plot for DRz8 model is omitted. It can be seen that the constraints obtained with NFW and Burkert profiles are nearly equal for all diffusive re-acceleration models. While for JF12 GMF model, the difference between the constraints obtained with these two DM profiles are less obvious, thus the results are not shown.
Figure 14. Constraints on fPBH for different DM profiles and for different CR propagation models (as labeled), assuming monochromatic mass functions and SUNE GMF model. The plot for DRz8 model is omitted. |
In figure 15, we show constraints on fPBH for different GMF models and for different CR propagation models (as labeled in the figure), assuming monochromatic mass functions and NFW DM profile. The plot for DRz8 model is omitted. It can be seen that the constraints obtained with SUNE and JF12 GMF models are nearly equal for all diffusive re-acceleration models. The results for Burkert DM profile are quite similar, thus are not shown.
Figure 15. Constraints on fPBH for different GMF models and for different CR propagation models (as labeled), assuming monochromatic mass functions and NFW DM profile. The plot for DRz8 model is omitted. |
In figure 16, we show the ratio of synchrotron signals to observational limits (intensity + 2 × error) at 22, 150 and 408 MHz (top to bottom) for SUNE and JF12 GMF models, respectively, assuming the monochromatic mass function with Mc = 5 × 1014 g, GH CR propagation model, and Burkert DM profile. Synchrotron signals are scaled by specific values of the constraints on fPBH derived from individual observational sky maps such that ratio = 1 corresponds to the pixel where the constraint is obtained. Observational sky maps at corresponding frequencies are also shown for comparison. The most constraining sky regions are similar to those for NFW profile, while the outer-disc region becomes more prominent, especially for the SUNE model at frequencies of ≳150 MHz. However, the most constraining sky region in the 22 MHz sky map for the SUNE model still matches the dark area enclosed by the North Polar Spur, as the outer-disc region absents in the 22 MHz observational sky map.
Figure 16. Same as figure 9 but for Burkert DM profile. |