The power spectrum on small scales: Robust constraints and comparing PBH methodologies

We compare primordial black hole (PBH) constraints on the power spectrum and mass distributions using the traditional Press Schechter formalism, peaks theory, and a recently developed version of peaks theory relevant to PBHs. We show that, provided the PBH formation criteria and the power spectrum smoothing are treated consistently, the constraints only vary by $\sim$10\% between methods (a difference that will become increasingly important with better data). Our robust constraints from PBHs take into account the effects of critical collapse, the non-linear relation between $\zeta$ and $\delta$, and the shift from the PBH mass to the power spectrum peak scale. We show that these constraints are remarkably similar to the pulsar timing array (PTA) constraints impacting the black hole masses detected by the LIGO and Virgo, but that the $\mu$-distortion constraints rule out supermassive black hole (SMBH) formation and potentially even the much lighter mass range of $\sim$(1-100) $\mathrm{M}_\odot$ that LIGO/Virgo probes.

horizon entry during radiation domination, there is an approximate one-to-one relation between the scale at which the primordial power spectrum has a large amplitude peak and the mass of PBHs that form. See [6][7][8] for reviews.
In order for PBHs to form, the amplitude of the power spectrum must become orders of magnitude larger than the value of 2 × 10 −9 detected on large scales, e.g. via observations of the cosmic microwave background (CMB) [9]. Precisely how much larger it must become is a matter of active research, with significantly differing values being quoted in the literature, typically varying between O(10 −3 ) and O(10 −2 ) with values at the lower end quoted in e.g. [10,11]. O(10 −1 ) values have also been considered in e.g. [12]. Since the power spectrum amplitude is only logarithmically sensitive to the allowed energy density fraction of PBHs, this variation has little to do with the different PBH masses or constraints being considered and instead is primarily due to differences in the theoretical techniques being used to relate the power spectrum amplitude to the abundance of PBHs. Primordial non-Gaussianity also has an important impact on the required power spectrum amplitude, see e.g. [13][14][15][16][17][18][19], but we will not consider that issue further in this paper. However, we do include an accurate approximation for the significant correction arising due to the non-linear relation between the density contrast and curvature perturbation, the importance of which has only recently been quantified [20][21][22][23][24].
In this paper we make the first detailed study of how the PBH mass distribution differs when using Press Schechter or peaks theory as well as a recently developed treatment of peaks theory [25], which solves a problem for PBHs related to the cloud-in-cloud problem. When a PBH forms, the final mass depends on both the amplitude and scale of the perturbation from which it forms [26], and the new treatment of peaks theory ensures that the amplitude of peaks are evaluated at the correct scale, giving the correct mass distribution and abundance. We also consider the sensitivity to the choice of the window function. We show that, provided that all quantities are calculated in a self-consistent way -for example, the choice of window function must be reflected in the collapse threshold δ c -all techniques and window functions lead to quite consistent results whereby the uncertainty in the power spectrum amplitude is only of order 10%. This is a much smaller variation than [27] found even due to just the choice of the window function alone, consistent with the corrections accounted for in [28]. We also note that, throughout this paper, we assume a fixed value for the collapse threshold of primordial perturbations. In reality, the exact value of the collapse threshold depends on the specific shape of each individual perturbation (see e.g. [29,30] for a recent discussion), and neglecting this gives an additional uncertainty of order a few percent.
The uncertainty in the initial conditions required to generate a required number of PBHs has important implications for relating observations of PBHs to observations of the associated enhanced amplitude of the primordial perturbations. This can be done, for example, via the observation of a stochastic background of gravitational waves measurable by pulsar timing arrays (PTAs) which measure frequencies corresponding to a horizon scale which could have formed the black holes observed by LIGO and Virgo. In general, understanding how to map from a PBH abundance to a power spectrum constraint is important for our understanding of the initial conditions of the universe and the constraints on models of inflation.
In the next section we introduce the calculation of the PBH mass distribution. In section III we discuss how the result depends on the calculation technique and window function and we use these results to calculate robust constraints on the primordial power spectrum in section IV, in particular showing that the pulsar timing array constraints are not inconsistent with the formation of LIGO mass PBHs. We conclude in section V, and some technical details of the observational constraints and the non-linear mapping from the curvature perturbation to the density contrast are contained in the appendices.

II. OBTAINING THE PBH MASS DISTRIBUTION
The procedure for obtaining the mass distribution from the power spectrum is similar for all three methods considered, and is based on connecting the PBH abundance Ω PBH to β(R) (the fraction of mass in the universe that is contained within PBHs at their formation time) and then to the power spectrum. In every case, the PBH abundance is calculated from the mass fraction using where R is the horizon scale at the time the PBH is forming, R eq is the horizon scale at matterradiation equality and the ratio takes into account the relative growth of the PBH fraction during radiation domination. The abundance is then related to the PBH mass function f (m) through which satisfies the normalisation condition This can then be related to the mass distribution ψ(m) through which is a PDF and hence satisfies the normalisation condition dm ψ(m) = 1.
The relation between β(R) and the power spectrum then depends on the method used. In this paper, three methods are considered: a Press-Schechter-like calculation (PS), the traditional peaks theory method (TP) described in the classic BBKS paper [31], and a modified peaks theory derived by Young and Musso (YM) [25]. We note that other variations of peaks theory have been developed in [23,27,32,33]. In the Press-Schechter formalism, the mass fraction is related to a probability distribution in the compaction function C by where the compaction is a smoothed version of the density contrast δ (see eq. (C4)). The probability density function is given by and the mass ratio m/M H takes into account the effect of critical collapse. In traditional peaks theory, the mass fraction is related to the number density of peaks, n, through where the number density is a function of ν = C/σ 0 , given by [31] n(ν) = 1 The modified peaks theory developed in [25] also has β related to n in a similar way to eq. (8), but with a factor of R 4 rather than R 3 , i.e.
This is required to counteract an extra inverse spatial dimension in the number density, given by [25] n(ν) = 16 where γ 0,2 and α are related to the width parameters σ n (R). These width parameters relate the probability density (in Press-Schechter) or number density (in the peaks theories) to the power spectrum through the relation where P δ R (k) is the compaction power spectrum, related to the power spectrum for ζ through W (k, R) is a window function applied to the power spectrum. In this paper, two window functions are considered: a real-space top-hat 1 , given in Fourier-space by and a Gaussian window function modified by a factor of 2 in the exponent as suggested in [28], It should be noted that, in the case of the modified Gaussian window function, the compaction referred to by C above is not technically the compaction, but is rather a "compaction-like" function.
The compaction (or compaction-like function) is related to the PBH mass through the critical collapse equation, where K, C c , and γ are numerical factors that depend on the window function used to smooth the power spectrum, as well as the shape of the density perturbation [24,29,34]. The values  [28]. We will take the values stated in [25]: K = 4 and C c = 0.55 for the top-hat window function, and K = 10 and C c = 0.25 for the modified Gaussian window function. For both window functions we take γ = 0.36.
In this paper we will frequently consider a power spectrum with a lognormal peak, as a simple parametrisation of a peaked power spectrum with a position and width that can be easily tuned.
The form is which has been appropriately normalised such that the constraint on A becomes independent of ∆ in the limit of a narrow peak, and it matches the delta function power spectrum A δ(ln(k/k p )) in this limit. We show this later in table II. The integral of this power spectrum over ln k is A, independently of the value of ∆. The width ∆ is a free parameter, and we will normally choose two representative values for the width, ∆ = 0.3 as a narrow peak which results in a PBH mass distribution not very different from that due to a delta-function power spectrum, and ∆ = 1 as a broad peak which is roughly what one would expect if the inflaton field dynamics change over a time-scale of 1-efolding during inflation. We note that such a peak should not be extrapolated to values of k very different in magnitude from k p (and of course the power spectrum needs to match the quasi scale-invariant spectrum observed on CMB scales), but in practice we have checked that both the power spectrum constraints and the PBH mass distribution do not depend on the shape of the peak when sufficiently far from the peak position (where the power spectrum amplitude is significantly smaller than the peak value). We are therefore not concerned (for the values of ∆ we focus on) that a lognormal peak exhibits a growth steeper than k 4 on scales far from the peak, even though this is the approximate maximum growth rate of the power spectrum in canonical single-field inflation [37][38][39]. A steeper growth can be achieved in e.g. multifield inflation [40,41].
It is convenient to state the peak scale k p in terms of the horizon mass it corresponds to, using the relation [42] M H = 17 g * 10.75 where g * is the number of relativistic degrees of freedom. We define the horizon mass at the peak of the power spectrum as

A. Effect of the calculation method and window function
Constraints on the PBH abundance can be used to place constraints on the amplitude of the primordial power spectrum. If the black holes in the LIGO merger events are considered to be primordial in origin, a fit of the masses and number of events can be used to constrain the PBH mass distribution, and hence the power spectrum. Recent studies have shown that, in this case, f PBH would have to lie between 10 −2 and 10 −3 , and would be closer to the lower of these two values [43][44][45][46]. See however recent papers [47][48][49][50] discussing the effect of interactions between binary and single PBHs, which suggests that a much larger value for f PBH is possible provided that PBH It can be seen that, when being careful with the combination of the window function and the corresponding critical collapse values, all the amplitudes are of the same order. When changing either the method or the window function while keeping the other fixed, the difference in the required amplitude is 20%. The biggest difference when taking both the window function and the calculation method into account is ∼ 32%. We note that the maximum value of the power spectrum does not vary nearly as much when ∆ changes as suggested by table I due to our parametrisation of the power spectrum definition (17). Choosing a different normalisation by leaving out the division by ∆ would instead lead to a divergent value of the power spectrum amplitude in the limit ∆ → 0, instead of a value which matches the delta function power spectrum.
We can also examine the amount of variability in the shapes of the mass distribution generated with different methods/window functions. The effect of changing the method is shown in fig. 1  We have shown that the different calculation methods result in an O(10%) shift in the required power spectrum amplitude, and a small difference in the shape and position of the mass distribution.
We expect the BBKS peaks method (TP) to provide a more accurate result than the Press-Schechter (PS) case [51], and that the modified version (YM) be better than TP, since it is a direct extension.
Although the differences are small, they will become important in the future as experiments that can probe the PBH mass distribution become more accurate. For the remainder of this work, we will use the modified Gaussian window function in eq. (15) and the traditional peaks theory (TP) method. This allows comparisons between other works that use the TP method and the results in this paper, which can then be compared between the different methods based on the differences highlighted here.
B. Effect of the peak width ∆ As shown in section III A, the calculated mass distributions have a shift in the peak position which depends on the width of the power spectrum peak used. Additionally, we expect the width of the mass distribution to increase. We can demonstrate these effects by calculating the mass distributions for a range of values of ∆ between zero (i.e. a delta function peak) and two. The It is immediately apparent that, even for a delta function peak in the power spectrum, there is a minimum width in the mass distribution, associated with the critical collapse effect described in section II. It can also be seen that for very narrow peaks in the primordial power spectrum, the resulting mass distribution hardly varies until ∆ 0.1. Beyond that point, the shift of the peak and the increased width become apparent. This means that whilst a monochromatic mass spectrum is unrealistic, studying a mass distribution with the minimum width due to critical collapse and a delta function power spectrum may be a good approximation to a physically realisable PBH mass distribution. The increasing width is also obvious, and can be quantified by fitting a lognormal mass distribution (the shape expected for PBHs arising from a smooth, symmetric peak) to data generated from the curves, and comparing the widths of these lognormals. The lognormal mass distribution is given by where m c is the mean of the distribution and σ ψ is the width (note the subscript to avoid confusion with the σ n (R) parameters appearing in section II). The resulting lognormal parameters are shown in table II, and show that, as expected, the width of the calculated mass distribution increases with the peak width, as well as the amplitude required to keep f PBH fixed. This minimum width appears to be much larger than is required in order for PBH decay to result in a sufficiently rapid transition from an early matter dominated era (caused by low mass PBHs) to radiation domination to generate an observable stochastic background of gravitational waves [52]. A noteworthy point here is that the typical mass of a PBH is actually significantly larger than the horizon mass corresponding to the scale at which the power spectrum peaks, m c /M H,P > 1.
At first glance, this statement may seem to be in disagreement with previous works where the expected PBH mass has been shown to be smaller than the horizon mass at re-entry. Physically, this apparent discrepancy is due to the fact that, if there is a narrow peak in the ζ power spectrum at a scale k p , the resultant perturbations will, on average, have a significantly larger characteristic scale r m . In the calculation presented here, this manifests itself in the fact that the variance σ 2 0 (R) peaks at a larger value of R than that corresponding to the scale k p (as calculated in [34] for example). Thus, the final mass of PBHs is smaller than the horizon mass corresponding to r m , but larger than the horizon mass corresponding to k p . The important conclusion drawn from this is that constraints on the PBH abundance for a given mass of PBH correspond to constraints on the primordial power spectrum at a larger value of k than have previously been calculated. Now we have a clear picture of how the different method and window function choices affect the mass distribution ψ and the amplitude required to generate a fixed f PBH , we can calculate the constraints on the power spectrum from PBHs, being careful about the consistency of our window function and critical collapse choices. We show the procedure for obtaining these constraints, and the final constraint plots, in the next section.

IV. THE CONSTRAINTS ON THE POWER SPECTRUM
A. Relevant constraints and how they are calculated Whilst calculating the PBH abundance with different methods has a huge effect on the calculated abundance and mass distribution, we have shown that the resultant uncertainty in constraints on the power spectrum is relatively small. We will now consider how observational limits on the PBH abundance, as well as a swathe of other observational probes, constrain the amplitude of the primordial power spectrum. The key additional constraints on small scales come from cosmic µdistortions [53] and a stochastic background of gravitational waves, which could be generated with a large amplitude due to the non-linear coupling between the scalar and tensor perturbations around the time of horizon entry. The calculation of many of these constraints follows closely the procedure presented in [37], and we therefore relegate the details to appendix B. However, we describe the constraints from PBHs in detail here, and we also highlight that constraints from PTAs have been updated to use the improved analysis of the NANOGrav 11 year data set [54]. There are additional small-scale constraints on the power spectrum, including for example those from y-distortions [55,56], 21cm observations [57][58][59][60][61] and the non-detection of ultra-compact minihaloes [12,[62][63][64][65]. We do not display the former because the combination of CMB constraints and µ-distortion constraints are more competitive on commensurate scales, and we do not display either of the latter because they depend on the dark matter model. Big bang nucleosynthesis constraints are discussed in e.g. [66][67][68]. Translating this power spectrum to Ω GW h 2 with eq. (B4), we can then compare the predicted signal with PTA constraints from the NANOGrav 11 year data set. We choose this data set because the new analysis takes errors in the modelling of the solar system ephemeris into account. This can have a large effect on the constraints which will need to be factored into the previous NANOGrav 9 year constraints, as well as those from other arrays such as the European Pulsar Timing Array (EPTA) which have previously been used to constrain the primordial power spectrum with induced gravitational waves. Those constraints should now be revised upwards, but the analysis would need to be redone in each case to quantify by exactly how much. Since the NANOGrav data set has pulsar timing data for 11 years of observations, it does not extend to quite as large scales as does the EPTA data, which is from 18 years of observations. This means that our constraints do not span as wide a range of scales (and hence PBH masses) as previous constraints in the literature show, but the constraints we do show are more robust to errors in solar system ephemeris modelling. We also avoid confusion over different analyses from different data sets, and are able to use the free spectrum constraints on Ω GW h 2 consistently throughout.
These constraints (taken from the bottom panel of fig. 3 in [54]) are the 2-σ constraints derived as a function of frequency so as to represent the sensitivity to monochromatic signals. This means that we will construct our constraints based on finding the limiting amplitude of the lognormal power spectrum to which the NANOGrav constraints would be sensitive. One could do a more sophisticated analysis, taking into account the fact that confidence in a detection would become even stronger if there are also weaker detections of a given signal on larger or smaller frequencies than where the strongest detection would come. We choose to just show the 2-σ constraints for clarity. We convert from frequency to scale with k = 2πc/f and then find the minimum value of A for each k p . The minimum value of A for each k p is found by scanning over all values of k for which NANOGrav has sensitivity. We plot the results in figs. 5 and 6 for ∆ = 0.3 and ∆ = 1, where again to be clear, the constraint on P R at a given k represents the maximum amplitude A for a lognormal power spectrum centred at k = k p such that the induced second-order gravitational waves would not be in conflict with the PTA constraints from the NANOGrav 11 year data set.

C. Constraints from PBHs
Constraints on primordial black holes are normally presented in terms of either f PBH or the mass fraction β, so a method is required to relate these to the power spectrum amplitude. A relation between f PBH (or equivalently Ω PBH ) is complicated by the fact that the redshifting factor in eq. (1) means that the required amplitude to generate a fixed f PBH varies with the peak positions (as demonstrated in sec III A). In general, the best way to overcome this would be to produce a relation for A as a function of both f PBH and the relevant mass scale. However, this is computationally expensive, and so a simplified approach is necessary. We can find an approximation by relating the power spectrum amplitude to a parameter that does not vary with the peak position, which we achieve by modifying eq. (1), adjusting the redshift factor by introducing a new scale R * , such that If R * is chosen to be close enough to the peak scale in the power spectrum, then the relation between this quantity and the power spectrum amplitude will be independent of the peak position.
This quantity cannot be treated exactly as the abundance, because the abundance is calculated in the super-horizon regime before PBHs form, whereas this is at some later time, corresponding to when the horizon scale is R * . This quantity can be related to the constraints for PBHs using The relation between the power spectrum amplitude and Ω PBH * for all three methods is shown in fig. 4  To obtain constraints on the power spectrum, Ω PBH * must be related to constraints on either β or f PBH . The most recent PBH constraints are those on β stated in [8], which is a version of the mass fraction β with common parameters normalised out. These constraints are calculated assuming that all the PBHs form at the same time (or equivalently, the same scale R), but it is possible to relate the constraints to Ω PBH * , and hence determine the constraints on the amplitude for the calculation used throughout this paper, where PBHs form over a range of different scales.
We obtain this relation from eqs. (6) and (8) where the monochromatic PBH mass M in [8] has been substituted for the mean lognormal mass m c (the constraints do not change significantly when considering a reasonably narrow PBH mass distribution [69,70]). It can immediately be seen that, combining eqs. (24) and (25) Since solar mass PBHs are of special interest, it is sensible to rescale the mass fraction to be in terms of solar masses. The relation for this is We can then be relate this to the quantity Ω PBH * using eq. (23) to give For convenience, we have chosen R * such that the corresponding mass scale M * is approximately Substituting in the value of the horizon mass at matter-radiation equality, M eq = 2.8 × 10 17 M , the relation becomes Recent papers [20][21][22][23][24] have discussed the effect of the non-linear relation between the curvature perturbation ζ and the density contrast δ on the PBH abundance. The point is that, even if the level of primordial non-Gaussianity of ζ is taken to be zero, δ will not have a Gaussian distribution, and subsequently nor will the compaction. The non-linearity is difficult to account for, especially if window functions other than a top-hat are considered. This is discussed in some detail in appendix C, with the conclusion that constraints on the power spectrum will be approximately 1.98 times weaker once the non-linearity is included in the calculation. We include this factor in the PBH lines in figs. 5 and 6.
By applying the method described in this section, we are taking into account the effects of critical collapse (making sure it is treated consistently with the choice of window function), the shift between the PBH mass and the peak scale k p , and the non-linear relation between ζ and δ.
This is the first time that all of these effects have been captured simultaneously.

D. Summarising all the constraints
In fig. 5 we put together the key observational constraints to show the principal current constraints on the primordial power spectrum. The power spectrum has been accurately measured on large scales whilst PBHs constrain -albeit weakly -a far larger range of scales. We do not show PBH constraints on masses close to matter-radiation equality because we always assume PBHs form during radiation domination, and the smallest scale constrained corresponds to a PBH with m c ∼ 10 −24 M , which evaporates around the time of big bang nucleosynthesis.
By coincidence the PTA measurements constrain the power spectrum amplitude to almost the same amplitude as the non-detection of PBHs, meaning that there is a potential tension between the PTA bounds and any claim that LIGO detected PBHs (see fig. 5). This has been studied by various groups [37,[71][72][73][74][75][76][77][78][79], with no consensus reached on how severe the tension is. The impact of the PBH density profile was studied in depth in [24] but the PTA constraint was not varied to reflect changes in the shape of the primordial power spectrum. For example [75] claim that f PBH < 10 −6 over a significant range of PBH masses and the power spectrum constraint plots in [37] appear to show a significant tension. By making a careful study of the power spectrum amplitude required to generate PBHs, including the important reduction in the PBH constraining power due to the non-linear relation between ζ and δ, and using improved NANOGrav constraints, we have shown that there is no significant tension between generating LIGO mass PBHs and the PTA constraints.
We note that the slight overlap between the PBH and PTA constraint lines is not significant  [83]. It seems plausible that the associated stochastic background could be detectable with current PTA data if a dedicated search was made by using specific GW templates generated by power spectra that cause LIGO mass PBHs to form. The cosmic µ-distortion places an upper limit on the maximum PBH mass which can be generated by the collapse of large amplitude perturbations shortly after horizon reentry. The maximum mass decreases as the power spectrum width ∆ increases, but even for a narrow peak with ∆ = 0.3 the initial PBH mass cannot be much greater than 10 4 M , which is much smaller than the supermassive BHs seen in the centre of most galaxies even at high redshift, with masses 10 6 -10 9 M , whose origin remains a mystery. However, such large PBHs could still act as a seed to the SMBHs although one then needs to evade the strong Planck constraints on dark matter isocurvature modes [85,86]. For even broader power spectra the µ-distortion constraints rule out an ever greater range of PBH masses, and for ∆ = 2 they extend as far as the peak PTA constraint and thereby even rule out LIGO mass PBHs. Since such a wide peak in the primordial power spectrum provides the preferred PBH mass distribution width when fitting to LIGO data, it appears that the µ-distortions may surprisingly provide a stronger constraint on models in which all LIGO black holes are PBHs than the PTA constraints. Of course this conclusion may also depend on the assumed shape of the power spectrum peak.
Future constraints from µ-distortions and the gravitational wave background will significantly affect the PBH landscape. To examine the maximum extent of these future constraints, we calculate the PBH lines in the case that zero PBHs form in the observable universe. This is done using the method described in [87], particularly eq. (7) of that paper, but with β replaced with the Ω PBH * parameter used in this paper. For reasons summarised in [87], these extreme constraints might actually apply to the case of evaporated PBHs. Extremely tight constraints on f PBH for M PBH 10 −6 M are also possible if the majority of dark matter consists of "standard" WIMPs [4,[88][89][90][91][92].
We show these constraints in fig. 6, as well as future µ-distortion constraints from a detector like the Primordial Inflation Explorer (PIXIE) [93], and future gravitational wave background constraints from the Square Kilometre Array (SKA), the Laser Interferometer Space Antenna (LISA), and the Einstein Telescope (ET) 2 . The SKA constraints are derived from the sensitivity curve calculated in [94], the LISA constraints are derived from the most optimistic sensitivity curve in fig. 1 of [95], and the ET constraints are derived from fig. 13 of [96].
It can be seen that the SKA constraints are so tight that a non-detection will indicate that no PBHs can exist in the LIGO range of masses, and hence that the LIGO merger events cannot possibly be explained with a primordial origin. Additionally, the combined effect of the µ-distortion,

V. CONCLUSIONS
We have made the first detailed analysis of how the PBH mass distribution shape and amplitude varies between three different techniques to calculate the primordial mass distribution: Press Schechter, traditional peaks theory and a newly developed peaks theory variation. We also consider two choices of the window function, a real-space top-hat and a modified Gaussian. We show that the amplitude of the primordial power spectrum only varies by O(10%) for different choices, far smaller than may have been expected based on the large range of values of the power spectrum amplitude considered in the literature. A substantial variation remains depending on the shape of the peak in the primordial power spectrum, but this reflects a change in the physical theory rather than a change in methodology. The results are summarised in table I while fig. 1 shows that the mass distribution shape hardly changes depending on the calculation technique. These differences, while not significant now, will be important for future data that probes the PBH mass distribution accurately, at which point an improvement of the TP method, such as the Young-Musso technique, should be used. We also show that the PBH mass distribution becomes broader as the power spectrum peak becomes broader, as highlighted in fig. 3. In the limit of a narrow lognormal peak (∆ 0.3) the mass distribution tends to a constant width which is set by critical collapse, making a peak of this width a well-motivated choice.
We have also calculated robust constraints on the primordial power spectrum from PBHs, taking into account the effects of critical collapse and the non-linear relation between ζ and δ, as well as the choice of window function and the relation between the PBH mass scale and the peak power spectrum scale. This leads to tighter constraints that are shifted to different values of k compared to those presented in [8]. We show a summary of all of the key bounds on the amplitude of the primordial power spectrum in fig. 5. We stress that all the constraints must be recalculated when the shape of the primordial power spectrum peak is varied, and in the figure we choose ∆ = 0.3 as a representative narrow peak and ∆ = 1 as a broader peak. In both cases the PTA constraints (we use a recently improved data set from the NANOGrav collaboration) are almost identical to those from PBHs in the mass range that LIGO also probes. This interesting coincidence means that it is premature to rule out the possibility that LIGO detected PBHs that formed from large amplitude density perturbations during radiation domination, but if that is the case then there is a realistic hope that the PTA measurements will detect a stochastic background of gravitational waves in the near future and a dedicated analysis should be made. We note that the non-linear relation between ζ and δ weakens the PBH constraints by about a factor of 2, and had we not taken this into account (and normally it is not taken into account) we would have erroneously concluded that the PTA constraints do not come close to ruling out the formation of LIGO mass PBHs. However, we caution that if all BH binaries detected by LIGO were due to PBHs then the PBH mass distribution should be so broad (σ ψ 0.8 corresponding to ∆ = 2) that the cosmic µ-distortion constraints spread to relatively small masses and alternative shapes of the primordial power spectrum which are more "top-hat"-like than the lognormal power spectrum studied here should be considered.
In fig. 6 we show constraints on the primordial power spectrum that could be achieved in the foreseeable future (assuming there is no detection) from a PIXIE-like experiment measuring µ- Here we explain our procedure to produce constraints when using a real-space top-hat window function, which corresponds to a rapidly oscillating window function in Fourier-space, with consequent convergence issues. The width parameter σ 0 (R) is shown in fig. A.1 for a delta function peak (left) and the lognormal widths ∆ = 0.3 (middle) and ∆ = 1 (right). It can be seen that the oscillatory nature of the top-hat window function leads to a ringing effect in the width parameter σ 0 (R). For broader peaks in the power spectrum, this ringing effect merges into a constant height for large values of R. This leads to a divergent integral when evaluating eq. (1), and so the mass distribution cannot be calculated with this window function without some form of adjustment. It is common to supress the large-R constant effect using a transfer function, but this method is not compatible with other parts of our calculation (i.e. [28]). Therefore, we take an alternative approach, which is to adjust the calculation of σ n (R) in eq. (12) with a large-k cutoff. This is placed at the point where the window function reaches its first trough, which is at 4.49/R. This solves the divergence problem and removes the ringing/constant effect, but it must be noted that the window function is technically not a true top-hat any more.
Appendix B: Observational constraints 1. Constraints due to spectral distortions of the CMB Spectral distortions of the energy spectrum of the CMB are able to constrain the primordial power spectrum on small scales. They quantify deviations from the black-body temperature distribution of the CMB, caused by energy injection and removal from the plasma in the early universe.
A large boost in the primordial power spectrum at a particular scale or over a range of scales will lead to fluctuations in the density of the baryons and photons as a function of scale after reheating.
This will mean that the photon distributions on different scales will be described by different blackbodies, and as those photons mix via Thomson scattering, a spectral distortion will be induced if Compton scattering, Double Compton scattering and Bremsstrahlung processes aren't efficient enough to bring them into equilibrium. So-called y-distortions quantify late-time processes and place constraints on larger modes k < 3 Mpc −1 , whilst µ-distortions quantify earlier energy injection and removal and hence constrain the smaller scales, up to k ∼ 10 4 Mpc −1 which will be most interesting for PBH production. The final µ-distortions induced by the scalar perturbations can be approximated by [55] wherek = k/1 Mpc −1 and k min 1 Mpc −1 . Given a particular form for the power spectrum, this can be used to compute the total induced µ or y-distortion. Comparing this with observations then results in constraints on the primordial power spectrum.
The Far-InfraRed Absolute Spectrophotometer (FIRAS) instrument on board the COsmic Background Explorer (COBE) satellite measured spectral distortions to be smaller than ∆ρ γ /ρ γ < 6 × 10 −5 [97], and a proposed future detector such as the Primordial Inflation Explorer (PIXIE) [98], or a more recent proposal [99] aims for constraints of ∆ρ γ /ρ γ < 8 × 10 −9 . To calculate the constraints on the amplitude of the power spectrum due to the COBE/FIRAS observations, we insert eq. (17) into eq. (B1) and set µ = 9×10 −5 which is the 2-σ constraint. We can then rearrange for A and compute the integral over k, plotting the constraint on A for each k p . Our results for lognormal power spectra of widths ∆ = 0.3 and ∆ = 1 are shown in fig. 5. For complete clarity, the constraint on P R at a given k represents the maximum amplitude A for a lognormal power spectrum centred at k = k p so as not to induce µ-distortions that would be in conflict with the COBE/FIRAS constraint of µ < 9 × 10 −5 .

The stochastic gravitational wave background
Here we summarise how the GW background can be calculated given a primordial power spectrum, adding more details to section IV B. The contribution to the tensor power spectrum from the square of the scalar power spectrum is given by [100] where u = |k−k|/k, v =k/k andk is the wavelength corresponding to the scalar source. I(v, u, kτ ) is a highly oscillatory function which contains the source information. We solve this integral numerically but note that it can be solved analytically in some regimes [101]. The observational quantity related to this power spectrum is the energy density of gravitational waves given by If we assume that the entire contribution to any stochastic background detection is from the tensor power spectrum in eq. (B3), then constraints on the stochastic background can be translated to constraints on the scalar power spectrum. This is a conservative constraint, as there may be other unresolved astrophysical contributions to the signal. If a detection is made, as opposed to an upper limit on the amplitude from non-detection, spectral information of the signal will be required to distinguish between the possible sources. To calculate the constraints on the primordial power spectrum, we first calculate Ω GW h 2 today as a function of k by inserting the lognormal power spectrum in eq. (17) with given k p and ∆ into eq. (B3), pulling out the amplitude A which is the quantity that we aim to constrain. We perform this integral numerically once for each value of ∆, and the results can be shifted post-integration for any value of k p .
Appendix C: The non-linear relationship between ζ and δ In recent years, there has been a large amount of literature discussing the fact that, even if the curvature perturbation ζ is Gaussian, the density contrast will not be [20][21][22][23][24], due to the non-linear relationship between the 2 parameters. In the super-horizon limit, the relationship between the 2 parameters can be calculated with a gradient-expansion approach. At first order in gradients, the full non-linear relationship, in polar coordinates and assuming spherical symmetry, is given by whilst the linear relation is For simplicity, we will set the equation-of-state parameter w = 1/3 from here on.
We can define a time-independent component of the density contrast, where R is taken to be the scale of the perturbation. The compaction function C(x, R) is obtained by calculating the mass excess δM within a sphere of radius R, and dividing by R, which corresponds to smoothing the time-independent component of the density contrast with a top-hat smoothing function, C(x, R) = δM R = d 3 y δ TI (x − y)W (y, R).
Performing this integral gives an expression for the compaction function at the centre of spherically symmetric peaks: where C L is the expression one would obtain using the linear relation above, where the prime denotes a derivative with respect to the smoothing scale R.
The rare, large-amplitude peaks from which PBHs form are well approximated by sphericallysymmetric peaks [31], and so the above equation can be used to relate relevant peaks in C L to peaks in the compaction C. We note that the compaction has a maximum value, C max = 2/3, corresponding to C L = 4/3. For higher values of C L , the compaction decreases -and perturbations of this type correspond to a case for which PBH formation has not been simulated. For this reason, only perturbations with C L < 4/3 are typically considered -although in practice this has little effect on the PBH abundance since such large values of C L are exponentially suppressed.
If we then wish to calculate parameters related to the PBH abundance, we can simply replace the equation for the PBH mass, eq. (16), with a corresponding equation which relates the PBH mass to the linear, Gaussian component of the compaction instead In order to make an analytic estimate for how constraints on the power spectrum are affected by this non-linearity, we can make a simple assumption that all peaks which form PBHs are close to the critical amplitude (since the abundance of significantly larger peaks is exponentially suppressed).
In this simple case, and assuming C c = 0.55 (the case for the top-hat window function, see eq. (14)), the critical amplitude for the linear component of the compaction is C c,L ≈ 0.77, i.e. we can assume that peaks in the linear field need to have an amplitude 1.41 times larger than if we assumed a linear relation between ζ and δ, as in eq. (C2). Therefore, the power spectrum (which is proportional to the variance of perturbations) should be approximately 1.41 2 = 1.98 times greater. We can test this approximation by comparing the full calculation of the amplitude required to generate a fixed abundance f PBH = 2×10 −3 in the linear and non-linear cases. For the two lognormal power spectra considered in this paper, with widths ∆ = 1 and 0.3, the approximation holds to the precision of two decimal places stated above. Although this validity may vary with the position of the peak, we assume it holds globally for the results shown in figs. 5 and 6.
For the top-hat window function, there is a relatively simple analytic relationship relating the compaction function to the curvature perturbation (which we assume to be Gaussian). However, we note that if one instead uses a Gaussian window function, as we have considered in this paper, there