Primordial black hole isocurvature modes from non-Gaussianity

Primordial black holes (PBHs) are black holes that might have formed in high density regions in the early universe. The presence of local-type non-Gaussianity can lead to large-scale fluctuations in the PBH formation rate. If PBHs make up a non-negligible fraction of dark matter, these fluctuations can appear as isocurvature modes, and be used to constrain the amplitude of non-Gaussianity. Assuming that the parameters of non-Gaussianity are constant over all scales, we build upon the results of previous work by extending the calculation to include peaks theory and making use of the compaction $C$ for the formation criteria, accounting for non-linearities between $C$ and the curvature perturbation $\zeta$. For quadratic models of non-Gaussianity, our updated calculation gives constraints that are largely unaltered compared to those previously found, while for cubic models the constraints worsen significantly. In case all of the DM is made up of PBHs, the parameters of non-Gaussianity are $-2.9\cdot10^{-4}


Introduction
Primordial Black Holes (PBHs) are black holes that could have formed in the early universe in the radiation dominated regime after inflation.PBHs are capable of providing an explanation for numerous observed and unexplained cosmological phenomena [1].Most relevantly, PBHs could (fully or partially) make up dark matter (DM) (see e.g.ref. [2] for an overview), be the seeds for supermassive black holes in the centres of galaxies [3] and be the source of observed LIGO/Virgo events [4].
Regions where the density exceeds a critical value as the density perturbation enters the horizon collapse and form a PBH.This threshold value was first studied by Carr and Hawking [29] who considered that gravity dominates over pressure when the length scale of the density fluctuation exceeds the Jeans length.The jeans length R J and the horizon length R H are related by R J /R H = √ w, with w the equation of state parameter relating the pressure p and density ρ by p = wρc 2 .It takes the value w = 1/3 in the radiation dominated regime.It was found in ref. [30] that for the region to collapse, the density contrast δ ≡ δρ/ρ, with δρ ≡ ρ − ρ b the difference between the density and the background density ρ b , should exceed a critical value of δ c ≈ w at horizon crossing.
Many papers have since considered PBH formation (e.g.[31][32][33][34][35][36][37][38]), making use of both numerical and analytical methods, finding that the critical value is closer to δ c ≈ 0.5 (where the exact value depends on the specific profile shape of the perturbation), and finding also that the PBH mass follows a scaling law given by with κ ≈ 4 and γ ≈ 0.36 (where, again, the exact values can vary with the profile shape).
The PBH mass therefore depends on the scale of the perturbation (described by the horizon mass, M H ) and amplitude δ of the perturbation which forms a PBH.The PBH mass is therefore of the same order as the horizon mass (see e.g.ref. [39]) where we have substituted the Hubble horizon mass in the radiation dominated regime in terms of cosmic time t, the speed of light c and Newton's gravitational constant G. Thus, depending on the time at which the PBH forms, the mass can be arbitrarily small or large.
The abundance of PBHs is described by the parameter with ρ PBH the PBH density and ρ the density of the universe.Both are evaluated at the time of PBH formation t f .Cosmological observations place constraints on the abundance of PBHs (see e.g.refs.[2,39] and references therein).These constraints usually give an upper bound on the fraction of DM that is made up of PBHs of a single mass where Ω PBH ≡ ρ PBH /ρ and Ω DM ≡ ρ DM /ρ the density parameters of PBHs and DM, respectively.
The abundance of PBHs can be strongly affected by primordial non-Gaussianity in the density distribution [40][41][42][43][44].It is therefore important to take non-Gaussianity into account when calculating PBH abundances, but we can also use this tight dependency to find constraints on the parameters that describe non-Gaussianity.Previous work [45] has done this by studying the effects of a peak-background split on modal coupling and constraints on the iso-curvature perturbations of the Planck Collaboration [46], which gave constraints on the relation between the fraction f PBH and the parameters of non-Gaussianity.Integral to this calculation was the use of Press-Schechter theory to calculate the abundance of PBHs that form in regions where the density perturbation exceeds the critical value.
In this paper, we will take into account recent developments in the field to calculate constraints more accurately (as described in e.g.ref. [47]).These developments suggest the use of peaks theory and the compaction (instead of Press-Schechter theory and the density contrast), as well as accounting for the mass-scaling relationship of PBHs.
This paper is structured as follows: in section 2 we give a brief overview of non-Gaussianity and the calculation of PBH abundance.In section 3 we describe how the abundance of PBHs can be evaluated.In section 4 we derive up to date constraints on the abundance of PBHs, before summarising our findings in section 5.

Non-Gaussianity
Primordial fluctuations created during inflation are predicted and observed to follow a distribution that is very close to Gaussian on the CMB scales [48].However, many models of inflation predict the emergence of some levels of non-Gaussianity (see e.g.ref. [49]).The detection of Non-Gaussianity therefore offers important insights into the workings of inflation.The presence of even small amounts of primordial non-Gaussianity can also have a large impact on the amount of PBHs that are formed [40][41][42][43][44].
The curvature perturbation ζ appears as a perturbative quantity in the FLRW metric with a the scale factor.Local-type non-Gaussianity can be modeled by expanding the curvature perturbation as a polynomial of a Gaussian distributed curvature perturbation where f = 3f NL,local /5 and g = 9g NL,local /25 describe the amplitude of non-Gaussianity.In eq. ( 2.2) the variance σ 2 = ζ 2 G is subtracted to ensure that the average value ⟨ζ⟩ vanishes.Because the abundance of PBHs is sensitive to non-Gaussianity, it is important to take non-Gaussianity into account when calculating PBH abundances.However, as we will see later, this can also result in strong constraints on the non-Gaussianity parameters if PBHs make up a non-negligible fraction of dark matter.

Primordial black hole formation and the effect of non-Gaussianity
PBHs do not form at a single time and the PBH density parameter can be expressed as [37,50] with M min and M max the smallest and largest horizon masses at which PBHs form and M eq the horizon mass at matter-radiation equality.Note that this equation assumes pure radiation domination right up until matter-radiation equality.The term (M eq /M H ) 1/2 ∝ 1/a H accounts for the redshift of the PBH density parameter during radiation domination, where a H is the scale factor at the time of PBH formation.With the horizon mass serving as a parameter of time, this integral integrates the PBH abundance over the whole period in which PBHs form.
Assuming that PBHs form over a short time-interval, parameterised by a single horizon mass, a good estimate of the PBH density parameter can be found as This assumption can be valid in the case that the power spectrum peaks sharply at this scale (which will be assumed throughout this paper) -meaning that PBH formation at other scales is negligible.This also coincides with the condition that the local-type expansion, equation 2.2, gives a valid description of the statistics of the compaction (see ref. [43] for more discussion).We use M eq = 2.8 • 10 17 M ⊙ [37,51] and choose M H = M ⊙ throughout this paper, and note that the final constraints on the non-Gaussianity parameters depend only very weakly on this choice (see also the discussion in ref. [45]), justifying also the approximation done in eq.(2.4).Using the definition of the fraction f PBH , eq. (1.4), and Ω CDM = 0.26 [52], we then find Ω CDM f PBH ∼ 10 −10 f PBH . (2.5) In addition to the already mentioned effect on the total abundance of PBHs [40][41][42][43][44], non-Gaussianity can also lead to the production of isocurvature modes in the early universe.The modal coupling which can arise as a result of the non-Gaussianity means that the amplitude of small-scale perturbations can be coupled to a long wavelength perturbation.This means that PBH formation can be enhanced (reduced) in regions where the small-scale perturbations are larger (smaller).These fluctuations in the PBH formation rate therefore appear as dark matter isocurvature perturbations.Refs.[45,53] calculated the amplitude of these isocurvature perturbations, and used costraints from the Planck Collaboration [46] on isocurvature modes to find constraints on the non-Gaussianity parameters (as a function of the PBH abundance).The authors used a Press-Schechter approach, and applied the statistics of the curvature perturbation, derived from eq. (2.2), (see e.g.ref. [41]).
It has since been argued that using peaks theory is more suitable for calculating the PBH abundance [54][55][56].As with Press-Schechter theory, peaks theory assumes that PBHs form in regions where perturbations exceed a critical value, but it also introduces the condition that PBHs form in regions where the perturbations are at a maximum.Additionally, instead of using the curvature perturbation or density contrast to describe when a region forms a PBH, refs.[57,58] argue that the compaction is a more appropriate parameter use.In the above, M (x, r, t) is the Misner-Sharp mass, with M b (x, r, t) its background value.The Misner-Sharp mass gives the mass within a sphere of areal radius R(x, r, t) = a(t) exp(ζ(x))r with spherical coordinate radius r, centred around position x and evaluated at time t.Refs.[43,47] have recently considered the effect of localtype non-Gaussianity on the statistics of the compaction, and we will apply their methods here.
As with the density contrast, PBHs form in regions where the compaction exceed a critical value.Furthermore, like the density contrast, the compaction directly measures the overabundance of mass in a region and is therefore better suited than the curvature perturbation for determining when a region collapses.The compaction is time-independent on super-horizon scales, and can also be expressed as the time-independent component of the top-hat smoothed density contrast δ TH [58] where ϵ = r/r H , the ratio between the perturbation scale r and the Hubble scale r H . See ref. [58] for more discussion on the compaction and its use as the formation criterion.
In light of these developments, ref. [47] has evaluated the abundance of PBHs in the presence of local-type non-Gaussianity using peaks theory and the compaction, in addition to using the correct mass scaling in eq.(1.1).In this work, we will use these same methods to calculate the PBH abundance and use these methods to update the constraints on non-Gaussianity parameters previously found in ref. [45] from the isocurvature modes.For clarity, we will refer to the results of ref. [45] as the results from Young and Byrnes.

Compaction and peaks theory
In this section we will discuss how using the compaction, eq.(2.6), can be used to determine the abundance of PBHs.We will follow the discussion in ref. [47] in this section.

Compaction
The compaction has a form similar to the density contrast.Indeed, the compaction can be written in terms of the density contrast as where the integral of the density contrast over volume evaluates precisely the mass difference in eq.(2.6).On super-horizon scales, the density contrast δ is related to the curvature perturbation in real space by the non-linear relation We shall take the equation of state parameter to be w = 1/3 for the radiation dominated regime, and will take the high peak limit throughout.The high-peak limit assumes that PBHs form in the high positive tail of the density distribution (or correspondingly, that peaks in the density fluctuation must have a high amplitude for PBH formation).In the high peak limit, perturbations can be approximated as spherically symmetric [54], and we find where ζ ′ ≡ dζ/dr.Eq. (3.1) then becomes having used R = a exp(ζ)r.Taking the linear component C 1 = −4rζ ′ (r)/3, this gives the relation This quadratic equation is sketched in figure 2 and has a maximum at C 1 = 4/3 ≡ C 1,to , which we define as a turnover point, denoted by the subscript 'to'.The corresponding value of C is C to ≡ C(C 1,to ) = 2/3.Perturbations with C 1 < C 1,to are referred to as type I perturbations, while perturbations with C 1 > C 1,to are referred to as type II perturbations.Type II perturbations are exponentially suppressed compared to type I perturbations as we will later see in eq.(3.35) and section 3.3, and the formation mechanism of PBHs from type II perturbations is not well understood [47,59].We will therefore limit our discussion to PBHs formed from type I perturbations.

The effect of non-Gaussianity on the compaction
The effect of non-Gaussianity on the compaction can be described by consulting eq.(2.2) to find that in the presence of non-Gaussianity, the linear term of the compaction becomes The −4rζ ′ G (r)/3 term arises from the smoothing of the density contrast over a top-hat smoothing function, whilst the ζ G (r) term arises from the surface term, corresponding to a surface smoothing function.The top-hat and surface smoothing functions are defined, respectively, as with θ H (x) the Heaviside step function and δ D (x) the Dirac delta function.Their respective Fourier transforms are We can use the smoothing functions to construct the following quantities which appear in the compaction, in eq.(3.6): In each case, we have made the assumption of spherical symmetry for the second equality.We also define the following correlation functions; with To describe the PDF of C G and ζ r , we define the following quantities, which are normalised to their variance, ) The PDF is then given by the two-variate Gaussian distribution with Y = (ν, ν r ) and Σ the covariance matrix.To diagonalise this PDF, we introduce the variable where is the correlation function of ν and ν r .Like ν and ν r , the parameter z r follows a Gaussian distribution This gives the diagonalised PDF The high peak limit implies γ 0r ν ≫ 1 − γ 2 0r .This implies that the Gaussian distribution in eq.(3.21) has a mean that is much greater than its variance.Such distributions can be estimated as Dirac delta functions (see figure 1), and we therefore approximate Integrating the PDF over ν r therefore results in the substitution ν r = γ 0r ν, or equivalently For brevity we define γ ≡ γ 0r σ r /σ 0 .From eq. (3.6) we then find where In this paper, to ensure that the local-type expansion considered is valid, we will consider only narrowly peaked power spectra, which can be well approximated by a Dirac delta function with A the amplitude.Although such a power law is unphysical, it justifies the use of the high peak limit and thus also spherical symmetry and is an accurate approximation of the lognormal power law often used in the literature [47,58,60].Using this expression of the power law and eqs.(3.13) and (3.15) we can calculate having used r = 2.74/k p , which relates the peak of the power spectrum to the amplitude scale (and maximises σ 0 ) [58].Using eq.(3.20), we then also find which is actually independent of σ r .

Primordial black hole abundance in peaks theory
Young and Byrnes used Press-Schechter theory to calculate the abundance of PBHs, which assumes that PBHs form at points where the perturbation exceeds a critical value.Peaks theory expends on this condition by adding the requirement that PBHs form at the point where the perturbation is at a local maximum.
The number density of peaks can be found to be [37] PBHs form in regions where the compaction exceeds a critical value C c .The abundance of PBHs can be found as We will use the mass scaling in eq.(1.1), which in terms of the compaction reads with K = 4, C c = 0.5 and γ = 0.36 [61].Cf. eq.(1.1), C c is the critical value of PBH formation, i.e.PBHs form in regions where C > C c .Combining these terms together with σ 1 /σ 0 = k p for the Dirac delta power spectrum (3.28) in eq.(3.13) gives1 taking only solutions where C 1 < C 1,to , corresponding to type I perturbations.For C = C c = 0.5, this gives C 1 = 0.67 ≡ C 1,c . 2 This means that PBHs form when This corresponds to the range for C 1 where the compaction exceeds the critical value, but does not exceed the turnover point where perturbations become type II perturbations (see figure 2).Because eq.(3.35) contains an integral over C G we need to invert eq.(3.25) to find the range of C G that corresponds with the above range of C 1 .We will do this for the quadratic and cubic models separately.

Cubic expansion
For a cubic expansion (i.e.f = 0), eq.(3.25) becomes which we can invert to find three solutions with Similar to before, the range of integration depends on the value of g.The general behaviour again depends on the square root term.However, because Q i could also be complex, we no longer necessarily require 4 g 3 + 27C 2 1 g 4 > 0. As before, we indicate all ranges of C G with Greek letters and include them in figure 4.
4 Constraints on the primordial black hole abundance

Peak-background split
Under a peak-background split, perturbations could be split into a small-scale "peak" component ζ s and a large-scale "background" component Up to first order in Fourier space, eq.(3.2) becomes Because large-scale fluctuations are suppressed by a factor k 2 , long wavelength perturbations do not contribute directly to PBH formation [62].In the presence of local-type non-Gaussianity, they do however contribute to PBH formation indirectly.The amplitude of small wavelengths will be boosted around peaks of long wavelengths, increasing the probability that the total perturbation exceeds the critical perturbation [45] (see figure 5).Small-scale perturbations need to be much larger than those with CMB scale wavelengths in order for a significant number of PBHs to be formed.Following the approach in Young and Byrnes, we will assume ζ l ≪ 1 and analyse the effect that this has on the abundance of PBHs under quadratic and cubic expansion of non-Gaussianity.The effect of modal coupling.Normally, the amplitude of a short wavelength mode would not exceed the critical value of PBH formation.However, when a long wavelength mode is present in the background, the amplitude of short wavelength modes could be boosted to above the critical value.In this figure, we have subtracted terms that depend exclusively on long wavelength modes, similar to eq. ( 4.3) and (4.9).

Effects of peak-background split on the compaction
We will now consider the effects of the peak background split, eq.(4.1), on the compaction.

Quadratic expansion
Using eqs.(3.6) and (2.2) for the linear component of the compaction for a quadratic model (i.e.g = 0), we can add in the peak-background split using eq.(4.1) where in the final line we have neglected derivatives in ζ l , substituted C G = −4rζ ′ s /3 and we have used eq.(3.24), which implies that up to first order Eq. ( 4.3) can now be solved for C G to give For ζ l ≪ 1 the range of PBH formation for C G will be very similar to those calculated in the previous section: PBHs form in the ranges PBHs form in the range PBHs do not form.

Cubic expansion
When including only the linear and cubic terms in eqs.(3.6) and (2.2) (i.e.f = 0), the linear component of the compaction transforms as having again neglected derivatives in ζ l , substituted C G = −4rζ ′ s /3 and used eq.( 3.24) to find Similar to before, we can now solve eq. ( 4.9) for C G to find the range of integration for C G corresponding to C 1,c < C 1 < C 1,to .Doing so gives three solutions C G,i (C 1 ), with i ∈ {1, 2, 3} similar to those discussed in section 3.3.2.However, even though the additional quadratic term in the above equation does not bring any new free parameters, it still brings significant complexity to the algebraic form of these solutions.This makes it significantly harder to identify the range of integration.However, since we can assume ζ l ≪ 1, we can simplify the analysis and ignore two of the three solutions.The reason for this depends on the sign of g: with the other two solutions C G,k̸ =j (C 1 ) being complex.
• For g < 0 there will always be 3 real solutions C G,i (C 1 ), as long as C 1,c < C 1 < C 1,to and g is sufficiently small in magnitude (g ≳ −1).However, the smallest of these solutions C G,j (C 1 ) is smaller than the other two C G,k̸ =j (C 1 ) by a factor ∼ few.Because the integrand in eq.(3.35) depends exponentially on |C G |, contributions from larger values of |C G | will be exponentially suppressed.We can therefore ignore these larger solutions.
In conclusion, irrespective of the sign of g, we can then integrate over the range

Bias factor and non-Gaussianity parameters
Having properly found the bounds of the integral in eq. ( 3.35), we can calculate the PBH abundance β.Eqs.(4.3) and (4.9) offer perturbations compared to eq. (3.6) in terms of ζ l .We can use this to express the relative change in PBH abundance under the peak-background split as Figure 6 shows δ β as a function of ζ l , and we see that, for the range of values considered, we can express δ β as a linear function of ζ l δ β = bζ l , (4.12) with b the bias factor [53].
The DM density parameter Ω DM up to first order in ζ is related to the background density parameter Ω DM by As the universe expands over time, matter perturbations evolve by a factor (a exp(ζ)) −3 , which means that the 3ζ term is the adiabatic mode to first order in ζ.The f PBH bζ term is an isocurvature mode and forms a deviation from the adiabatic mode.The isocurvature mode is either fully correlated or anti-correlated, depending on the sign of the non-Gaussianity parameters [45].The Planck Collaboration [46] has found constraints on isocurvature modes.On CMB scales these are with fully correlated modes corresponding to b > 0 and fully anti-correlated modes corresponding to b < 0. We can express β iso in the above as where P iso is the isocurvature power spectrum and P ζ the perturbation power spectrum.These are related as [45] If a fraction f PBH of DM is made up of PBHs, the constraints (4.14) can then be used to constrain b as −0.028 < bf PBH < 0.036.( With b being independent on ζ l , there are three free parameters in this inequality; f or g for respectively quadratic and cubic expansion, σ 0 , as it appears in eq.(3.35), and f PBH .Together with eq.(2.5), we can then find a relation between f PBH and f or g by the system of equations bf The sign of b is the same as that of f or g.Thus, the condition of the sign of b in eq. ( 4. 19), can also be understood as a condition on the sign of f or g.In eqs.(4.18) and (4.19), both β and b depend on σ 0 and f or g.Thus, eq. ( 4.18) can be solved for σ 0 as a function of f or g and f PBH .We can then insert this relation in eq.(4.19) for the dependency of b on σ 0 and solve to find a relation between f PBH and f or g.This relation gives constraints on f or g for a given value of f PBH .
The resulting relation between f PBH and f or g is shown in figure 7. We compare these to the constraints that were found in Young and Byrnes.To plot the results from Young and Byrnes, we followed the same approach as the one presented in their paper.
From the results in figure 7 we find that if we set f PBH = 1, which corresponds to all of the DM being made up of PBHs, we find the bounds for the quadratic and cubic models respectively.The results for the quadratic model are similar to those found by Young and Byrnes, but broader for the cubic model by order O (10).
We see that the constraints on f found in this work remain mostly unaltered in both shape and magnitude compared to those found by Young and Byrnes.For a cubic expansion on the other hand the constraints look significantly different than before in two different ways.Firstly, our results display significant broadening at all ranges of g, leading to much weaker constraints, and secondly, the constraints become a lot less symmetric.Both of these effects are explained below.with that of C 1 , given by eq.(4.9) We see that in eq.(4.23) the linear term of C 1 is perturbed by O(ζ 2 l ), the quadratic term by O(ζ l ), while the cubic term remains unperturbed.Because b is determined from perturbations in ζ l ≪ 1, we can conclude that the quadratic term in eq.(4.23) has the most dominant effect on b, and therefore also on f PBH by eq.(4.19).The same could be said for the quadratic term in eq.(4.22).However, while the prefactor of the quadratic term in eq.(4.22) is 3gζ l , that of the quadratic term in eq.(4.23) is 6gγζ l and is thus smaller by a factor 2γ = O(0.1).
Because the constraints in our work relied on perturbations in C 1 , while those found in Young and Byrnes relied on perturbations in ζ, this leads to a broadening of roughly magnitude O(10) for the constraints found in this work as we need g to be larger by the same order O(10) to compensate for this smaller prefactor.From figure 7 we can indeed see that for lower values of g, the constraints indeed worsen by a factor ∼ 10.For larger values of g non-linearity starts to have a more pronounced effect and constraints worsen even more.This effect does not occur for the quadratic expansion.The reason for that is that in a peak-background split, ζ transformed as [45] and C 1 as eq.( 4.3) In both cases, the perturbation lies only in the linear term with identical prefactor 2f ζ l and suppression by a factor γ is absent here.Therefore we expect that the constraints on f PBH found in this work would be similar to those found by Young and Byrnes.

Symmetrisation of constraints for cubic models
To explain why the constraints for cubic expansions are more symmetric for positive and negative values of g, we instead need to consider that the abundance of PBHs was found by integration.In this work, this integration was performed over C G , while in the work of Young and Byrnes this was done over ζ s .The range of integration in this work was determined from C 1 = C 1,c and C 1 = C 1,to , while that of Young and Byrnes was determined from ζ = ζ c .This is the critical value of ζ above which PBHs form, equivalent to the critical value of the density perturbation δ c .It was found in ref. [62,63], that this takes the value ζ c = 1.Because ζ l ≪ 1, we can ignore all terms that depend on ζ l in order to find the solutions for ζ s and C G of eqs.
Although these equations have the same form, there is again suppression in the equation of C 1 since g is smaller than g by a factor 3γ 2 = O(0.01).Therefore, for |g| ≲ 1, the solution to C G will be dominated by the linear term of eq.(4.27), with the cubic term offering only a slight perturbation.This means that for a given magnitude |g|, solutions for C G for g > 0 will be very close to those for g < 0, which symmetrises the constraints.The main difference between g > 0 and g < 0 in the constraints then comes from the different magnitude of the bounds of b in eq.(4.17) for positive and negative values of b, or equivalently positive and negative values of g.Meanwhile, the cubic term of eq. ( 4.26) only depends on g.Thus, for |g| ≲ 1 the cubic term could have a more appreciable effect on the solutions C G , and solutions for g > 0 could deviate significantly from those for g < 0, which asymmetrises the constraints.Given that our work relies on perturbations in C 1 , while that of Young and Byrnes relies on perturbations in ζ, we find that our constraints on cubic expansions is more symmetric between positive and negative values of g.
The same symmetrisation can actually also be observed in the quadratic model in figure 7, albeit with a lesser effect.Solutions to ζ and C 1 for quadratic models can roughly be found by solving the equations Similar to before, f is smaller than f by a factor 2γ = O(0.1),which means that the quadratic term has a more appreciable effect on the solutions to ζ s than it has on C G , which makes our results more symmetrical compared to those of Young and Byrnes.

Constraints for higher order terms
With the knowledge of the broadening of the cubic expansion we are now also in the position to understand what the effect of adding higher order terms to the non-Gaussian expansion of C 1 would be.Following the patterns of eq.(2.2) and (3.6) we find that ζ and C 1 transform for an n-th order term as with h the prefactor of the n-th order non-Gaussianity term and we ahve subtracted ⟨ζ n G ⟩ to ensure that ⟨ζ⟩ = 0. Expanding ζ G = ζ s + ζ l in short and long wavelength modes, Taylor expanding the resulting equations for small ζ l and ignoring again any terms that depend exclusively on ζ l gives for these two equations ) where h ′ = nhγ n−1 .In the first of these equations we have subtracted ζ n−1

G
to again ensure that ⟨ζ⟩ = 0, in case n is even.
The bias factor b is determined by the leading order term in ζ l .Having ignored all higher order terms in ζ l , we see that due to the factor γ n−2 , the second term in eq. ( 4.33) has an increasingly less dominating effect for orders n.As with the cubic term, the term γ n−2 ensures that higher order terms will be subject to more intense broadening as larger values of h are needed to produce the same value for b, leading in turn to weaker constraints on the non-Gaussianity parameter h.This effect leads to particularly significant deviation from the constraints found using the approach of Young and Byrnes as the second term in eq.(4.32) is larger than that in eq.(4.33) by an order γ 2−n = O(10 n−2 ).

Conclusion
PBHs offer very tight constraints on non-Gaussianity even if they make up only a small portion of the DM.Because non-Gaussianity is an important prediction of many models of inflation, PBHs are a great window to study the early universe.
Previous work by Young and Byrnes [45] and Tada and Yokoyama [53] considered how the local-type non-Gaussianity parameters could be constrained from isocurvature modes appearing from fluctuations in the formation rate of PBHs, if PBHs make up a non-negligible fraction of dark matter.This was done by assuming that PBHs form in regions where the curvature perturbation exceeds a critical value and using Press-Schechter theory to calculate the abundance of PBHs.By considering modal coupling arising from non-Gaussianity, the peak-background split was applied to calculate the isocurvature modes appearing from fluctuations in the PBH formation rate.By then using constraints from the Planck Collaboration [46] on isocurvature modes they were able to constrain the non-Gaussianity parameters f = 3f NL,local /5 and g = 9g NL,local /25 for a given fraction of DM that is made up of PBHs, f PBH .
In this work, we have updated these constraints by using recent developments in the field, improving the calculation of the PBH abundance, and providing more accurate results.These include the use of peaks theory instead of Press-Schechter theory, the use of the compaction instead of the curvature perturbation to describe when PBHs form and the correct mass scaling of the PBH mass (which depends on the scale and amplitude of the perturbation which formed the PBH).Applying again non-Gaussianity when calculating the PBH abundance and applying the peak-background split to calculate the isocurvature modes, we found updated constraints on the non-Gaussianity parameters by using again the Planck data [46].Based on refs.[43,47], we consider only a narrow peak in the power spectrum to be responsible for PBH formation, ensuring the validity of the local-type expansion considered.The calculation could be extended to predict the isocurvature modes for specific models, such as the curvaton model considered in [42,43].
Similar to the results found by Young and Byrnes and Tada and Yokoyama, our updated calculation for quadratic models of non-Gaussianity offers tight constraints of |f | ≲ 10 −4 for f PBH = 1, when all of the DM is made up of PBHs.This means that the constraints earlier found by Young and Byrnes for quadratic models through Press-Schechter theory still provide a reasonable approximation in case the easier approach of Press-Schechter theory is desired.
For cubic models of non-Gaussianity the constraints in this work weaken significantly compared to those found by Young and Byrnes, however.This is due to additional suppression of cubic order non-Gaussian terms by a factor γ ∼ 0.1.For smaller values of g this broadens constraints by a factor O(10), while for larger values of g the equations start to become non-linear, and constraints worsen even more.Even so, for large values of f PBH , even cubic models of non-Gaussianity can still be tightly constrained by the updated calculation, with |g| ≲ 10 −3 for f PBH = 1.
We emphasise that the constraints calculated in this paper assume that the local model of non-Gaussianity is accurate over a large range of scales; from the large scales visible in the CMB to the small scales at which PBHs form.The constraints therefore apply specifically to the local-type non-Gaussianity parameters.A local-type bispectrum peaks in the squeezedlimit, implying a strong correlation between large and small scales -giving strong constraints on the non-Gaussianity parameter.However, other bispectrum shapes (such as equilateral and orthogonal) do not peak in this limit, implying a weak(er) correlation -which would result in much weaker constraints.
In the case of models predicting varying levels of non-Gaussianity on different scales, the constraints calculated here would not be accurate.For example, if the curvaton model is considered, then the small-scale perturbations can be sourced by the curvaton, which can be decoupled from the inflaton, which sources the large-scale perturbations.This would imply a much weaker correlation between the small and large scales, and correspondingly weaker constraints on the non-Gaussianity parameters.Such constraints could be computed on a model-by-model basis using the methods presented here, although such consideration is beyond the scope of this paper.However, the constraints on the non-Gaussianity parameters from PBH isocurvature modes could still be competitive with those available from other observations (such as the CMB) even if weaker by several orders of magnitude.
We also briefly considered the constraints on higher-order non-Gaussianity parameters.Young and Byrnes previously found that the isocurvature constraints applied (almost) equally at all orders of non-Gaussianity and would essentially rule out non-Gaussianity at all orders if PBHs are detected.However, we have demonstrated here that this is not the case, and that the constraints become significantly weaker as higher order terms are considered.

Figure 1 :
Figure1: The PDF of a Gaussian distribution where the average value is much larger than the standard deviation.Such distributions approach Dirac delta distributions.

. 35 )
Because C 1 takes larger values for type II perturbations than its does for type I perturbations, C G will also take larger values for type II perturbations by eq.(3.25) for |f | , |g| ≲ 0.1.The exponent in this integral then allows us to justify our earlier statement that type II perturbations are exponentially suppressed and therefore less relevant than type I perturbations.Using eq.(3.5) and (3.25) for either a quadratic or cubic model of non-Gaussianity we can solve the above equation if the range of C G is known.PBHs form when the compaction exceeds the critical value, C > C c .To find the value of C 1 corresponding to C = C c , we can invert eq.(3.5) to find the solution

. 39 )Figure 2 : 2 Figure 3 :
Figure 2: Perturbations with C 1 < 4/3 are referred to as type I perturbations, while perturbations with C 1 > 4/3 are referred to as type II perturbations.PBHs form when C > C c and we only consider type I perturbations.This means that we only consider the range 0.67 = C 1,c < C 1 < C 1,to = 4/3 in our calculations.

. 11 )Figure 6 :
Figure 6: Graph of the change in δ β as a function of small ζ l for different values of f in the quadratic expansion (left panel) and g for the cubic expansion (right panel).We use σ 0 = 0.15 for reference here.

Figure 7 :+ gζ 3 s
Figure 7: Constraints on f and g for a given value of f PBH as found in this work compared to those found by Young and Byrnes for quadratic expansions (left panel) and cubic expansions (right panel).Indicated are the allowed values for f and g for a given value of f pbh .
(4.22) and (4.23) that correspond to these values.We find that for a given ζ = O(1) and C 1 = O(1), the solution for ζ s and C G are found respectively by good approximation as ζ ≈ ζ s + gζ 3 s , (4.26)