Induced gravitational waves and baryon asymmetry fluctuations from primordial black hole formation

We consider black hole formation due to the gravitational collapse produced by large density fluctuations during an epoch of reheating with a stiff equation of state and calculate the induced gravitational wave spectrum. By considering the existing bounds on the total energy density of gravitational waves today, we find constraints on the parameter space of this scenario. We also calculate the lepton asymmetry generated by metric perturbations via the chiral gravitational anomaly present in the Standard Model and find that, once the electroweak sphaleron processes have taken place, the large spectrum of scalar perturbations responsible for black hole formation induces a peak in the baryon asymmetry fluctuations on small scales.


Introduction
The baryon-to-photon ratio observed in CMB experiments [1] as well as in measurements of the abundance of light elements produced during Big Bang nucleosynthesis (see e.g.[2]), leads to the requirement of some mechanism, known as baryogenesis, to produce an asymmetry between matter and antimatter in the early Universe.One possibility that has been extensively considered is that of leptogenesis [3], which consists in generating a lepton asymmetry at early stages which is later converted to a baryon asymmetry via sphaleron transitions, non-perturbative processes that take place in the electroweak sector of the Standard Model.Typical models require extending the particle content of the theory by adding, for instance, right-handed neutrinos [4].An alternative possibility is gravitational leptogenesis, which consists in exploiting the chiral gravitational anomaly [5,6] present in the Standard Model (that is, in the presence of only left-handed neutrinos 1 ), where R R = (1/2)ϵ µν ρσ R µναβ R ρσαβ denotes the contraction between the Riemann tensor and its dual, and the factor N L−R = 3 arises from the difference between the number of left-and righthanded neutrinos in the theory. 2By expanding the right-hand side of the above equation to quadratic order in perturbations, one obtains a term proportional to the product h ij ϕ, as well as terms quadratic in h ij , where ϕ denotes the Newtonian potential and h ij the transverse, traceless tensor perturbation of the metric.From the structure of the resulting terms one can check that, in order for the mechanism to work, a chiral gravitational wave background is required, as one would expect from the Sakharov conditions [9].Such a spectrum can be generated, for instance, during inflation, in models in which the inflaton φ contains a CP-odd component and couples to gravity through a term of the form f (φ)R R [7]. 3  It was recently suggested in [11] that, due to the fact that inflation is a stochastic process, there is actually no need to invoke the presence of these couplings to produce a lepton asymmetry via the above mechanism.Indeed, due to the dependence of the anomaly on the stochastic variables ϕ and h ij , we expect a non-vanishing variance for the lepton number density ⟨|n L | 2 ⟩ to be generated in different patches as inflation progresses.The average of the root-mean-square variance over some particular region therefore quantifies the expected asymmetry.Since these fluctuations are generated during inflation (out of thermal equilibrium) and the anomaly violates C, CP, and L (with the electroweak sphalerons providing the required B violation later on), the Sakharov conditions are fulfilled.As shown in [11], however, when this average is computed over the entire observable Universe today, the resulting asymmetry turns out to be extremely suppressed, rendering the proposal incapable of producing baryogenesis.Although unable to explain the observed background value of the baryon asymmetry, the presence of the chiral gravitational anomaly would nevertheless 1 Even if right-handed neutrinos are added to the Standard Model, this anomaly is still present once they are integrated out, at sufficiently low energies [7]. 2 In general, all chiral fermions in the theory contribute to NL−R, with the contribution of each particle weighted by its corresponding B − L factor (see e.g.[8]), but since the Standard Model contains an equal amount of left-and right-handed quarks, their contribution vanishes and B drops out.The same argument applies to charged leptons.
3 See [10] and the references therein for a discussion about the caveats of this mechanism, in particular the presence of ghost modes, and some proposed solutions.
lead to unavoidable fluctuations in the baryon number density on sufficiently large scales, 4 which could potentially be used as an observable to probe different models of inflation.We anticipate, however, that in the particular scenario discussed here the fluctuations are much smaller than the observed background value, so measuring them would likely require some additional enhancement mechanism.
The calculation presented in [11] makes use of the h ij ϕ term that arises after expanding eq. ( 1) in perturbations, so that the size of the fluctuations in the baryon asymmetry is proportional to the amplitude of both the tensor and scalar power spectrum, and becomes negligible if the tensor-toscalar ratio r is sufficiently small.A natural question is therefore whether an enhancement in either of these quantities could significantly increase the size of these fluctuations.Such an enhancement can be obtained, for instance, in single-field models of inflation that aim to generate a significant population of primordial black holes by introducing an inflection point in the potential (see e.g.[12][13][14] for particular implementations of this mechanism).In this class of models, the presence of the inflection point leads to a phase of ultra-slow-roll that enhances the curvature power spectrum, producing large energy overdensities once the perturbations re-enter the horizon after the end of inflation, which in turn induce gravitational collapse, generating the black holes.These black holes could account for the entirety of the observed dark matter provided their masses lie in the range where the upper bound comes from microlensing observations [15] and the lower one from their Hawking evaporation [16][17][18][19].Such a large scalar power spectrum would also source gravitational waves at second order in perturbations after the end of inflation [20][21][22], leading to a peaked gravitational wave signal which would also contribute to the variance of the lepton number density and which is completely independent of the primordial tensor spectrum arising from inflation.The spectrum of gravitational waves induced by large scalar perturbations during a radiation era has been extensively studied in the literature [21][22][23][24].However, since the resulting signal depends on the evolution of the scalar perturbations after they re-enter the horizon, we expect the result to change if the background fluid has a different equation of state.The case in which the Universe is dominated by non-relativistic matter before transitioning to the radiation era was studied in [24][25][26], and the scenario with a general background was considered in [27], where an enhancement of the spectrum for background fluids with stiff equations of state (that is, for w ≲ 1, with p = wρ) was reported.A similar setup with a scalar field rolling down an exponential potential was studied in [28].This paper has two objectives.The first is to connect the enhancement of the induced gravitational wave spectrum for stiff background equations of state to the abundance of primordial black holes. 5We consider a scenario in which primordial black holes form during a reheating stage 6 with an equation of state parameter 1/3 < w < 1 before the Universe transitions to the radiation era at some temperature T r .We find that, depending on the specific values of the equation of state, the transition temperature, and the scale at which the peak in the scalar power spectrum is located (therefore, on the mass and abundance of the black holes that form), the enhancement of the signal can be large enough to violate the existing bounds on the total gravitational wave energy density derived from CMB observations and the abundance of light elements produced during Big Bang nucleosynthesis [31,32], effectively reducing the parameter space of this scenario.We perform the calculation of the gravitational wave spectrum for two cases, one in which the transition between epochs is instantaneous, and one in which the stiff fluid gradually decays into radiation, and show that the resulting constraints on the parameter space for PBH formation depend only mildly on the smoothness of the transition.We improve upon the results of [27] by implementing a matching procedure for the transfer function of the scalar perturbations and the Green's function of the tensor modes in the sudden transition case akin to the one presented in [24] for the case of matter-domination, effectively taking into account the full time evolution of both quantities.We also perform a fully numerical calculation of the scalar transfer function and the tensor Green's functions in the gradual transition case, in contrast to [26]. 7The second objective is to calculate the baryon asymmetry fluctuations induced by the chiral gravitational anomaly in eq. ( 1) for the same scenario. 8We extend the results of [11] in two significant ways.The first is that we consider a peaked scalar power spectrum responsible for PBH formation, as opposed to a scale-invariant one, thereby enhancing the asymmetry reported there.The second one is that we also consider the purely scalar contribution to the asymmetry due to induced gravitational waves, effectively removing one of the essential ingredients in [11], namely, the need to have a non-vanishing gravitational wave background generated during inflation.We once again remark that the mechanism studied here is unable to produced the observed baryon asymmetry of the Universe, but it does allow us to predict a spectrum of fluctuations in this quantity which would be present in any model of PBH formation from single-field inflation using an inflection point in the potential, assuming only the matter content of the Standard Model.Moreover, the machinery presented here could be used to compute these fluctuations in other models of gravitational leptogenesis, such as the one in [7].
The paper is structured as follows.In Section 1 we expand eq. ( 1) in perturbations and derive the corresponding expressions for the lepton number density at each order.In Section 2 we derive the expressions for the black hole mass and abundance assuming that they form during a reheating stage with a stiff equation of state, and discuss the relevant constraints on the parameter space.In Section 3 we derive the induced gravitational wave spectrum when the transition between the reheating stage and the radiation era is instantaneous, and determine the constraint on the PBH masses due to the aforementioned bounds on the gravitational wave energy density.In Section 4 we repeat the calculation by assuming that the stiff fluid gradually decays into radiation.Finally, in Section 5 we determine the baryon asymmetry produced by the chiral gravitational anomaly in this scenario.

Gravitational lepton anomaly
In this section we expand eq. ( 1) to third order in perturbations and find expressions for the lepton number density n L at each order.Throughout the paper we denote second-order perturbations of the metric with bold symbols.The perturbed FLRW metric is, in conformal time dη = dt/a, where we have fixed the Newtonian gauge, so that E = B = E = B = 0, we have assumed that the first-order vector perturbations vanish, E i = B i = 0, and we have used a second-order vector gauge-transformation to set B i = 0.The second-order vector perturbation E i cannot be set to zero, since it is sourced by terms quadratic in first-order scalar perturbations by virtue of Einstein's equations.We have additionally assumed that no anisotropic stress is present, so that the two first-order Newtonian potentials are equal to each other, ϕ = ψ.The difference between the two Newtonian potentials at second order Φ − Ψ does not vanish in the absence of anisotropic stress, since it is sourced by terms quadratic in first-order scalars.
The quantity R R which appears on the right-hand side of eq. ( 1) can be easily checked to vanish at the background level.This quantity can be calculated at each order in perturbations by direct expansion using the above metric.However, the result can also be obtained by noting that every term in the expansion must contain four derivatives (with two coming from each factor of R) and the indices can only be saturated in a limited number of ways.For instance, this quantity vanishes at leading order9 (LO), since the only perturbations available are ϕ and h ij , and every possible contraction (for instance, ϵ ijk ∂ i ∂ j ∂ k ϕ) vanishes due to the antisymmetry of ϵ ijk .At the next-toleading order (NLO) we have only one possible term mixing scalar and tensor perturbations, where primes denote derivatives with respect to conformal time (′ = d/dη), as well as several possible terms mixing two tensor perturbations, such as We generically expect these terms to be suppressed with respect to the scalar-tensor ones (since the stochastic gravitational wave background produced during inflation is much smaller than the scalar power spectrum) and thus we do not consider them.
At NNLO, we focus on the situation in which the first-order tensor perturbation h ij is negligible and the only relevant contribution comes from the scalar-induced h ij , which is the case in models of PBH formation from single-field inflation.In this case, all terms containing only scalar modes (such as ϵ ijk ∂ i ϕ∂ j ϕ∂ k ϕ ′ or ϵ ijk ∂ i ∂ j ϕ∂ k Φ ′ ) can be easily seen to vanish due to the antisymmetry of ϵ ijk .Similarly, the scalar-vector terms vanish due to the fact that terms such as , since only the symmetric combination ∂ i E j + ∂ j E i appears in the metric.Thus, the only relevant term at this order is the scalar-tensor one since terms quadratic in h ij are of higher order, and we neglect terms mixing h ij and h ij , which are subdominant with respect to the above contribution.We conclude that, at NLO, the only relevant term is the one in eq. ( 4), assuming that the firstorder gravitational wave background produced during inflation is suppressed with respect to the scalar one at the scales of interest.The numerical prefactor can be obtained by explicit calculation, and the result is On the other hand, at NNLO only the term in eq. ( 6) contributes, assuming vanishing first-order tensor modes.The prefactor can be obtained by simply substituting h ij → 1 2 h ij in eq. ( 7), We remark that the NNLO terminology in no way implies that the term in eq. ( 8) is smaller than the one in eq. ( 7).The reason is that, since h ij is sourced by scalar perturbations, its amplitude will depend on their initial conditions and time evolution, which can be very different from those of the unsourced h ij .The left-hand side of eq. ( 1) can be expanded by writing J µ = (a −1 n L , 0) and using the following identity, We therefore obtain, for the NLO term, We can expand the perturbations in Fourier modes where e s ij denotes the two transverse, traceless polarization tensors, defined via where v and v are two unit vectors satisfying k We can therefore regard e s ij as a k-dependent quantity.We then find The mean value of this quantity clearly vanishes, as can be easily checked by expanding ϕ k and h s k in terms of creation and annihilation operators.The variance, however, does not, and will be computed in Section 5.
To find the expression for the NNLO term in Fourier space it is necessary to solve the equation of motion for the second-order tensor modes induced by scalar fluctuations, given, in momentum space and in the absence of anisotropic stress, by (see e.g.[22][23][24][25]) where H = a ′ /a denotes the conformal Hubble factor and the source term is, in the Newtonian gauge, where p and ρ are the background pressure and energy density, and we have used the boldface e s (k) to denote the matrix e s ij (k).Explicitly, we have p • e s (k) • p ≡ p i e s ij (k)p j .The solution to eq. ( 15) is given by where G k (η, η ′ ) is the Green's function of the differential equation, and we have assumed that inflation ends at η = 0 and the gravitational waves induced up to that point are negligible.This solution can be rewritten as where with Here we have used the definition of the scalar transfer function ϕ k (η) ≡ T ϕ (kη)ϕ k (0), obtained by solving the equation of motion for ϕ k , as well as the fact that the initial value of the Newtonian potential ϕ k (0) is related to the frozen curvature perturbation R k on superhorizon scales via [34] where w is the equation of state of the Universe at the time at which the initial conditions are imposed (that is, shortly after the end of inflation).In Section 3 we will give explicit expressions for the Green's function G k (η, η ′ ) and the transfer function.
At NNLO, the lepton number density is given by where we reiterate that we have neglected all terms involving the first-order tensor modes h ij .Using eq. ( 18) we can write this quantity as As we will see in Section 5, the mean value of this quantity once again vanishes, but the variance does not.In general, n L is given by the sum of both the NLO and NNLO contributions.

Black hole masses and abundance
Black holes form when sufficiently large perturbations re-enter the horizon after the end of inflation.Their abundance can be described using the Press-Schechter formalism, which states that collapse occurs only when the density contrast δ = δρ/ρ is above a critical value δ c , where β = ρ PBH /(γρ) is the ratio of the energy density in a Hubble patch that ends up in the form of PBHs to the total energy density at the time of collapse, with γ an O(1) factor encoding the efficiency of the collapse. 10In the above expression we have assumed the probability distribution of the fluctuations to be Gaussian and the smoothed variance of the density contrast σ can be related to the dimensionless primordial spectrum of curvature perturbations P R by means of the gradient expansion [36] where W (x) = e −x 2 /2 2/π is a window function, which we take as a Gaussian.The dependence of the critical value δ c in eq. ( 24) on the equation of state parameter w was analytically estimated in [37] to be Note that this formula is not valid for w ≪ 1/3, since in this case the non-sphericity and angular momentum of the collapsing cloud, neglected in [37], become relevant, see [38,39]. Figure 1: Left panel: contours depicting the black hole masses as a function of the reheating temperature T r and the scale k ♯ at which the peak in the power spectrum is located for w = 1, together with the amplitude of the power spectrum A ♯ necessary to obtain f PBH = 1 (labeled solid lines).The horizontal and vertical shaded regions represent the constraints on the PBH masses in eq. ( 2).The region with bottom-right to top-left shading represents the re-entry constraint in eq. ( 32), and the region with bottom-left to top-right shading represents the GW constraint in eq. ( 54).Right panel: same as left panel, but leaving w as a free parameter, for k ♯ = 5 × 10 12 Mpc −1 .
Recent studies have focused on using the peak theory formalism instead of Press-Schechter to calculate the PBH abundance, see e.g.[40].Similarly, the fact that in models of PBH formation the probability distribution deviates significantly from the Gaussian approximation is well-understood [41].However, as noted in [40,41], the change in abundance induced by taking into account these modifications can always be compensated by multiplying the power spectrum by an O(1) factor, since eq.( 24) depends exponentially on P R .Since the induced gravitational wave spectrum depends quadratically on P R , these corrections will not heavily impact the results presented in this work, and in any case can be incorporated by a simple rescaling of the power spectrum.The same argument applies to the threshold δ c , which strictly speaking should be calculated by finding the local maxima of the compaction function, a quantity that measures how close the overdensity is to fulfilling the hoop conjecture, see e.g.[42].The specific value of δ c will not be particularly important for our arguments, since the only issue relevant for us is the value of P R necessary to obtain Ω 0 PBH = Ω 0 DM (that is, for PBHs to form all of the dark matter today), which is always of order P R ∼ 10 −2 .We take eqs.(24,26) simply as a convenient benchmark, and we remark that the main results of this paper do not rely on the use of the Press-Schechter formalism at all.
In what follows we consider a situation in which the Universe goes through a long phase of reheating after inflation with an arbitrary equation of state w > 1/3 during which PBHs form, and then transitions instantaneously into a radiation-dominated era at temperature T r .By following a procedure analogous to the one presented in [43] for the case w = 0, we can write the present fraction of dark matter in the form of PBHs f PBH ≡ Ω 0 PBH /Ω 0 DM as where T 0 is the temperature of the Universe today, Ω γ the energy density in radiation, g ⋆ denotes the effective number of relativistic degrees of freedom, and g ⋆s the effective number of degrees of freedom in entropy. 11We have normalized the scale factor as a 0 = 1.In deriving the above equation we have assumed that collapse occurs immediately after perturbations with wavenumber k re-enter the horizon, at a time t k defined by k = a(t k )H(t k ).We assume entropy is conserved only between the transition time t r and today, but not necessarily between the time of collapse and t r .The mass of the PBHs is proportional to the mass contained in a Hubble patch at the time of formation, M PBH = 4πγM 2 p /H.By following a similar procedure to the one we used for f PBH , we can write As a cross-check, note that, for w = 1/3, eqs.( 27) and ( 28) become independent of T r and, in particular, we recover the well-known scaling relation It is clear from eqs. ( 24) and ( 25) that to produce a large PBH population one requires a spectrum of curvature perturbations P R which is peaked at a particular scale k so as not to overproduce black holes with masses outside of the unconstrained range in eq. ( 2).Such a power spectrum can be obtained from single-field models of inflation in which the potential has an approximate inflection point leading to a short phase of ultra-slow-roll, see e.g.[13,14].In what follows we will consider two possibilities, a flat scale-invariant spectrum with an amplitude A ♭ that matches the value measured on CMB scales of O(10 −9 ), and a sharp Dirac delta spectrum with an amplitude A ♯ of O(10 −2 ) in order to produce an O(1) fraction of PBHs as dark matter, where we have assumed that the peak occurs at a scale k ♯ .This form of the power spectrum will allow us to derive several results analytically later on.Putting everything together and considering the sharp power spectrum, we find the following expression for β at the peak of the distribution, which is located at k = k ♯ / √ 2 (as can be checked by setting dσ/dk = 0), Note that the total PBH abundance is obtained by integrating f PBH over k, but since we are considering essentially monochromatic mass distributions, we can approximate the abundance by the value of f PBH at its peak, . The black hole mass and abundance are plotted in Fig. 1.As anticipated, the power spectrum required to obtain f PBH = 1 is of order P R ∼ 10 −2 over the entire parameter space, so the wdependent threshold δ c does not have a particularly strong effect.In the Figures we also show the GW constraint derived in the next Section, in eq. ( 54), for a fixed value A ♯ = 0.05.Choosing the value of A ♯ necessary to obtain f PBH = 1 at each point in parameter space would only marginally change this constraint.We have also included the constraint imposed by the requirement that collapse occurs in the reheating stage.That is, the mode with comoving wavenumber k ♯ that induces collapse must re-enter the horizon at most at the time of the transition, We remark that this constraint does not impede the formation of black holes, but if it is not satisfied then collapse occurs during the radiation era, and one should use the expressions for the mass and abundance for w = 1/3 instead.We conclude that a population of PBHs with unconstrained masses which is large enough to explain the observed dark matter abundance can only form during an early stiff era (with w = 1) in a narrow region of parameter space, for 10 3 GeV < T r < 10 7 GeV, depending on the location of the peak in the power spectrum.We therefore find that the new GW constraint derived in this work (see the following Section) significantly restricts the parameter space for this mechanism.

Induced gravitational waves
In this Section we solve the equation of motion (15) for the tensor modes induced by scalar perturbations at second order.The solution, shown in eq. ( 18) can be written in terms of the Green's function G k (η, η ′ ) of the differential equation and the scalar transfer function T ϕ (kη), which we now determine.
The Green's function can be obtained by solving the homogeneous equation Since we assume that the transition between the reheating stage with equation of state parameter w and the radiation epoch occurs instantaneously at η r , the evolution of H is given by where we have chosen the integration constant η w ≡ (1 − 3w)η r /2 to make the function continuous at η r .We remark that throughout the paper w always denotes the equation of state parameter during the stage of reheating after inflation.We never use w to denote the equation of state parameter during the radiation era, and instead write explicitly p/ρ = 1/3 wherever necessary.The general solution to eq. ( 33) is where the ℓ subscript denotes quantities in the reheating stage, and r refers to quantities in the radiation era.We have also defined The Green's function is given by where h 1 and h 2 are any two linearly independent solutions to eq. ( 33) and we have suppressed the k and s indices in h for readability.Since this expression is independent of which pair of solutions is chosen, we can simply fix the integration constants in one of the regions by hand, and obtain the constants in the other epoch by imposing continuity of the solutions and their derivatives at η r .We therefore set A 1 hr = 1 and B 1 hr = 0 for the first solution, together with A 2 hr = 0 and B 2 hr = 1 for the second one.
For the calculation of the lepton number density we need the tensor transfer function, defined by We remind the reader that initial conditions are imposed at the end of inflation, which we choose as η = 0. Assuming that tensor modes are initially frozen outside the horizon (so that we can impose the initial conditions T h (0) = 0 and T ′ h (0) = 1 for superhorizon modes after inflation ends), we find the following expression for the tensor transfer function We remark that the transfer function is unrelated to the Green's function, and in this case the constants during the reheating stage are fixed by the initial conditions mentioned above, whereas the constants A hr and B hr must be determined by imposing continuity of the solutions and their derivatives at η r .This transfer function does not enter the calculation of the induced GW spectrum, but it is necessary for the calculations of Section 5, so we write it here for completeness.The equation of motion for the Newtonian potential is where we have assumed that the background is dominated by a perfect fluid with δp = (p/ρ)δρ.We remind the reader that we are working in the Newtonian gauge and in the absence of anisotropic stress.Assuming once again that the modes are initially frozen outside the horizon, we obtain the following expression for the transfer function, valid for w ̸ = 0, where and the constants A ϕr and B ϕr can once again be found by matching the solutions and their derivatives at η r .With these ingredients in hand we can calculate the I k function in eq. ( 19).The GW energy density is given by [44] where the brackets ⟨• • • ⟩ denote a time average, which must be taken due to the stochasticity of the signal, and P h denotes the power spectrum of h ij .Starting from eq. ( 18), we find the following expression for the power spectrum (see e.g.[23,24]) Since we are interested in measuring this quantity at late times, we can take η → ∞ in eq. ( 19).After pulling the solutions h 1 (η) and h 2 (η) in the Green's function out of the integral and splitting the limits of integration into the two contributions (0, η r ) and (η r , ∞), we obtain the following expression, valid late into the RD era, where and Analytical expressions can be found for the last two integrals, whereas the first two must in general be performed numerically.
Squaring I k and taking the average as required by eq. ( 43), we find where we have used ⟨sin 2 (kη)⟩ = ⟨cos 2 (kη)⟩ = 1/2 and ⟨sin(kη) cos(kη)⟩ = 0 and highlighted the dependence of J on the dimensionless parameter kη r .We can evaluate the integrals over y and z in eq. ( 44) by using the sharp power spectrum in eq. ( 30) necessary to obtain a significant PBH population.The result is where Θ H denotes the Heaviside function, and we have used the fact that these gravitational waves behave as radiation at late times to evolve the result from some late time deep into the radiation era until today.This expression only depends on the two dimensionless quantities k ♯ /k and k ♯ η r .The transition time is, in terms of temperature, In Fig. 2 we show the GW energy density from eq. ( 51) for different choices of the two parameters w and k ♯ η r .We find that the GW abundance grows for stiff equations of state, as reported in [27], due to the fact that the tensor modes redshift more slowly than the background.We also find that the energy density grows as k ♯ η r increases.The reason for this is that the largest contribution to the momentum integral in eq. ( 44) comes from scales around the narrow peak in the scalar power spectrum, and therefore we should only expect the signal to be affected if the relevant modes reenter the horizon before the transition to radiation has taken place, so that k ♯ > H r ∝ η −1 r .For k ♯ η r ≪ 1, the relevant modes enter during the radiation era and the effect of the stiff epoch on the signal is washed out, as can be clearly seen in the Figure .The total energy density can be found by integrating As we anticipated earlier, there is a bound on this quantity arising from CMB observations and the abundance of light elements produced during Big Bang nucleosynthesis [31,32], Since the signal grows as both w and k ♯ η r increase, this bound is eventually violated.In Fig. 3 we perform a numerical scan over the parameter space, showing the forbidden region.12By using eq.( 52), this bound can be written in terms of the transition temperature, and translated to a constraint on the PBH mass and abundance, shown in Fig. 1 for w = 1 (left panel) and for k ♯ = 5 × 10 12 Mpc −1 (right panel).These constraints are one of the main results of this paper.As anticipated at the end of the previous Section, this bound significantly limits the available parameter space for PBH formation during a stiff epoch.The constraint disappears completely for the standard radiation scenario with w = 1/3, where the bound is satisfied.In these Figures we have taken A ♯ = 5 × 10 −2 (thereby assuming that f PBH ≃ 1) and neglected the mild dependence of this number on w depicted in Fig. 1.
Before moving on, we remark that, although we have restricted our attention to stiff epochs due to the fact that PBH formation during a matter-dominated era is not particularly well understood at the time of writing (since the effect of angular momentum and non-sphericity of the collapsing cloud becomes important [38,39]), the calculation of the GW spectrum presented in this Section is also valid for the case w = 0, with the only modification being the scalar transfer function in eq. ( 41), which takes the simpler form T ϕ (kη) = 1 for η < η r .

Gradual decay into radiation
For the calculation in the previous Section, we assumed that the transition between the reheating stage and the radiation era occurred instantaneously.We now study the case in which the initial fluid decays gradually, thereby making the transition between epochs smooth.This scenario was studied in [26] for the particular case w = 0 and for a scale-invariant scalar spectrum.Here we extend these results for a generic fluid with constant equation of state parameter w and using a peaked scalar spectrum relevant for PBH formation.It was reported in [26] that, for a scaleinvariant scalar spectrum, the induced gravitational wave spectrum in the smooth transition case is suppressed with respect to its sudden counterpart.As we will see, although the same conclusion remains true for a peaked spectrum, the suppression is very mild and therefore does not strongly affect the bounds derived in the previous Section.
In what follows we denote quantities related to the radiation fluid by a subscript γ, and quantities related to the decaying fluid by a subscript w.We model the transition by considering a rate of energy transfer of the form Γρ w between the two fluids, with Γ constant, such that the total stress-energy tensor is conserved, but the separate components are not, so for some vector Q µ [45].The background equations then become, using the number of e-folds dN ≡ Hdη as the time variable, At the perturbative level, we work in Newtonian gauge and assuming no anisotropic stress.We use the following dimensionless set of variables where δq b = a(ρ b + p b )δv b , and δv b denotes the velocity perturbation for each fluid component. 13he Einstein equation that determines the evolution of the Newtonian potential ϕ is The continuity and Euler equations for the fluid are [45] where we have introduced the additional assumption that δp w = wδρ w , and similarly for the radiation component.Together, these five equations form a complete system that can be solved numerically.
To solve the above system of equations, we choose adiabatic initial conditions [34].On superhorizon scales k ≪ H and at sufficiently early times so that Γ ≪ H and ρ γ ≪ ρ w , the equation for the Newtonian potential ϕ becomes Since there is essentially only one fluid present at this time, and we have assumed δp w = wδρ w , the Newtonian potential can be shown to satisfy eq. ( 40) and is therefore frozen on superhorizon scales.Thus, on these scales, we have from the above equation the relation Figure 4: Left panel: evolution of the background energy density of each fluid, as well as the total energy density ρ and the equation of state p/ρ, for the case w = 1.By normalizing all quantities with respect to their value at the time N eq defined by ρ γ = ρ w , the plot becomes independent of Γ.Right panel: calculation of the GW energy density for the sudden transition case (solid) and the gradual case (dashed) for w = 1 with k ♯ η r = 1 (blue) and k ♯ η r = 100 (red).The dotted lines correspond to the gravitational reheating case with Γ = 0.The time η r in the gradual transition scenario is defined through k eq in eq. ( 69).The integrated energy density is very similar for both cases, leaving the bounds of the previous Section essentially unaltered.
The adiabatic initial conditions for the energy density perturbations are [34] Finally, if we assume that θ γ and θ w are also frozen on superhorizon scales, we obtain, from the last two equations, By using the relation ( 21), we can therefore fix all initial conditions in terms of R.
The evolution of the background energy density of each fluid is shown in Fig. 4 for w = 1.Note that, as explained in [26], the evolution of the background quantities (and in fact, that of the perturbations) is completely independent of the choice of Γ once they are normalized by their values at the time N eq defined by ρ γ (N eq ) = ρ w (N eq ) ≡ ρ eq , as can be checked from eqs. (56-58).The left panel of Fig. 5 shows the time evolution of the perturbations for particular values of the parameters.We define the transition time η r in this case via where k eq is the mode that re-enters the horizon at N eq .This is the relation between k eq and η r Left panel: evolution of the perturbations with adiabatic initial conditions for w = 0.5 and k ♯ η r = 10, with η r defined via eq.( 69).The vertical dashed line denotes the horizon crossing time for this mode.Right panel: evolution of the two independent solutions to the homogeneous equation for the tensor modes used to construct the Green's function (solid lines) for w = 1, k ♯ η r = 1, and k/k ♯ = 1.At late times the two solutions approach the asymptotic values in eq.(72) (dashed lines).
in the sudden transition case of the previous Section, so this definition allows us to make a fair comparison between both scenarios.
For the calculation of the gravitational wave spectrum we need both the scalar transfer function T ϕ , obtained by solving the system in eqs.(60-64) with initial condition T ϕ = 1 (the condition d dN T ϕ = 0 is guaranteed by the choice of adiabatic initial conditions for the other perturbations), and the Green's function for the tensor modes, which we now turn our attention to.The integral in eq. ( 19) can be rewritten in terms of the number of e-folds as The equation of motion for the tensor modes is, in terms of the number of e-folds, To calculate the GW spectrum today, we must evaluate eq.(70) at late times, deep into the radiation era.At this stage the background pressure is p = ρ/3, and the Hubble scales as H = H c (a c /a), where the c subscript denotes some late time N c ≫ N eq .The product a c H c approaches a constant value, which can be found numerically, after the transition.Using this scaling in the above equation, we find the following two independent solutions at late times By using these two solutions we can take the time-average of I 2 k as we did in the previous Section, and we obtain, after performing the momentum integral over the Dirac delta power spectrum, where with Once again we see that the GW energy depends only on k/k eq = (k/k ♯ )(k ♯ /k eq ), as well as the two parameters k ♯ /k eq and w.The dependence on the parameter k eq therefore replaces η r from the previous section, with the relation between both given by eq. ( 69).We remark that, although the Green's function does not depend on which two linearly independent solutions are chosen to construct it, the fact that we perform the time average over I 2 k by using ⟨sin 2 (kη)⟩ = ⟨cos 2 (kη)⟩ = 1/2 means we need to project the two solutions that behave as eq.( 72) at late times.This can be accomplished by imposing the boundary conditions h 1,2 (N late ) = h late 1,2 (N late ) at some late time N late ≫ N eq , sufficiently deep into the radiation era (numerically, we find that a few e-folds after the transition suffice).The corresponding solutions are shown in the right panel of Fig. 5 together with their late time limits for w = 1, k ♯ η r = 1 and k/k ♯ = 1.
We show the resulting GW energy density for w = 1 in the right panel of Fig. 4, for two illustrative examples with k ♯ η r = 1 and k ♯ η r = 100.We find that, although the signal is very mildly suppressed with respect to the sudden transition case of the previous Section, the difference between both is essentially negligible the purpose of estimating the curve in Fig. 3, so we conclude that the results from the previous Section are robust.We nevertheless remark that this conclusion might change if the transition is modelled differently.In the same panel, we show the gravitational reheating case with Γ = 0. We find that in this case the oscillations are also washed out, and the amplitude of the signal is comparable to the gradual transition case.The transfer functions obtained numerically using the procedure in this Section are compared to the results using eq.( 41) in Fig. 6.
Having established that the gradual transition scenario is not significantly different from the sudden case, we go back to using the analytical formulas of the latter in the following Section.

Smoothed lepton number density
In this Section we estimate the size of the baryon asymmetry fluctuations induced by the gravitational chiral anomaly.The quantity of interest is the variance of the lepton number density14 n L computed in eqs.( 14) and ( 23), averaged over a region of size r σ , where W rσ (r) is a Window function that decays smoothly on scales r ≫ r σ , which we take to be a Gaussian for concreteness, Let us clarify why we focus on the variance of this quantity.We would naively expect the baryon asymmetry to be determined by the mean value of n L .However, in the absence of a chiral gravitational wave background (generated for instance, via some inflationary coupling of the form f (φ)R R such as the one considered in [7], where φ is the inflaton field), which would make the terms quadratic in h ij of the form shown in eq. ( 5) nonzero, the mean value ⟨n L ⟩ vanishes.This can be easily seen at NLO, since, schematically, taking the mean value of eq. ( 14) yields ⟨ϕh⟩ ∼ ⟨ϕ⟩⟨h⟩ = 0.
The fact that the mean value also vanishes at NNLO is less obvious.The dimensionless bispectrum B R of curvature perturbations is defined by where we have introduced the shorthand δ 3 p ≡ δ 3 (p) for later convenience.This quantity vanishes if the fluctuations are Gaussian and is therefore slow-roll suppressed in conventional inflationary scenarios, but this could change depending on the dynamics of the inflaton (particularly in models of PBH production featuring a potential with a near-inflection point, in which the slow-roll approximation breaks down [41]).After some manipulation, we obtain, from eq. ( 23), where we have used the Dirac delta function in the definition of the bispectrum to perform one of the momentum integrals, as well as the normalization of the window function The right-hand side of eq. ( 81) clearly vanishes, since p j p k is symmetric, but ϵ ijk is not.In what follows we assume that the curvature perturbation R follows a Gaussian distribution.Instead of focusing on the mean value of this quantity, we assume that the baryon asymmetry of the Universe is generated via some other mechanism which we remain agnostic about, and shift our attention to the variance of n L . 15As noted in [11], since inflation is a stochastic process, we expect the lepton number density to deviate from its mean value in different patches.This deviation will generically have a magnitude of order ∼ ⟨|n L | 2 ⟩ rσ in a region of size r σ , leading to fluctuations in the baryon asymmetry.Before moving on with the calculation however, let us note that, as discussed in [11], this asymmetry does not survive at late times on arbitrarily small patches, due to annihilation processes.Let us consider two neighboring patches, one with a matter excess, and another one with an antimatter excess.If particles are able to freely travel from one patch to the other, annihilation will take place, leading to a smaller asymmetry overall.On sufficiently large patches, however, the asymmetry always survives at late times.To see why, note that the maximum distance a particle in the radiation bath can travel between collisions is λ/a, where λ is its mean free path.In a time ∆t, the particle undergoes N = ∆t/λ collisions.The average displacement in a time ∆t for a random walk is therefore ∆r = (λ/a) √ N = (λ/a) ∆t/λ.Integrating this displacement yields the Silk length [47], which sets the limit below which annihilation can take place. 16On the other hand, at late times, after the electroweak sphaleron processes have taken place, the asymmetry is carried by quarks, which are relativistic.Once the temperature is low enough, quarks are confined into non-relativistic baryons and their mean free path drops significantly.The Silk length is therefore approximately given by where λ Q = 1/Γ Q denotes the mean free path of quarks in the plasma, and Γ Q is the interaction rate, which can be estimated as Γ Q ≃ T [48].For simplicity, let us neglect the effect of the equation of state on the evolution of a and assume that the Universe is always radiation-dominated.We also neglect the change in the relativistic degrees of freedom in entropy g ⋆s between the end of inflation and the QCD scale.We then have T ≃ T 0 /a and H = H 0 a −2 , which leads to We remark that this is only an order of magnitude estimate, 17 and an accurate description of the damping dynamics requires solving the Boltzmann equation.We conclude that the asymmetry is preserved on scales larger than the Silk length of quarks at the confinement scale, r σ ≳ r S ≃ 10 −17 Mpc.This is the domain of validity of the calculations in this Section.Below this scale, annihilation dynamics become important and we expect the amplitude of the baryon asymmetry fluctuations to drop significantly.
In what follows we will compute the fluctuations in three different cases, 1. for the NLO term in eq. ( 14) with a flat scalar spectrum (29) and with the tensor spectrum given by the relation [49,50] is the tensor-to-scalar ratio, 2. the same NLO term with a sharp scalar spectrum (30), but with a flat tensor spectrum given by the same relation as in the first case, and 3. the NNLO term in eq. ( 23) for a peaked scalar spectrum, and with the tensor modes induced by scalar perturbations.As we will see, each term dominates at different scales, and the total asymmetry fluctuations will in general be given by the sum of the three contributions, plus the subdominant mixed terms that we have neglected as per the arguments of Section 1.We also neglect the term involving the flat part of the scalar spectrum together with the induced GW piece, since we expect it to be suppressed with respect to the NNLO term mentioned above.Before moving on with the technical details of the calculation, let us anticipate that, since the variance of n L is proportional to both the tensor and scalar power spectrum, from dimensional analysis we can guess the result in all three cases to be |n L | 2 1/2 rσ ∼ √ P R P h /(ar σ ) 3 , since r σ is the only dimensionful parameter.Dividing by the entropy density s ∼ g ⋆s (T )T 3 and using entropy conservation, we find up to an overall numerical factor.The NNLO term will be highly peaked around r σ ∼ k −1 ♯ , whereas the NLO term due to the flat part of the power spectrum will instead grow simply as r −3 σ and reach its highest value at the Silk scale r S , as per the previous discussion.
Straightforward evaluation of the correlation functions yields where we have introduced the shorthand δ 3 (k) ≡ δ 3 k and assumed that h Using the Dirac delta functions to perform two of the momentum integrals, we find where we have introduced the Fourier transform of the window function, By expanding the Levi-Civita symbols in terms of Kronecker deltas, we find the following expression for the term in brackets s ϵ ijk ϵ abc q j q ℓ p k e s iℓ (p)q b q d p c e s ad (p) = q 2 p 2 − (q If we choose the coordinate system in such a way that the z axis is aligned with p, then the vectors in the definition of the polarization tensors in eq. ( 13) are simply v = x and v = y.Using θ to denote the angle between q and the z axis, and φ for the corresponding azimuthal angle, we obtain, after some straightforward algebra, s ϵ ijk ϵ abc q j q ℓ p k e s iℓ (p)q b q d p c e s ad (p) = 1 2 p 2 q 4 sin 4 θ.
The window function is obtained by Fourier-transforming eq. ( 79), After switching to spherical coordinates, all of the angular integrals can be performed explicitly, and we arrive at the expression As explained earlier, the observable quantity of interest is the root mean square of the lepton number density.In order to compare with the baryon asymmetry after the electroweak sphaleron processes have taken place, we need to divide by the entropy density The main contribution to the time integrals in eq. ( 94) occurs during the reheating stage, since the transfer functions for both tensors and scalars decay quickly during the radiation era.Let us suppose that the integrals freeze at some time t f shortly after reheating ends.Then, all of the time dependence of the lepton number density is in the a 3 factor, and we can relate the value of n L /s to its value today by using entropy conservation, Moreover, the upper limit in the time integrals can then be taken as η → ∞.  99) as a function of the parameter η r /r σ for the different equations of state w = 0 (black, dashed), w = 0.5 (black, dot-dashed), and w = 1 (black, dotted), together with the analytical estimate in eq. ( 104) (red, solid).
So far, we have not made any explicit choice for the power spectra.Let us set We then obtain with where we have introduced the dimensionless variables as well as the following notation for the transfer functions which simply makes explicit the dependence on η r .There is, of course, also an implicit dependence on the equation of state w during the reheating stage.In fact, G ♭ NLO is a dimensionless numerical factor that depends only on the parameters w and η r /r σ , which can be found from eq. ( 52), The dependence on both of these quantities is quite mild, as illustrated by the numerical results presented in Fig. 8.It is instructive to obtain an analytical result for G ♭ NLO by introducing some approximations, following the procedure in [11].The first is that since the time integrals freeze around the end of reheating, we can take the upper limit as η r instead of ∞, which allows us to use the expressions for the transfer functions valid during the reheating stage.Moreover, since the dependence on the equation of state w is quite mild, we can simply set w = 0, so that T ϕ (kη) = 1.The time integral can then be performed immediately, yielding where we have suppressed the second arguments in the transfer functions, since we are using the expressions during the reheating stage, and no matching of coefficients is involved.We can finally make the assumption that η r ≫ r σ , so that T h (p σ η r /r σ ) → 0. The remaining integrals over p σ and q σ can then be performed analytically, and we find This result is very close to those obtained by numerically calculating the integral in eq. ( 99), as shown in Fig. 8, even for η r ≲ r σ and w ̸ = 0. Our result differs from the one in [11] in two ways.The first is that we obtain a different expression for the integral in eq.(99).The momentum integrals in [11] diverge in both the IR and UV, so the authors have calculated it by imposing cutoffs on both limits.Our integral, in contrast, converges without the need for cutoffs, so we obtain a finite result for all parameter choices.The second difference is that the evolution of the transfer functions during the radiation era was not considered in [11], and instead an analytical estimate similar to the one performed above was presented.We have taken the effect of this evolution into account in our calculation, and computed the integrals in eq.(99) numerically by varying the two parameters η r /r σ and w, confirming that the dependence of the result on these quantities is very mild.
We conclude that, apart from the mild dependence of G ♭ NLO on r σ , the resulting n L /s scales as r −3  σ .As discussed earlier, this result is only valid up to the Silk scale r S , where the annihilation damping becomes relevant.The maximum value it reaches is therefore, taking r = 0.01 [49,50] and G ♭ NLO = 2/3 for simplicity, 3 + 3w 5 + 3w   110) (solid black), together with the asymptotic limits in (111) (dashed and dotted red).Right panel: numerical calculation of the full integral in eq. ( 109) for k ♯ r σ = 3 (blue) and k ♯ r σ = 10 (orange) for w = 0 (dashed), w = 0.5 (dot-dashed), and w = 1 (dotted), together with the approximate expression depicted in the left panel for each case (solid red and dashed red).
The asymmetry over the entire Hubble patch today is much smaller, 3 + 3w 5 + 3w so we confirm, as claimed in [11], that these fluctuations alone cannot explain the observed value of n L /s ≃ 10 −10 .The same conclusion holds for the two following cases, but, as we will see momentarily, the resulting spectrum of fluctuations has a rich structure that could, in principle, allow us to probe different inflationary models if it were observable on small scales.
Case 2: NLO, sharp scalar spectrum (NLO♯) We now turn our attention to the cases in which the power spectrum is sharply peaked at some particular scale (not considered in [11]), and compute the corresponding enhancement to the baryon asymmetry fluctuations.The calculation in this case is exactly the same as the previous one up to the choice of the power spectra in eq.(97).We now set and obtain with Now the dimensionless factor G ♯ NLO depends not only on w and η r /r σ , but also on the parameter k ♯ r σ .It is once again instructive to obtain an analytical result using the approximations of the previous case.Setting w = 0 and following the same procedure, we obtain We remark that, as in the previous case, although we are interested in w > 1/3, the choice w = 0 is convenient because it allows us to obtain an analytical result.The numerical result for w ̸ = 0 is shown in Fig. 9, where the dependence on w can be seen to be very mild.This function is highly peaked around k ♯ r σ ≃ 3, for which we find G ♯ NLO ≃ 0.38.This quantity has the asymptotic limits The function and its asymptotic limits are shown in the left panel of Fig. 9.The numerical results for the full integral in eq. ( 109) are shown in the right panel of the same Figure for k ♯ r σ = 3 and k ♯ r σ = 10 as a function of η r /r σ and w.We once again find that the dependence of the result on these two parameters is very mild, so only the scaling with k ♯ r σ is relevant for our purposes.The dependence of the result on r σ is shown in Fig. 11 for values of k ♯ of interest for PBH dark matter, together with the NLO♭ and NNLO♯ contributions, see the discussion below.
Case 3: NNLO, sharp scalar spectrum (NNLO♯) We now turn our attention to the NNLO contribution to the baryon asymmetry fluctuations.The amplitude of this term is entirely determined by the initial conditions and evolution of the scalar perturbations, and is therefore present even if the tensor-to-scalar ratio r is vanishingly small.The calculation in this case is more involved than the previous ones, so we perform it in Appendix A. The resulting expression is 3 + 3w 5 + 3w where G ♯ NNLO is computed in eq. ( 121).The coefficient G ♯ NNLO is, in this case, a function of w and k ♯ η r , as well as k ♯ r σ .The numerical calculation of this quantity is shown in Fig. 10 for different values of the parameters.We once again find that only the dependence on k ♯ r σ is relevant, and for w = 1 and k ♯ η r = 1 the function can be approximated by G ♯ NNLO ≃ 10 −2.8 exp − k 2 ♯ r 2 σ .We adopt this expression for the following discussion.
The sum of the contributions in all three cases is shown on the left panel of Fig. 11 for the representative value k ♯ = 10 14 Mpc −1 , and on the right panel of the same Figure as a function of k ♯ .We have rescaled n L /s by the factor (k ♯ r σ ) 3 to make the peaked structure of the spectrum more clear.This is the main result of this Section.We find that the NLO♯ contribution dominates for r −1 σ ≪ k ♯ , and the NLO♭ contribution becomes relevant only for r −1 σ ≫ k ♯ .For r −1 σ ≃ k ♯ , the distribution is heavily peaked due to the NNLO♯ term, which yields an enhancement over the NLO♭ result of O(10 6 ).As the right panel shows, the amplitude of the fluctuations at the peak of the distribution increases with k ♯ .We note, however, that for k ♯ ≃ r −1 S , where the enhancement is largest, the PBH masses obtained are far too small, as can be checked from eq. ( 28), yielding a distribution of black holes that would have evaporated today [51] (though having a peak in the scalar spectrum at this scale is still possible, provided its amplitude is small enough so that PBHs are not overproduced).It is also clear from the Figure that the NLO♯ and NNLO♯ contributions do not change the distribution significantly on large scales, so the estimate in eq. ( 106) of the asymmetry on a Hubble-sized patch today holds in the presence of the new terms.

Conclusions
We have calculated the expressions for the mass and abundance of PBHs produced during an early epoch of reheating with a stiff equation of state.We find that the parameters that determine the black hole distribution are the scale at which the peak in the power spectrum is located k ♯ , the temperature at which the transition to radiation occurs T r , and the equation of state w.There are three relevant constraints in this scenario.The first is that, in order to reproduce the observed dark matter abundance, the black hole masses must be in the range (2).The second is that for collapse to occur before the radiation era, the scale at which the peak in the spectrum is located must re-enter the horizon before the transition occurs, see eq. (32).Finally, since tensor modes are enhanced in the presence of a stiff epoch, we must also ask that GWs are not overproduced so that the bound on their energy density today, which arises from BBN and CMB observations, is not violated, see eq. ( 54).The allowed region of parameter space is shown in Fig. 1.
In order to calculate the induced GW signal we have implemented a matching procedure for both the Green's function of the tensor modes and the transfer function of the scalar perturbations, thereby extending the results of [27] by taking into account the full time evolution of these quantities.We confirm that, in the presence of a stiff epoch, the induced GW signal is enhanced.We have explicitly checked that the smoothness of the transition does not significantly change our results, and provided a procedure to compute the signal in these gradual decay scenarios exactly, without resorting to approximate analytical expressions of the Green's and transfer functions, extending the results of [25] to decaying fluids with stiff equations of state.As mentioned above, we have used these results to translate the bound on the GW abundance to a constraint on the PBH masses formed in this scenario.The results are shown in Figs. 1 and 3.
Finally, we have computed the chiral gravitational anomaly to third order in perturbations, and we find that the large scalar spectrum responsible for PBH formation induces a peak in the baryon asymmetry fluctuations on small scales.These results are shown in Fig. 11.We have shown that this spectrum is essentially independent of the cosmological history (i.e. the reheating scale and equation of state) and is therefore a generic prediction present in every model of PBH formation from collapse induced by large density perturbations, even in the standard scenario where the black holes form during radiation domination, assuming only the matter content of the Standard Model.
These fluctuations could, in principle, be used as an observable to probe not only the existence of PBHs, but also different models of inflation.Assessing whether these fluctuations can be measured by any future experiments is beyond the scope of this paper, but we point out that the main obstacle to this end is the fact that they are much smaller than the observed background value, and therefore some sort of enhancement mechanism would likely be required in the particular scenario studied here.We remark, however, that the machinery we have developed could also be applied to other models of gravitational leptogenesis, such as the one in [7], in which the fluctuations might be less suppressed.
The possibility of generating baryon isocurvature perturbations via the chiral gravitational anomaly was briefly mentioned in [11].Since we remain completely agnostic about the inflationary model giving rise to the large curvature perturbations responsible for the peak in the power spectrum in our scenario, no definitive statements can be made about this possibility.We remark, however, that these fluctuations would only arise in multi-field models of inflation (or, potentially, in models in which the inflaton couples directly to R R, such as [7]).In the minimal single-field inflation scenario, the arguments in [52][53][54] apply, and isocurvature would be generated neither during inflation, nor during the subsequent reheating era. 18The possibility of using different models of inflation to generate isocurvature would have to be studied on a case-by-case basis.

NNLO
In this Appendix we compute the quantity G ♯ NNLO relevant for the calculation of the lepton number density variance in the third case of Section 5.By putting together eqs.( 23) and ( 78) we obtain We can evaluate the six-point function by using Wick's theorem 2(2π) where we have used the symmetry of the integrand under k → p − k and k → p − k to simplify the result. 19We can use the Dirac delta functions to perform the momentum integrals over k, k and q.After using the symmetry of I k under the exchange of the two momenta in the argument, we obtain where we have also renamed the remaining dummy variable p → k.By choosing the sharp power spectrum in eq. ( 30) and following the same procedure as for the two NLO cases in Section 5, we find eq.( 112), with st q • e s (p) • q (p + q) • e t (k) • (p + q) F st (p, q, k) W rσ (|p + q|) 2 + F ts (−k, q, −p) where F st (p, q, k) ≡ ϵ ijk q j q ℓ p k e s iℓ (p) ϵ abc (p b + q b − k b )(p d + q d )k c e t ad (k) .
To obtain this expression from eq. ( 115) we have renamed the dummy variables k → −p and p → −k, as well as s ↔ t, in the second term inside the brackets.Let us orient our coordinate system in such a way that k coincides with the z axis, and denote the angle between p and q by Θ, and the angle between p + q and k by Φ.We can then use the following property of the Dirac delta function, We now switch to spherical coordinates in all of the momentum integrals in eq. ( 116).We denote the polar and azimuthal angles of p by θ p and ϕ p , respectively, and similarly for q and k.Since k coincides with the z axis, we can immediately perform the integrals over θ k and ϕ k .We can also perform the integrals over the moduli by using the Dirac delta functions and including a Heaviside function for each to take into account the fact that the moduli must be positive.We find G ♯ NNLO = 1 16π 2 dθ p dϕ p dθ q dϕ q sin θ p sin θ q Θ H (− cos Θ)Θ H (cos Φ) q • e s ( p) • q (−2 cos Θ p + q) • e t ( k) • (−2 cos Θ p + q) F st ( p, q, k) W k ♯ rσ (| q − 2 cos Θ p|) 2 + F ts (− k, q, − p) W k ♯ rσ (| q − 2 cos Φ k|) 2 , (121) where x ≡ k ♯ η, the hatted vectors are normalized to unity (and bear no relation to the hatted dummy variables in eq. ( 113)), and The window functions can be explicitly written as We can also write the Θ and Φ angles explicitly as follows, cos Θ = sin θ p cos ϕ p sin θ q cos ϕ q + sin θ p sin ϕ p sin θ q sin ϕ q + cos θ p cos θ q , (124) cos Φ = cos θ q − 2 cos Θ cos θ p .
With these matrices in hand, the dot products in eqs.(126, 127) can be written in terms of θ p , θ q , ϕ p , and ϕ q .Let us turn our attention to the time integrals in eq. ( 121), The I q function is given by eq. ( 19), so that with y = qη ′ , together with and In these equations we have highlighted the dependence of the transfer function T ϕ on the parameter k ♯ η r .Similarly, the factor q/H in eq. ( 134) is a function of y and qη r = (q/k ♯ )k ♯ η r .For the Green's function, we think of independent solutions h 1,2 and the functions G 1,2 defined in eqs.(76, 77) as functions of (qη, qη r ).We therefore find that I q is a function of (x, q/k ♯ , k ♯ η r , w).
Computing the 8-dimensional integral in eq.(121) numerically is difficult, even with Monte Carlo methods.The strategy we adopt is to first calculate the integral I for different choices of the parameters w and k ♯ η r as a function of q/k ♯ .The results are shown on the right panel of Fig. 12.Since q/k ♯ lies between 0 and 2 in eq. ( 121), we restrict our calculation to this range.The integrand of I is shown on the left panel of the same Figure .With these integrals in hand, we can then calculate the remaining four integrals over the angles in eq. ( 121) numerically by varying the remaining parameter, k ♯ r σ .The results are shown in Fig. 10.

Figure 2 :
Figure 2: Left panel: gravitational wave energy density as a function of k/k ♯ for w = 0.5 varying the dimensionless parameter k ♯ η r .The energy density grows as this parameter increases.Right panel: same as left panel, for w = 1.

Figure 3 :
Figure 3: Constraint in eq.(54) as a function of the two parameters w and k ♯ η r .The top-right region is forbidden due to the fact that the redshift of the tensor modes during an epoch with a stiff equation of state makes the energy density grow, as illustrated in Fig.2.
Figure 5:Left panel: evolution of the perturbations with adiabatic initial conditions for w = 0.5 and k ♯ η r = 10, with η r defined via eq.(69).The vertical dashed line denotes the horizon crossing time for this mode.Right panel: evolution of the two independent solutions to the homogeneous equation for the tensor modes used to construct the Green's function (solid lines) for w = 1, k ♯ η r = 1, and k/k ♯ = 1.At late times the two solutions approach the asymptotic values in eq.(72) (dashed lines).

Figure 7 :
Figure 7:Schematic depiction of the late-time behaviour of the baryon asymmetry fluctuations.Small neighbouring patches with matter and antimatter excesses annihilate, homogenizing the distribution on small scales.Larger regions remain unaffected, since quarks cannot travel from one patch to the other if the distance separating them is larger than the Silk length (see text).This length scale sets the range of validity for the calculations in this Section.

Figure 9 :
Figure 9: Left panel: approximate expression for G ♯ NLO in eq.(110) (solid black), together with the asymptotic limits in (111) (dashed and dotted red).Right panel: numerical calculation of the full integral in eq.(109) for k ♯ r σ = 3 (blue) and k ♯ r σ = 10 (orange) for w = 0 (dashed), w = 0.5 (dot-dashed), and w = 1 (dotted), together with the approximate expression depicted in the left panel for each case (solid red and dashed red).

Figure 11 :
Figure 11: Left panel: sum of the NLO♭ (solid, blue), NLO♯ (solid, orange), and NNLO♯ (solid, red) contributions to the baryon asymmetry fluctuations, represented by the dashed black line and scaled by the factor (k ♯ r σ ) 3 .The plot is cut at the Silk scale r S , where the calculation stops being valid.We have taken k ♯ = 10 14 Mpc −1 as a representative value for PBH dark matter (see Fig. 1), together with A ♯ = 5 × 10 −2 , r = 0.01, A ♭ = 2 × 10 −9 and w = 1.Right panel: sum of the three contributions to the baryon asymmetry for the same parameters as in the left panel, but varying k ♯ .

Figure 12 :
Figure 12: Left panel: integrand of I (I int ) for the specific parameters shown in the label.The integral converges quickly for x ≳ 4. Right panel: integral I as a function of q/k ♯ for w = 1 (black) and w = 0.5 (red) for k ♯ η r = 1, 10, and 100 (solid, dashed, and dotted lines, respectively).