# Limits on primordial black holes from $\mu$ distortions in cosmic microwave background

Tomohiro Nakama, Bernard Carr, Joseph Silk

► To cite this version:

Tomohiro Nakama, Bernard Carr, Joseph Silk. Limits on primordial black holes from $\mu$ distortions in cosmic microwave background. Physical Review D, 2018, 97 (4), pp.043525. 10.1103/Phys-RevD.97.043525. hal-01719926

HAL Id: hal-01719926

https://hal.science/hal-01719926v1

Submitted on 2 Oct 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Limits on primordial black holes from $\mu$ distortions
in cosmic microwave background

Tomohiro Nakama,¹ Bernard Carr,²,³ and Joseph Silk¹,⁴,⁵

¹Department of Physics and Astronomy, Johns Hopkins University,
3400 North Charles Street, Baltimore, Maryland 21218, USA

²School of Physics and Astronomy, Queen Mary University of London,
Mile End Road, London E1 4NS, United Kingdom

³Research Center for the Early Universe (RESCEU), Graduate School of Science,
University of Tokyo, Tokyo 113-0033, Japan

⁴Institut d'Astrophysique de Paris, UMR 7095 CNRS, Sorbonne Université,
98 bis Boulevard Arago, 75014 Paris, France

⁵BIPAC, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, United Kingdom

<figure><img src="image_2.png" /></figure>

(Received 25 October 2017; published 21 February 2018)

If primordial black holes (PBHs) form directly from inhomogeneities in the early Universe, then the number in the mass range $10^5 - 10^{12} M_\odot$ is severely constrained by upper limits to the $\mu$ distortion in the cosmic microwave background (CMB). This is because inhomogeneities on these scales will be dissipated by Silk damping in the redshift interval $5 \times 10^4 \lesssim z \lesssim 2 \times 10^6$. If the primordial fluctuations on a given mass scale have a Gaussian distribution and PBHs form on the high-$\sigma$ tail, as in the simplest scenarios, then the $\mu$ constraints exclude PBHs in this mass range from playing any interesting cosmological role. Only if the fluctuations are highly non-Gaussian, or form through some mechanism unrelated to the primordial fluctuations, can this conclusion be obviated.

DOI: 10.1103/PhysRevD.97.043525

# I. INTRODUCTION

Primordial black holes (PBHs) have been a focus of great interest for nearly 50 years [1–3], despite there still being no definite evidence for them. One reason for this is that only PBHs could be small enough for Hawking radiation to be important [4], those smaller than about $10^{15}$ g having evaporated by now with many interesting cosmological consequences [5]. Recently, however, attention has shifted to PBHs larger than $10^{15}$ g, which are unaffected by Hawking radiation. This is because of the possibility that they provide the dark matter, an idea that goes back to the earliest days of PBH research [6] and has been explored in numerous subsequent works [7–10]. Since PBHs formed in the radiation-dominated era, they are not subject to the well-known big bang nucleosynthesis (BBNS) constraint that baryons can have at most 5% of the critical density [11], which is well below the 25% associated with the dark matter. They should therefore be classed as nonbaryonic and, from a dynamical perspective, behave like any other cold dark matter (CDM) candidate. There is no compelling evidence that PBHs provide the dark matter, but nor is there evidence for any of the more traditional CDM candidates, either from direct searches with underground detectors and particle accelerators or from indirect searches for the expected gamma-ray, neutrino or positron signatures [12].

Even if nonevaporating PBHs do not provide all the dark matter, they could still have interesting cosmological effects.

For example, they have been invoked to explain the heating of the stars in our Galactic disc [13], the seeding of the supermassive black holes in galactic nuclei [14–16], the generation of large-scale structure through Poisson fluctuations [16,17] and the associated generation of an infrared background [18], the reheating and ionization of the Universe [19,20], and the production of r-process elements [21]. More recently, it has been proposed that coalescing PBHs could explain the LIGO gravitational wave bursts [22], although this may only require a small fraction of the dark matter to be in PBHs [23]. The detection of four black holes with mass around $30 M_\odot$ has come as a surprise to stellar evolution modelers, so it is natural to consider more exotic types of black holes. The suggestion that LIGO could detect gravitational waves from a population of binary intermediate-mass black holes was originally proposed in the context of the Population III scenario by Bond and Carr [24], and—rather remarkably—a paper in 2014 predicted a Population III coalescence peak at $30 M_\odot$ [25]. Since Population III stars are baryonic, such black holes could not provide the dark matter, but this would not preclude intermediate-mass PBHs from doing so. There have been a large number of recent papers on this topic, but the suggestion that there could be a stochastic background of gravitational waves from PBHs goes back a long way [26,27].

There are other possible explanations for these effects, so they do not necessarily require the existence of PBHs.
However, all these effects can be used to place interesting constraints on the number of PBHs, and this in turn constrains the cosmological models which generate them. The constraints are most usefully expressed as limits on the fraction $f(m)$ of the dark matter in PBHs of mass $m$ and have recently been summarized in Refs. [5,9]. Taken together, they suggest that there are only a few mass windows where PBHs could provide all the dark matter ($f=1$): the intermediate-mass range ($10-100 M_{\odot}$), the lunar-mass range ($10^{20}-10^{24} \text{ g}$), the asteroid-mass range ($10^{16}-10^{17} \text{ g}$) and Planck mass relics of evaporation ($10^{-5} \text{ g}$). Even some of these windows may be excluded for a monochromatic PBH mass function, so recently it has been suggested that the PBH dark matter proposal may require an extended mass function [9]. However, there is some dispute over whether this helps or hinders the proposal for the various mass windows [28–30]. In any case, it seems clear that all the dark matter could comprise PBHs only if they were smaller than about $10^2 M_{\odot}$.

PBHs larger than this might still have an appreciable cosmological density, so it is important to consider whether there is a maximum possible mass for a PBH. Since a PBH forming at a time $t$ after the big bang is expected to have a mass of order the particle horizon size ~$10^5 (t/s) M_{\odot}$, this depends on how late a PBH can form [31]. It is sometimes argued that this should be before weak freeze-out at 1 s, corresponding to a maximum mass of $10^5 M_{\odot}$, since otherwise BBNS would be affected. This is because PBH production usually requires large inhomogeneities, and this might be expected to disturb the usual BBNS scenario [32]. However, this argument is not clear cut because the fraction of the Universe in PBHs at a time $t$ after the big bang is only ~$10^{-6} \Omega_{\text{PBH}}(t/s)^{1/2}$, where $\Omega_{\text{PBH}}$ is the current PBH density in units of the critical density [31], so this would be at most $10^{-6}$ at weak freeze-out.

Various limits can be imposed on the density of such large PBHs. As reviewed in Ref. [9], the microlensing of stars has been sought in a wide variety of contexts, and this limits $f(M)$ over various mass ranges below about $10 M_{\odot}$. Above this mass, numerous dynamical effects come into play. In particular, PBHs in the intermediate-mass range would disrupt wide binaries in the Galactic disc. It was originally claimed that this would exclude objects above 400 $M_{\odot}$ [33]. However, more recent studies may reduce this mass [34], so the narrow window between the microlensing and wide-binary bounds is shrinking. Nevertheless, this suggestion is topical because PBHs in the intermediate-mass range might also explain the sort of massive black hole mergers observed by LIGO. More recent studies have examined the induced gravitational wave background from the fluctuations with amplitude below the PBH threshold [35,36].

Could such huge PBHs be *expected* to form? Only a very low production efficiency is needed to generate a significant PBH density today, and they may be generated by three mechanisms [37]: (i) the collapse of large-amplitude

inhomogeneities [31]; (ii) a temporary softening of the equation of state (even if the inhomogeneities are small) [38–41]; (iii) some form of cosmological phase transition [5]. The last two mechanisms are unlikely to be relevant after 1 s, but the first mechanism could be. For example, the inflationary models discussed in Refs. [42–47] could produce a spike in the power spectrum of density fluctuations at a mass scale which is not uniquely specified. It could be tuned to the intermediate-mass range ($10-10^3 M_{\odot}$) if one wants to explain the dark matter [48], but it could, in principle, be much larger. Another interesting formation mechanism, involving the collapse of inflationary bubbles, has recently been explored in Refs. [49,50].

What would be the possible cosmological consequences of a relatively small density of large PBHs? It has been claimed that some observational anomalies may require the existence nonlinear structures early in the history of the Universe [51]. For example, it is now known that most galactic nuclei contain supermassive black holes (SMBHs), extending from around $10^6 M_{\odot}$ to $10^{10} M_{\odot}$ and already in place by a redshift of about 10 [52]. However, it is hard to understand how such enormous black holes could have formed so early unless there were already large seed black holes well before galaxy formation, with these subsequently growing through accretion. Indeed, pregalactic SMBHs might help to *explain* galaxy formation, either by acting as condensation nuclei on account of their gravitational Coulomb effect [53] or through the Poisson fluctuations in their number density [54]. As discussed by Carr and Silk [16], these seeds might well be PBHs.

Another interesting consequence—and the main focus of this paper—arises because the formation of PBHs through mechanism (i) requires density fluctuations, and the dissipation of these fluctuations after $10^6 \text{ s}$ by Silk damping generates a $\mu$ distortion in the CMB spectrum [55]. This applies to fluctuations which fall within the photon diffusion scale during the redshift interval $5 \times 10^4 < z < 2 \times 10^6$, and gives an upper limit $\delta(M) < \sqrt{\mu} \sim 10^{-2}$ over the mass range $10^4 < M/M_{\odot} < 10^{13}$. When the PBH formation probability is relatively large, the dispersion of primordial fluctuations is also expected to be large, so, in principle, the observational limits on the $\mu$ distortions translate into upper limits on the PBH abundance. This constraint was first mentioned in Ref. [56], based on a result in Ref. [57], but the limit on $\mu$ is now much stronger [58].

It should be stressed that this is a limit on the fluctuations from which the PBHs derive, and it can only be translated into a limit on the PBHs themselves if one assumes a model for their formation. If the fluctuations are Gaussian and the PBHs form on the high-$\sigma$ tail, as in the simplest scenario [31], one can infer a constraint on $f(M)$ over a wide mass range. Indeed, a few years ago Kohri et al. [59] obtained an upper limit of $10^5 M_{\odot}$ based on this assumption. However, the Gaussian assumption may be incorrect. For example, Nakama et al. [60] have proposed a “patch” model, in which the relationship between the background inhomogeneities
and the overdensity in the tiny fraction of the volume which collapses to PBHs is modified. Also, the production mechanism advocated in Refs. [49,50] may entail no density perturbations outside the PBHs at all. The $\mu$-distortion constraint could thus be much weaker, depending on the degree of non-Gaussianity of the primordial fluctuations. A phenomenological description of such non-Gaussianity was introduced in Ref. [60] and involves a parameter $p$, such that—for a fixed PBH formation probability—the dispersion of the primordial fluctuations becomes smaller as $p$ is reduced from its Gaussian value of 2, thereby reducing the $\mu$ distortion.

None of these previous works used the $\mu$ limits to directly constrain the PBH mass fraction, so it is not clear whether this precludes the low PBH density required in some cosmological proposals.<fnref n="1" /> In this paper we address this problem by calculating the constraints on $f(M)$ explicitly, using both the FIRAS limit of $\mu = 9 \times 10^{-5}$ [62] and the projected upper limit of $\mu < 3.6 \times 10^{-7}$ from PIXIE [63]. We find that the $\mu$ distortion is predominantly determined by fluctuations with $30 \text{ Mpc}^{-1} < k < 5000 \text{ Mpc}^{-1}$, and this corresponds to the PBH mass range of $10^5 M_{\odot} < M < 10^{10} M_{\odot}$. The mass range around $10^9 M_{\odot}$ is especially restricted, and we would need huge non-Gaussianity if such massive PBHs were to evade the $\mu$-distortion constraints. Alternatively, one could argue that the PBHs formed with initial mass below $10^5 M_{\odot}$ and then underwent substantial accretion. The $\mu$ constraint would then be avoided altogether.

The plan of this paper is as follows: In Sec. II we derive the form of the $\mu$ constraint on the fraction of the dark matter in PBHs on the assumption that they form from Gaussian or non-Gaussian primordial fluctuations with a monochromatic power spectrum. In Sec. III we extend the analysis to include the $y$-distortion limit, and we discuss how our conclusions depend on the type of non-Gaussianity. We also discuss briefly how our conclusions are modified if the fluctuations are nonmonochromatic or the PBHs accrete, both of which would be expected in a more realistic scenario. We summarize our conclusions in Sec. IV.

## II. DEPENDENCE OF PBH LIMITS FROM $\mu$ ON NON-GAUSSIANITY

Primordial fluctuations over a wide range of mass scales will be dissipated by Silk damping when they fall within the photon diffusion scale. The diffusion scale at any epoch is the geometric mean of the horizon mass and the mass of unit optical depth to Thompton scattering:

$$ M_D \sim \sqrt{M_{\tau} M_H} \sim \begin{cases} 10^{10} (t/t_{eq})^{7/4} M_{\odot} & (t < t_{eq}) \\ 10^{13} (t/t_{\text{dec}})^{11/6} M_{\odot} & (t_{eq} < t < t_{\text{dec}}). \end{cases} \quad (1) $$

The masses here refer to the radiation content, which is much larger than the matter content before the time of matter-radiation equality ($t_{eq} \sim 10^{10}$ s). The value of $10^{13} M_{\odot}$ at decoupling ($t_{\text{dec}} \sim 10^{12}$ s) corresponds roughly to the Silk mass, the matter and radiation densities then being comparable.

The associated heat production could affect the CMB in various ways, depending on the epoch, so we first discuss these more general limits. Photons generated by the dissipation of fluctuations before $t_1 = 7 \times 10^6$ s will be completely thermalized, leading to an increase in the photon-to-baryon ratio $S$ of the CMB through acoustic reheating [64]. In principle, this places constraints on fluctuations on scales below $10^5 M_{\odot}$ since $S$ cannot increase much after BBNS, but these limits are relatively weak. Photons generated by the dissipation of fluctuations between $t_1$ and $t_2 \approx 3 \times 10^9$ s will lead to a $\mu$ distortion in the CMB. Photons generated by the dissipation of fluctuations between $t_2$ and $t_3 \approx 3 \times 10^{12}$ s (decoupling) will lead to a $y$ distortion. Note that these effects can also be used to limit the number of evaporating PBHs [65], those in the mass ranges $10^{11} - 10^{12}$ g and $10^{12} - 10^{13}$ g evaporating in the periods $t_1 < t < t_2$ and $t_2 < t < t_3$, respectively. The source of the photons is different—Hawking radiation versus dissipation of fluctuations—but the interaction with the CMB is the same.

Fluctuations on a given comoving scale produce PBHs when that scale enters the horizon, and only much later are fluctuations on the same scale dissipated. To find the PBH mass ranges relevant to the three types of distortion mentioned above, we must find the PBH mass $M$ corresponding to the diffusion mass $M_D$. These masses are different because a PBH forms when the Universe is radiation dominated. Thereafter its mass remains constant, but the radiation mass of a comoving region outside the black hole decreases as $a^{-1}$, which scales as $t^{-1/2}$ for $t < t_{eq}$ and $t^{-2/3}$ for $t > t_{eq}$. One can show that the PBH mass associated with the diffusion mass $M_D$ is

$$ M \sim \begin{cases} 10^2 (M_D/M_{\odot})^{6/7} M_{\odot} & (t < t_{eq}) \\ 10 (M_D/M_{\odot})^{10/11} M_{\odot} & (t_{eq} < t < t_{\text{dec}}) \end{cases} \quad (2) $$

This shows that the PBH mass ranges constrained by the $S$, $\mu$ and $y$ observations are $M/M_{\odot} < 10^5$, $10^5 < M/M_{\odot} < 10^{11}$ and $10^{11} < M/M_{\odot} < 10^{13}$, respectively.

We now focus on the $\mu$ distortions. In Ref. [60], the following phenomenological description of the non-Gaussian probability density function (PDF) of the curvature perturbation $\zeta$ was introduced:

$$ P(\zeta) = \frac{1}{2\sqrt{2\tilde{\sigma}\Gamma(1+1/p)}} \exp \left[ - \left( \frac{|\zeta|}{\sqrt{2\tilde{\sigma}}} \right)^p \right]. \quad (3) $$

This satisfies $\int_{-\infty}^{\infty} P(\zeta) d\zeta = 1$ and reduces to the Gaussian distribution of Ref. [31] when $p = 2$. The dispersion is

<footnote n="1">¹The mass range associated with the $\mu$ limit is indicated in Fig. 1 of Ref. [61] but not the limit on $f(M)$ itself.</footnote>
$$ \sigma^2 = \int_{-\infty}^{\infty} \zeta^2 P(\zeta) d\zeta = \frac{2\Gamma(1+3/p)}{3\Gamma(1+1/p)} \tilde{\sigma}^2, \quad (4) $$

where $\Gamma(a)$ is the gamma function. In particular, $\sigma = \tilde{\sigma}$ when $p=2$, as expected.

We estimate the fraction of the Universe collapsing into PBHs to be

$$ \beta = \int_{\zeta_c}^{\infty} P(\zeta) d\zeta = \frac{\Gamma(1/p, 2^{-p/2}(\zeta_c/\tilde{\sigma})^p)}{2p\Gamma(1+1/p)}, \quad (5) $$

where $\zeta_c$ is the threshold for PBH formation and $\Gamma(a, z)$ is the incomplete gamma function. For $p=2$ this reduces to

$$ \beta = 2^{-1} \text{erfc}(2^{-1/2} \zeta_c / \sigma). \quad (6) $$

This shows that increasing $\mu$ sensitivity is equivalent to increasing $p$. Our results are very sensitive to the value of $\zeta_c$, but there is some ambiguity about this since it depends upon the perturbation profile, pressure gradients playing an important role [66]. Therefore the threshold can only be specified in terms of some range. If we use the value of $\zeta$ at the peak of the perturbation, $\zeta_c$ lies in the range $0.67-1.05$ according to numerical simulations [67]. But this peak value is subject to what is called the environmental effect [68,69]. Harada et al. [67] suggest the range $0.95-1.26$ by approximately converting the threshold in terms of the density perturbation in the comoving slice to $\zeta$. To make our limits as conservative as possible, we use the threshold $\zeta_c \approx 0.67$ obtained for the class of perturbation profiles considered in previous numerical simulations [67]. Equation (5) can be inverted to give

$$ \tilde{\sigma} = \frac{2^{-1/2} \zeta_c}{[Q^{-1}(1/p, 2\beta)]^{1/p}}, \quad (7) $$

where $Q^{-1}(a, z)$ is the inverse of the regularized incomplete gamma function $Q(a, z) = \Gamma(a, z)/\Gamma(a)$, so that $z = Q^{-1}(a, s)$ if $s = Q(a, z)$.

Let us consider the following dimensionless delta-function power spectrum:

$$ \mathcal{P}_{\zeta} = \sigma^2 k \delta(k - k_*) \quad (8) $$

which leads to the $\mu$ distortion [70]

$$ \mu \approx 2.2\sigma^2 \left[ \exp\left(-\frac{\hat{k}_*}{5400}\right) - \exp\left(-\left[\frac{\hat{k}_*}{31.6}\right]^2\right) \right], \quad (9) $$

where $\hat{k}_*$ is the wave number in units of Mpc$^{-1}$. The wave number and the PBH mass are related via [35]

$$ k \approx 7.5 \times 10^5 \gamma^{1/2} \text{ Mpc}^{-1} \left( \frac{g}{10.75} \right)^{-1/12} \left( \frac{M}{30 M_{\odot}} \right)^{-1/2}, \quad (10) $$

where $\gamma$ gives the size of the PBH in units of the horizon mass at formation and $g$ is the number of degrees of freedom of relativistic particles. We set $\gamma = 1$ and $g = 10.75$ for all masses for simplicity.<fnref n="2" /> The initial abundance $\beta$ is related to $f = \Omega_{\text{PBH}}/\Omega_{\text{DM}}$, where $\Omega_{\text{DM}} \approx 0.27$ [71], via [35]

$$ \beta \approx 1.1 \times 10^{-8} \gamma^{-1/2} \left( \frac{g}{10.75} \right)^{1/4} \left( \frac{\Omega_{\text{DM}}}{0.27} \right)^{-1} \left( \frac{M}{30 M_{\odot}} \right)^{1/2} f. \quad (11) $$

This is just Eq. (2.5) of Ref. [5] with $M$ normalized to the mass $30 M_{\odot}$ indicated by the LIGO events. Then, using both the FIRAS limit ($\mu < 9 \times 10^{-5}$) [62] and the projected upper limit of $\mu < 3.6 \times 10^{-7}$ for PIXIE [63], the above phenomenological model of non-Gaussianity can be used to calculate the $\mu$-distortion constraints on $\beta$ or $f$. The results for $f$ are presented in Fig. 1; the results for $\beta$ are not shown explicitly, but they have a similar form and can be inferred from Eq. (11). The figure shows that the limits are highly restrictive but much less so in the non-Gaussian cases than the Gaussian one. The PIXIE constraint may no longer be relevant since the project has not been approved. Therefore we give constraints for both PIXIE and a hypothetical future experiment that we dub HYPERPIXIE, assumed to give $\mu < 10^{-9}$.

The term in square brackets in Eq. (9) peaks with a value of $1$ at $\hat{k}_* = 80$, which corresponds to a mass $2.6 \times 10^9 M_{\odot}$, so the $\mu$ distortion is most sensitive to modes on this scale and has a value $\mu \approx 2.2\sigma^2$. For $\zeta_c = 0.67$, this implies $\zeta_c/\sigma \approx (100, 1.7 \times 10^3, 3.1 \times 10^4)$ for $\mu = (9 \times 10^{-5}, 3.6 \times 10^{-7}, 10^{-9})$. Hence the decimal logarithm of the PBH formation probability in the Gaussian case, given by Eq. (6), formally reaches $(-2.4 \times 10^3, -6.0 \times 10^5, -2.1 \times 10^8)$. Although this limit is very strong, there is also a value of $f(M)$ which corresponds to having just one PBH per current Hubble horizon. This has been described as the “incredulity limit” [72] and can be written as

$$ f_1 = \frac{M}{\Omega_{\text{DM}} M_H} = 8.2 \times 10^{-14} \left( \frac{\Omega_{\text{DM}}}{0.27} \right)^{-1} \left( \frac{h}{0.67} \right) \left( \frac{M}{10^9 M_{\odot}} \right), \quad (12) $$

where $M_H$ is the current horizon mass,

$$ M_H = \frac{c^3}{2GH} = 4.5 \times 10^{22} M_{\odot} \left( \frac{h}{0.67} \right)^{-1}. \quad (13) $$

The $\mu$ limit can be well below this, so we also plot $f_1$ in the figures for comparison.

<footnote n="2">²Strictly speaking, $g$ should be somewhat smaller for $M > 10^5 M_{\odot}$, but this simplification does not cause significant errors due to the weak dependence on $g$ in Eq. (10).</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Upper limits on $f = \rho_{PBH}/\rho_{DM}$ for different values of the non-Gaussianity parameter $p$ [Eq. (3)]. The solid curves are for FIRAS ($\mu = 9 \times 10^{-5}$) and correspond to $p = 0.5, 1, 2$ (from top to bottom). The dashed and dotted curves are for PIXIE ($\mu = 3.6 \times 10^{-7}$) and HYPERPIXIE ($\mu = 10^{-9}$), respectively, with the same values of $p$ from top to bottom. The SINGLE line corresponds to having one PBH per current Hubble volume.</figcaption></figure>

Figures 2 and 3 show the corresponding results for the quadratic and cubic non-Gaussianity investigated in the context of PBHs in Refs. [35,73]:

$$ \zeta = \zeta_G + \frac{3}{5} f_{\text{NL}} (\zeta_G^2 - \sigma_G^2), \quad (14) $$

$$ \zeta = \zeta_G + g \zeta_G^3, \quad g = \frac{9}{25} g_{\text{NL}}. \quad (15) $$

The dispersion $\sigma^2 = \langle \zeta^2 \rangle$ is then related to $\sigma_G^2 = \langle \zeta_G^2 \rangle$ via

$$ \sigma^2 = \sigma_G^2 + 2 \left( \frac{3}{5} f_{\text{NL}} \right)^2 \sigma_G^4, \quad (16) $$

$$ \sigma^2 = \sigma_G^2 + 6g\sigma_G^4 + 15g^2\sigma_G^6. \quad (17) $$

<figure><img src="image_4.png" /><figcaption>FIG. 2. Upper limits on $f = \rho_{PBH}/\rho_{DM}$ for quadratic non-Gaussianity with $f_{NL} = \infty, 10, 0$ from top to bottom. The FIRAS, PIXIE, HYPERPIXIE and SINGLE curves are as for Fig. 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Upper limits on $f = \rho_{PBH}/\rho_{DM}$ for cubic non-Gaussianity with $g_{NL} = \infty, 1000, 0$ from top to bottom and the same values of $g_{NL}$. The FIRAS, PIXIE, HYPERPIXIE and SINGLE curves are as for Fig. 1.</figcaption></figure>

For each $f_{NL}$ or $g_{NL}$ and mass $M$, which can be translated into $k_*$, the dispersion $\sigma^2$ (or equivalently $\sigma_G^2$) corresponding to the upper limit on $\mu$ is obtained. The corresponding value of $f$ can then be obtained using Refs. [35,73], which leads to Figs. 2 and 3. The limits for cubic non-Gaussianity are weaker than those for quadratic non-Gaussianity, and the limits for $p$-type non-Gaussianity are even weaker for sufficiently small $p$. As noted in Ref. [35], in the limit $f_{NL} = \infty$, we have

$$ \beta = \operatorname{erfc}[(\tilde{\zeta}_c/2)^{1/2}] \quad \text{with} \quad \tilde{\zeta}_c = \sqrt{2\zeta_c}/\sigma + 1 \quad (18) $$

where $\operatorname{erfc}(x)$ denotes the complementary error function. In the limit $g_{NL} = \infty$, we have

$$ \beta = \frac{1}{2} \operatorname{erfc} \left[ \frac{1}{\sqrt{2}} \left( \frac{\zeta_c}{\sigma/\sqrt{15}} \right)^{1/3} \right]. \quad (19) $$

These relations and Eqs. (9) and (10) lead to the curves of Figs. 2 and 3 with $f_{NL} = \infty$ and $g_{NL} = \infty$, which serves as a check.

### III. DISCUSSION

Although we have focused mainly on the $\mu$ distortions associated with PBH formation, we should also comment briefly on the $y$ distortions, which can be calculated using Eq. (9b) in Ref. [70] and assuming $y < 1.5 \times 10^{-5}$ [62]. It is clear from Fig. 4, which shows the $\mu$ and $y$ limits for FIRAS with different values of $p$, that there is a transition at around $10^9 M_{\odot}$ above which the $y$ constraint dominates. This scale is independent of the value of $p$ and can be understood from the qualitative discussion at the start of Sec. II. The $y$ curves are interesting, even if one does not expect PBHs larger than $10^9 M_{\odot}$ in practice. Note that including other sources of $y$ distortion would merely strengthen the PBH limits. Indeed, HYPERPIXIE will
<figure><img src="image_6.png" /><figcaption>FIG. 4. Upper limits on $f = \rho_{\text{PBH}}/\rho_{\text{DM}}$ from $\mu$ (blue) and $\gamma$ (red) FIRAS observations for different values of $p$. The Gaussian case corresponds to $p=2$. Also shown are the PIXIE and HYPERPIXIE $\mu$ limits and the line corresponding to one PBH per current Hubble volume.</figcaption></figure>

be looking for the $\gamma$ signal from astrophysical sources at the level of $10^{-6}$ to $10^{-8}$, and such a signal will be obscured if the $\gamma$ from Silk damping exceeds this. The total constraint for given $p$ comes from the combination of the $\mu$ and $\gamma$ limits. The value of $f$ at the intersect gives a local maximum in the constraint, but there are extended minima on either side of this, with the $\mu$ limit being somewhat weaker than the $\gamma$ limit. The value of $f$ at the maximum decreases as $p$ increases. For $p > 2$, the constraints would be even tighter than in the Gaussian case.

At the present epoch we require SMBHs with mass $10^6 - 10^{10} M_{\odot}$ in galactic nuclei, and the Magorrian scaling implies $f_{\text{SMBH}} \sim 10^{-4}$ [74]. If there were no accretion, the PBH scenario would be excluded by the above analysis unless $p$ were as low as 0.5 or $f_{\text{NL}}$ and $g_{\text{NL}}$ as high as 5000. However, the values of $M$ and $f$ indicated in Figs. 1–4 correspond to the formation epoch of the PBHs, and the production of SMBHs in galactic nuclei probably requires a very large accretion factor. For example, for Eddington-limited growth, the initial seed mass could be smaller by a factor of $10^5 - 10^6$, corresponding to an initial $f$ of $10^{-9} - 10^{-10}$. If the final mass is much larger than the initial mass, one expects $f \propto M$, so the value of $(f,M)$ in Figs. 1–4 evolves along a straight line. More generally, since the $\mu$ limit is only important above $10^5 M_{\odot}$, one could circumvent it altogether, provided the mass increases by a factor $(m/10^5 M_{\odot})^{-1}$. However, in this case the accretion of gas onto the PBHs may also contribute to temperature and polarization anisotropies and $\gamma$ distortions of the CMB [19,20,75], so accretion may not save the scenario.

The above analysis has considered the types of non-Gaussianity associated with the parameters $p, f_{\text{NL}}$ and $g_{\text{NL}}$, so it may be useful to comment on the connection between these. The parameter $p$ changes the exponential falloff of the PDF, so that it declines as $\exp(-|\zeta|^p)$ rather than

$\exp(-|\zeta|^2)$, as in the Gaussian case. This is slower than the Gaussian decline for $p < 2$, so for fixed $\sigma$ this widens the PDF but lowers the central height since area is conserved. For $p < 2$ and the same value of $\sigma$, the modified curve goes above the Gaussian curve at a value of $\zeta$ less than the critical value $\zeta_c = 0.67$ required for PBH formation, so PBH production will be increased. This explains why the constraints on $f$ become weaker as $p$ decreases below 2. Of course, this is a particularly simple type of non-Gaussianity, and we do not currently know what kind of physics in the early Universe would lead to it; thus, this is just a convenient phenomenological toy model. Nevertheless, this illustrates qualitatively why the deviations from Gaussianity that we have investigated weaken the PBH constraints.

The parameters $f_{\text{NL}}$ and $g_{\text{NL}}$ are more traditional measures of non-Gaussianity, whose significance is well understood. From Eq. (14), $f_{\text{NL}}$ couples with $\zeta_G^2 - \sigma_G^2$, so that the average $\langle\zeta\rangle$ is zero, while $g_{\text{NL}}$ couples with $\zeta_G$ from Eq. (15). In principle, one could combine Eqs. (14) and (15) into a single equation to allow $f_{\text{NL}}$ and $g_{\text{NL}}$ to vary together. Generally a pure $f_{\text{NL}}$ can yield either a weaker or stronger constraint than some mixture of $f_{\text{NL}}$ and $g_{\text{NL}}$, depending on their precise values. The relative effects of $p, f_{\text{NL}}$ and $g_{\text{NL}}$ are shown in Fig. 5, where the values of $\sigma$ as functions of the non-Gaussian parameters are compared for $f = 10^{-4}$, corresponding to $\beta \approx 2 \times 10^{-10}$ for $10^6 M_{\odot}$.

It should be stressed that there are many other types of non-Gaussianity, so the above analysis should only be regarded as demonstrating the importance of this property in principle for the form of the $\mu$ constraint on PBHs. In particular, Ref. [60] discusses an inflationary scenario—termed the patch model—which can produce a sufficient

<figure><img src="image_7.png" /><figcaption>FIG. 5. The values of $\sigma$ as a function of the non-Gaussian parameters $p, f_{\text{NL}}$ and $g_{\text{NL}}$ for $f = 10^{-4}$. The function $\sigma(p)$ can take a wider range of values than $\sigma(f_{\text{NL}})$ and $\sigma(g_{\text{NL}})$. The Gaussian PDF corresponds to $p = 2, f_{\text{NL}} = 0$ and $g_{\text{NL}} = 0$. Note that $\sigma(f_{\text{NL}} = \infty) \approx 0.024$, and the PDF with $p \approx 0.7$ gives a similar value. The PDF is not an even function for quadratic non-Gaussianity, unlike the PDFs specified by $p$ and $g_{\text{NL}}$. Also $\sigma(g_{\text{NL}} = \infty) \approx 0.01 \approx \sigma(p = 0.5)$, so the PDFs with $g_{\text{NL}} = \infty$ and $p = 0.5$ would look very similar.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Upper limits on $f$ for extended fluctuation spectra, assuming the FIRAS upper limit $\mu = 9 \times 10^{-5}$. The width of the spectrum is taken to be $\sigma_* = (0, 0.5, 1)$. Note that the spectrum is reduced to a delta-function spectrum for $\sigma_* = 0$.</figcaption></figure>

number of PBHs to account for the observed high-redshift quasars, while keeping $\sigma$ sufficiently small to evade constraints from $\mu$ distortions. This model is highly non-Gaussian and essentially introduces a spike in the PDF. This form of non-Gaussianity is very different from the types investigated here, so we have not tried to relate them.

Although we have assumed a delta-function primordial power spectrum, we can generalize our analysis to the more plausible situation in which the fluctuations have an extended power spectrum. As a first step, let us consider the following power spectrum:

$$P_{\zeta} = \sigma^2 k \frac{1}{\sqrt{2\pi}\sigma_* k_*} \exp\left(-\frac{(k-k_*)^2}{2\sigma_*^2 k_*^2}\right). \quad (20)$$

This reduces to the delta-function spectrum in the limit $\sigma_* \to 0$. The $\mu$ distortion can then be calculated using the formalism of Ref. [70] and, for each $\sigma_*$ and $k_*$, the value of $\sigma$ corresponding to the upper limit on $\mu$ can be obtained. In the Gaussian case, $\beta$ is calculated using Eq. (6), and we translate this to $f$ using the correspondence between $M$ and $k_*$ implied by Eq. (10). The upper limit on $f(M)$ so obtained is plotted in Fig. 6. The limits for $M \lesssim 10^5 M_{\odot}$ ($M \gtrsim 10^{12} M_{\odot}$) are tighter when $\sigma_*$ is larger, due to the additional presence of modes with $k < k_*$ ($k > k_*$) which contribute to $\mu$ relatively efficiently but are absent in the corresponding case for a delta-function spectrum.

This is distinct from the effect of an extended PBH mass function, which may arise even for a delta-function power spectrum if one has critical collapse [9]. In this case, one needs to consider which value of $M$ dominates the $\mu$ distortion. This depends upon the form of the PBH mass function; however, the limit will be important whenever

this extends above $10^5 M_{\odot}$, and it will be strongest if it extends up to $10^9 M_{\odot}$. The two most plausible mass functions are lognormal (which applies for a large class of inflationary scenarios) and a steeply rising power law with an exponential upper cutoff (which applies for critical collapse). In the first case the mass function only extends over a few decades of mass, so the above analysis is still a reasonable guide. In the second case, the mass function is essentially monochromatic since the number of PBHs on the power-law tail is very small. If the mass function has a power-law form, so that it is very extended, the situation is more complicated, but this is less likely to apply in any realistic scenario. As discussed in Ref. [76], all the constraints on $f(M)$—not just the $\mu$ constraints—are modified anyway for an extended PBH mass function, but we do not consider this complication further here.

## IV. SUMMARY

PBH formation is a plausible process in the very early Universe. It is motivated observationally by the need for dark matter and by LIGO observations, and theoretically by generic inflationary fluctuation considerations. An inevitable consequence of PBH formation is that the dissipation of primordial fluctuations is also expected to be large. We have developed alternative formalisms for studying primordial non-Gaussianities in these fluctuations and evaluated the associated $\mu$ limits on $f(M)$. Our formalism is equivalent to the more conventional $f_{\text{NL}}$ approach in other recent studies of PBH formation from non-Gaussian fluctuations [73]. We have demonstrated that Silk damping could produce unacceptably large $\mu$ distortions of the cosmic microwave background since the fluctuations with wave numbers corresponding to the PBH masses under consideration are expected to dissipate during the redshift interval $5 \times 10^4 \lesssim z \lesssim 2 \times 10^6$. We have explored current and future constraints on these $\mu$ distortions and shown how they can be translated into upper limits on the PBH abundance over a wide mass range.

## ACKNOWLEDGMENTS

We thank T. Harada, M. Kamionkowski, K. Kohri, V. Poulin, A. Vilenkin, J. Yokoyama and an anonymous referee for helpful input. T. N. and B. C. thank the Research Center for the Early Universe (RESCEU) at the University of Tokyo and the Institut d'Astrophysique de Paris for hospitality received during this work. T. N. was supported by JSPS Postdoctoral Fellowships for Research Abroad and J. S. by European Research Council Grant No. 267117.
[1] Ya. B. Zel'dovich and I. D. Novikov, Sov. Astron. **10**, 602 (1967).

[2] S. Hawking, Mon. Not. R. Astron. Soc. **152**, 75 (1971).

[3] B. J. Carr and S. W. Hawking, Mon. Not. R. Astron. Soc. **168**, 399 (1974).

[4] S. W. Hawking, Commun. Math. Phys. **43**, 199 (1975).

[5] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D **81**, 104019 (2010).

[6] G. F. Chapline, Nature (London) **253**, 251 (1975).

[7] A. Barrau, D. Blais, G. Boudoul, and D. Polarski, Ann. Phys. (Berlin) **13**, 115 (2004).

[8] S. Clesse and J. Garcia-Bellido, Phys. Rev. D **92**, 023524 (2015).

[9] B. J. Carr, F. Kuhnel, and M. Sandstat, Phys. Rev. D **94**, 044029 (2016).

[10] P. H. Frampton, Mod. Phys. Lett. A **31**, 1650093 (2016).

[11] R. H. Cyburt, B. D. Fields, and K. A. Olive, Phys. Lett. B **567**, 227 (2003).

[12] E. Di Valentino, E. Giusarma, M. Lattanzi, A. Melchiorri, and O. Mena, Phys. Rev. D **90**, 043534 (2014).

[13] C. G. Lacey and J. P. Ostriker, Astrophys. J. **299**, 633 (1985).

[14] R. Bean and J. Magueijo, Phys. Rev. D **66**, 063505 (2002).

[15] M. Kawasaki, A. Kusenko, and T. T. Yanagida, Phys. Lett. B **711**, 1 (2012).

[16] B. J. Carr and J. Silk, Primordial black holes as seeds for cosmic structure (unpublished).

[17] N. Afshordi, P. McDonald, and D. N. Spergel, Astrophys. J. Lett. **594**, L71 (2003).

[18] A. Kashlinksy, Astrophys. J. **823**, L25 (2016).

[19] M. Ricotti, J. P. Ostriker, and K. J. Mack, Astrophys. J. **680**, 829 (2008).

[20] Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D **95**, 043534 (2017).

[21] G. M. Fuller, A. Kusenko, and V. Takhistov, Phys. Rev. Lett. **119**, 061101 (2017).

[22] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. **116**, 201301 (2016).

[23] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. **117**, 061101 (2016).

[24] J. R. Bond and B. J. Carr, Mon. Not. R. Astron. Soc. **207**, 585 (1984).

[25] T. Kinugawa, K. Inayoshi, K. Hotokezaka, D. Nakauchi, and T. Nakamura, Mon. Not. R. Astron. Soc. **442**, 2963 (2014).

[26] B. J. Carr, Astron. Astrophys. **89**, 6 (1980).

[27] T. Nakamura, M. Sasaki, T. Tanaka, and K. S. Thorne, Astrophys. J. Lett. **487**, L139 (1997).

[28] A. M. Green, Phys. Rev. D **94**, 063530 (2016).

[29] F. Kuhnel and K. Freese, Phys. Rev. D **95**, 083508 (2017).

[30] B. J. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermäe, Phys. Rev. D **96**, 023514 (2017).

[31] B. J. Carr, Astrophys. J. **201**, 1 (1975).

[32] R. H. Cyburt and B. Davids, Phys. Rev. C **78**, 064614 (2008).

[33] D. P. Quinn, M. I. Wilkinson, M. J. Irwin, J. Marshall, A. Koch, and V. Belokurov, Mon. Not. R. Astron. Soc. Lett. **396**, L11 (2009).

[34] M. A. Monroy-Rodríguez and C. Allen, Astrophys. J. **790**, 159 (2014).

[35] T. Nakama, J. Silk, and M. Kamionkowski, Phys. Rev. D **95**, 043511 (2017).

[36] J. García-Bellido, M. Peloso, and C. Unal, J. Cosmol. Astropart. Phys. 09 (2017) 013.

[37] B. J. Carr, in *Inflating Horizons in Particle Astrophysics and Cosmology*, edited by H. Susuki et al. (Universal Academy Press, Tokyo, 2005), pp. 119–149.

[38] M. Y. Khlopov and A. G. Polnarev, Phys. Lett. **97B**, 383 (1980).

[39] K. Jedamzik, M. Lemoine, and J. Martin, J. Cosmol. Astropart. Phys. 09 (2010) 034.

[40] T. Harada, C. M. Yoo, K. Kohri, K. i. Nakao, and S. Jhingan, Astrophys. J. **833**, 61 (2016).

[41] B. J. Carr, T. Tenkanen, and V. Vaskonen, Phys. Rev. D **96**, 063507 (2017).

[42] P. Ivanov, P. Naselsky, and I. Novikov, Phys. Rev. D **50**, 7173 (1994).

[43] J. Yokoyama, Astron. Astrophys. **318**, 673 (1997).

[44] J. S. Bullock and J. R. Primack, Phys. Rev. D **55**, 7423 (1997).

[45] J. Yokoyama, Phys. Rev. D **58**, 083510 (1998).

[46] M. Kawasaki, A. Kusenko, Y. Tada, and T. T. Yanagida, Phys. Rev. D **94**, 083523 (2016).

[47] M. Kawasaki and Y. Tada, J. Cosmol. Astropart. Phys. 08 (2016) 041.

[48] P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida, J. Cosmol. Astropart. Phys. 04 (2010) 023.

[49] J. Garriga, A. Vilenkin, and J. Zhang, J. Cosmol. Astropart. Phys. 02 (2016) 064.

[50] H. Deng and A. Vilenkin, J. Cosmol. Astropart. Phys. 12 (2017) 044.

[51] A. D. Dolgov, Phys. At. Nucl. **80**, 987 (2017).

[52] J. Kormendy and D. Richstone, Annu. Rev. Astron. Astrophys. **33**, 581 (1995).

[53] F. Hoyle and J. V. Narlikar, Proc. R. Soc. A **290**, 177 (1966).

[54] P. Meszaros, Astron. Astrophys. **38**, 5 (1975).

[55] J. Silk, Astrophys. J. **151**, 459 (1968).

[56] B. J. Carr and J. E. Lidsey, Phys. Rev. D **48**, 543 (1993).

[57] J. D. Barrow and P. Coles, Mon. Not. R. Astron. Soc. **248**, 52 (1991).

[58] J. Chluba, J. Hamann, and S. P. Patil, Int. J. Mod. Phys. D **24**, 1530023 (2015).

[59] K. Kohri, T. Nakama, and T. Suyama, Phys. Rev. D **90**, 083514 (2014).

[60] T. Nakama, T. Suyama, and J. Yokoyama, Phys. Rev. D **94**, 103522 (2016).

[61] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D **95**, 123510 (2017).

[62] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, Astrophys. J. **473**, 576 (1996).

[63] M. H. Abitbol, J. Chluba, J. C. Hill, and B. R. Johnson, Mon. Not. R. Astron. Soc. **471**, 1126 (2017).

[64] T. Nakama, T. Suyama, and J. Yokoyama, Phys. Rev. Lett. **113**, 061302 (2014).

[65] H. Tahiro and N. Sugiyama, Phys. Rev. D **78**, 023004 (2008).
[66] T. Nakama, T. Harada, A. G. Polnarev, and J. Yokoyama, J. Cosmol. Astropart. Phys. 01 (2014) 037.

[67] T. Harada, C. M. Yoo, K. Kohri, and K. I. Nakao, Phys. Rev. D **96**, 083517 (2017).

[68] T. Harada, C. M. Yoo, T. Nakama, and Y. Koga, Phys. Rev. D **91**, 084057 (2015).

[69] S. Young, C. T. Byrnes, and M. Sasaki, J. Cosmol. Astro-part. Phys. 07 (2014) 045.

[70] J. Chluba, A. L. Erickcek, and I. Ben-Dayan, Astrophys. J. **758**, 76 (2012).

[71] P. A. R. Ade et al. (Planck Collaboration), Astron. Astrophys. **594**, A13 (2016).

[72] B. J. Carr and M. Sakellariadou, Astrophys. J. **516**, 195 (1999).

[73] C. T. Byrnes, E. J. Copeland, and A. M. Green, Phys. Rev. D **86**, 043512 (2012).

[74] N. J. McConnell and C.-P. Ma, Astrophys. J. **764**, 184 (2013).

[75] V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, Phys. Rev. D **96**, 083524 (2017).

[76] B. J. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veermae, Phys. Rev. D **96**, 023514 (2017).

[77] A. Soltan, Mon. Not. R. Astron. Soc. Lett. **200**, 115 (1982).