Primordial black holes and gravitational waves from dissipation during inflation

We study the generation of a localized peak in the primordial spectrum of curvature perturbations from a transient dissipative phase during inflation, leading to a large population of primordial black holes. The enhancement of the power spectrum occurs due to stochastic thermal noise sourcing curvature fluctuations. We solve the stochastic system of Einstein equations for many realizations of the noise and obtain the distribution for the curvature power spectrum. We then propose a method to find its expectation value using a deterministic system of differential equations. In addition, we find a single stochastic equation whose analytic solution helps to understand the main features of the spectrum. Finally, we derive a complete expression and a numerical estimate for the energy density of the stochastic background of gravitational waves induced at second order in perturbation theory. This includes the gravitational waves induced during inflation, during the subsequent radiation epoch and their mixing. Our scenario provides a novel way of generating primordial black hole dark matter with a peaked mass distribution and a detectable stochastic background of gravitational waves from inflation.


Introduction
Primordial black holes (PBHs) with masses 10 −16 to 10 −12 times lighter than the Sun may constitute the totality of the dark matter of the Universe [1][2][3][4][5].Several mechanisms have been proposed to explain how such objects could have formed in the early Universe.The most popular among them is the gravitational collapse, during the radiation epoch, of Hubble-sized regions with a density contrast above an O(1) threshold.These regions can originate from large curvature fluctuations produced during inflation.In the simplest approximation, which assumes that these fluctuations are Gaussian, the cosmological abundance of PBHs is exponentially sensitive to the primordial power spectrum of curvature fluctuations, P R .Within this approximation, it is required that P R ∼ O(10 −2 ) for PBHs to be all the dark matter.Interestingly, if such a large value of P R is reached at the comoving scales that correspond to the aforementioned PBH mass range, a stochastic background of gravitational waves (GWs) with frequency approximately peaking between 10 −3 Hz and 10 Hz is also produced [6].This range contains the frequencies at which LISA (to be launched in 2037) and other spacebased interferometers (currently being discussed) such as DECIGO and BBO, are planned to have their best sensitivities, offering an indirect handle on the possible existence of an abundant PBH population [7,8].
In this work we discuss the generation of a large P R , and its associated stochastic background of GWs, from a transient dissipative phase during inflation.This mechanism of producing PBHs from large density fluctuations of inflationary origin is fundamentally different from the others that have been previously proposed.The most popular among them assumes a brief phase of ultra slow-roll, during which the acceleration of the inflaton is rather abruptly diminished, leading to a temporarily-growing super-horizon curvature mode, see e.g.[9,10].In our case, the effect of dissipation on the homogeneous inflaton dynamics can be described under adequate conditions by introducing a dissipative coefficient Γ in its time evolution: where ϕ denotes the inflaton field, H (which is approximately constant) is the Hubble function describing the growth of the scale factor during inflation and V ϕ is the first derivative of the inflaton potential (with respect to the inflaton field).During ultra slow-roll (with Γ = 0) the equation above becomes φ + 3H φ ≃ 0. Instead, during a strongly dissipative regime (Γ ≫ H) it reads φ + Γ φ ≃ 0. One may naively think that both regimes are analogous and the deceleration of the inflaton must have similar effects on P R in the two cases, enhancing it at specific scales.This intuition is flawed.
In the absence of other effects, the presence of Γ alone in the equations for the fluctuations makes the primordial spectrum of curvature perturbations decrease [11].However, an enhanced P R can arise in a strongly dissipative regime due to a stochastic source for the fluctuations coming from a thermal bath originating from the couplings between the inflaton and other fields that are also responsible for Γ.The framework we use to describe the dynamics of dissipation during inflation is essentially the one of warm inflation [12].However, whereas in warm inflation dissipation is active throughout the complete duration of inflation (which ends when the energy density of the radiation bath overcomes that of the inflaton), we instead assume a brief duration (∼ O(1) e-folds) for the dissipative phase.PBH production in warm inflation has been considered in [13,14].In those models the primordial power spectrum P R (k) appears to grow towards the smallest comoving wavelengths (∼ 1/k) exiting the horizon at the end of inflation and re-entering shortly after.That leads to very light PBHs (below ∼ 10 −27 Solar masses) that cannot account for the dark matter of the Universe regardless of their abundance, as they would have evaporated by now due to Hawking radiation.Given that the mass of PBHs formed during radiation domination from the collapse of large overdensities scales like k −2 , PBHs in the appropriate window for dark matter may be obtained provided that dissipation occurs only at the right time during inflation, i.e. around ∼ 30 e-folds after CMB scales become super-horizon, producing a primordial spectrum with a well defined peak at k ∼ 10 14 Mpc −1 .
We model our scenario using a phenomenological approach, assuming a peaked Γ as a function of the inflaton background field (or, equivalently, any other clock), which is also proportional to the third power of the temperature of the radiation bath (as it is common in warm inflation).We provide an analytical understanding of the features of the primordial spectrum and solve numerically the stochastic differential equations for the perturbations of the inflaton field, the metric and the radiation density, using two different methods.One of them consists in a Langevin approach, in which we integrate the system of stochastic linearized equations for multiple realizations of the thermal noise, obtaining a histogram for P R (k) at each k.The other method reduces the stochastic system of differential equations to a deterministic one by focusing on the computation of the relevant two-point functions.We refer to this second method as the matrix formalism.The results of both methods agree with high precision and are in broad agreement with the analytical understanding of the system.Both numerical methods, which we discuss in detail, can be useful in other scenarios in which dissipation (and hence stochasticity) occurs during inflation, including in particular the case of "standard" warm inflation.
Having obtained the primordial power spectrum of curvature fluctuations, we then compute the stochastic background of GWs induced at second order in cosmological perturbation theory.We take into account not only the GWs generated during the radiation phase that we assume follows after inflation, but also the GWs sourced at second order during inflation itself (which happen to be subdominant).There is also an interference term between both contributions, whose expression we derive analytically and whose value we estimate numerically.
We start in the next section by presenting the key features of the framework and the scenario we consider, deferring to an appendix the introduction of a toy model that can lead to a peaked Γ like the one we use in our analysis.We leave for future work a more detailed search for concrete Lagrangians.
Throughout the paper we set c = ℏ = k B = 1.

Dissipation during inflation
We assume that for a period of time during inflation the Universe contains a non-negligible thermalized radiation component with energy density ρ r and temperature T related by where g ⋆ quantifies the effective number of relativistic degrees of freedom.The inflaton field ϕ acts as the source for this thermal component.At the background level, this dissipation is assumed to be adiabatic, according to the following equations: φ + (3H + Γ) φ + V ϕ = 0, ( ρr + 4Hρ r = Γ φ2 , (2.3) which satisfy the conservation of the total energy-momentum tensor.The dots denote derivatives with respect to cosmic time (• ≡ d/dt) and M p ≃ 2.45 × 10 18 GeV is the reduced Planck mass.As long as friction eventually comes to dominate the dynamics for sufficient time, the initial conditions for these equations are irrelevant due to the presence of an attractor characterized by the ratios We show these ratios on the right panel of Fig. 1, together with the slow-roll parameters where defines the number of e-folds of inflation.The figure corresponds to the benchmark example that we introduce in Section 2.1.The end of inflation, defined by ϵ = 1, would typically occur for a number of e-folds after the peak scale crosses the horizon of N − N peak ∼ 25.For this example, the dissipation coefficient Γ transiently grows in value for ∼ 4 e-folds, becoming even Γ ≫ H, as shown in the left panel of the same figure.This suppresses the kinetic energy of ϕ while it enhances the radiation energy density (see middle panel).
In order to describe the dynamics of cosmological fluctuations we use the Newtonian gauge, where the two scalar perturbations of the metric are identical (and denoted as ψ) by virtue of one of Einstein's equations and the absence of anisotropic stress.We denote by δϕ and δρ r the perturbations of the inflaton field and the radiation energy density, respectively.Defining the momentum fluctuation where δv r is the velocity perturbation of the radiation, the remaining Einstein's equations in Fourier space and at linear order are (see e.g.[15]) (2.12) Figure 1: Left: T , H and Γ in units of M p and as functions of the number of e-folds N , near the value N = N peak which can be approximately associated to the comoving wavenumber k peak at which the peak in the primordial curvature spectrum (shown in Fig. 4) occurs.The analytical expression for Γ as a function of the inflaton background field is given in Eq. (2.21).The parameters chosen for this benchmark example are given in Eq. (2.22).Center: The inflaton potential V = λϕ 4 /4 (green), its kinetic energy density, the radiation density and 3M 2 p H 2 = V + φ2 /2 + ρ r (black dashed).

The conservation of the energy-momentum tensor
= 0 is satisfied in agreement with the continuity equations for ϕ and ρ r [16]: where Q µ contains a stochastic piece induced by thermal fluctuations, whose form is determined by the so-called fluctuation-dissipation theorem (see e.g.[17,18] and Appendix F), In this expression u ν denotes the 4-velocity of the radiation component and ξ t ≡ dW t /dt, where dW t is a Wiener increment Here, ⟨• • • ⟩ S denotes a stochastic average over different realizations.The linearized equations for δϕ, δρ r and δv r are (see e.g.[18]) ) ) where Γ T ≡ ∂Γ/∂T .
1 See Appendix C for the definition of a Wiener process and a quick review on stochastic differential equations.
Eqs. (2.16)-(2.18),together with one of Einstein's equations, for instance (2.11), form a complete set for the four variables δϕ, δρ r , δq r and ψ.These equations can be further simplified using the following combination of Einstein's equations, Imposing this constraint allows to reduce the number of equations by one, so we can eliminate, for instance, Eq. (2.18).However, we find that not imposing this constraint can be more stable numerically, as we discuss in more detail in Section 3.1 and Appendix E. We use the initial conditions where the initial condition for δϕ comes from assuming that the field fluctuations are in the Bunch-Davies vacuum at early times.As we will show later, the choice of initial conditions is not very relevant, since the noise term leads to an attractor behaviour for the evolution of the perturbations.

A peaked dissipative coefficient
Our main goal is the description of a peak in the curvature power spectrum arising from transient dissipation.The perturbation equations are driven by a source of noise with amplitude ∼ √ ΓT , and so they are significantly enhanced whenever Γ is sufficiently large.If the peak of the spectrum of the curvature perturbation is localized around an adequate scale, the PBH mass function will be narrow enough so that the borders of the allowed window for PBH dark matter can be avoided.Therefore, we focus on modeling a dissipative coefficient Γ that satisfies Γ ≫ H only for a few e-folds at most.Rather than building a full model of the complete inflationary history we focus on the local description of the dynamics around the relevant region.Although we remain agnostic about the details of the microphysics that gives rise to such a peaked dissipative coefficient, in Appendix A we present a toy example of a Lagrangian with the necessary features that could potentially serve as a basis for future models (which should also fit the CMB and provide enough inflation).
We assume the following parameterization of the dissipative coefficient where m, M, Λ and ϕ ⋆ are free dimensionful parameters.As discussed in Appendix A, the T 3 dependence of Γ arises naturally in a specific low-temperature limit, which is common in warm inflation.The temperature dependence of Γ is not crucial for the stochastic noise to generate a peak in the primordial power spectrum.A temperature-independent Γ that is peaked as a function of ϕ also produces a similar effect.However, the parameterization (2.21) resembles more closely the actual Γ that may be expected from a concrete Lagrangian in which ϕ couples to other fields.For our benchmark example of Fig. 1 we choose the following set of parameters: where λ is the coupling of a quartic inflaton potential, We choose this potential for its simplicity and lack of features. 2Since our focus is on studying the phenomenology of the dynamics generated by dissipation (at the background and perturbation levels) we could have chosen any other potential compatible with slow-roll inflation that does not have any peculiarities that would introduce spurious effects beyond those we want to analyze.In addition, the potential V (ϕ) only needs to be valid a few e-folds before and after the region where Γ ≫ H because we are only concerned with describing the appearance of a large peak in the primordial power spectrum, which is a local feature.Nevertheless, the value of λ is chosen in such a way that we obtain the correct amplitude for the power spectrum, A s ≃ 2 × 10 −9 [20].Our choice of parameters leads to a P R with a peak value of ∼ 10 −2 .As seen in Fig. 1, the potential and the evolution of H are essentially featureless, whereas the kinetic energy of the inflaton does change significantly in the relevant region.
For the initial conditions of the background variables we choose although the last two are essentially irrelevant due to the presence of the background attractor discussed in the previous section.This choice makes the background quantities converge quickly to their attractor values.The time at which the localized growth in Γ occurs (and therefore the scale at which the peak in the power spectrum is located) can be controlled by varying ϕ ⋆ .Decreasing m or M makes the peak of P R larger.In particular, since we choose m ≪ M , decreasing m makes P R increase without changing the value of Γ far away from the wavenumbers associated to ϕ ⋆ , so that P R retains its normalization at small distance scales.Similarly, increasing Λ makes P R larger.Finally, decreasing g ⋆ makes the peak of P R larger.To understand why this is the case, let us determine how the coefficient in front of the thermal noise in the equation of motion for the perturbations, (2.16) and (2.17), scales with g ⋆ .This can be done by isolating the temperature dependence of this quantity.Let us define γ(ϕ) via Γ(ϕ, T ) = γ(ϕ)T 3 .Then, by assuming the system is in the attractor solution (2.5) and Γ ≫ H, we find . (2.25) We therefore find that decreasing g ⋆ makes the amplitude of the stochastic noise increase, thereby increasing the curvature power spectrum.The effect of varying m and g ⋆ on the spectrum is shown in Fig. 2.

The curvature power spectrum
In this section we present three different ways of computing the primordial power spectrum of the comoving curvature perturbation 2 Notice that previous works have shown that the presence of radiation can allow inflation with a quartic potential to be compatible with Planck constraints on the tensor-to-scalar ratio and scalar index [19].where p and ρ are the pressure and energy density of the total system (inflaton plus radiation).Due to the presence of the stochastic thermal noise, our main quantity of interest is the expectation value of the power spectrum at a given comoving scale, which we denote by ⟨P R (k)⟩ S .The most straightforward way to compute this quantity (though not necessarily the most economical) is to determine (3.1) for a large sample of stochastic realizations, and then calculate their average.Alternatively, one can bypass this noting that ⟨P R (k)⟩ S is determined by the correlation of the thermal noise, and it is therefore a deterministic quantity.The system of stochastic differential equations can then be rephrased as a deterministic system for the correlators of the scalar fluctuations, which is convenient to write in matrix form.As we show below, both approaches agree.Finally, under some approximations, the equations for the fluctuations can be solved analytically, allowing us to understand the qualitative features of the power spectrum better.

Matrix formalism
In this section we will develop a method to find the stochastic average of the power spectrum by solving a deterministic matrix differential equation instead of the full set of stochastic differential equations.We begin by noting that the equations of motion can be written, in Fourier space, as a system of linear first-order complex stochastic differential equations.Throughout this section we will work with the number of e-folds as time variable and we define the 'column vector'3 The equations of motion (2.11), (2.16), and (2.17) can then be conveniently written as a system of four first-order stochastic differential equations where the matrix A and the (column) vector B are real and independent of Φ. Explicit expressions will be given at the end of this section.We also assume that the constraint in Eq. (2.19) has been imposed to eliminate δq r from the system.In this equation, ξ N denotes the Wiener increment from Eqs. (2.16) and (2.17) written in terms of the number of e-folds 4 and satisfying, in Fourier space, We are interested in the power spectrum of the comoving curvature perturbation, which can be written as where the vector C can be read from Eq. (3.2) (see Eq. (3.13)).The corresponding power spectrum, averaged over stochastic realizations, can be expressed in terms of the correlation function matrix It makes no difference whether we work with the real and imaginary parts of Φ, or with Φ and its complex conjugate Φ ⋆ .We now choose the latter option.The probability density P (Φ, Φ ⋆ , N ) for the system to be in state {Φ, Φ ⋆ } at time N can be obtained by solving the Fokker-Planck equation 6

∂P ∂N
where the two-point statistical moments are defined as The equation of motion for Q can be found differentiating the previous expression with respect to time and using the Fokker-Planck equation. 7The resulting deterministic differential equation for the matrix Q is By solving this deterministic differential equation we can bypass solving the full system of stochastic differential equations for the perturbations as long as we are only interested in the stochastic average of the power spectrum, which is given by Eq. (3.7).
Let us give explicit expressions for each one of the matrices used in these equations with the number of e-folds as the time variable.The matrix A is given by where The vectors B and C are Finally, the matrix of initial conditions Q i ≡ Q(N ini ) is, in accordance with Eq. (2.20), where k i is the scale that crosses the horizon at some initial e-fold value N ini .In practice, we can start integrating at some time N ini such that k/k i ≃ 100, and terminate the integration a few e-folds after the strong dissipative phase (in which Γ ≫ H) ends and the mode being computed satisfies k ≪ aH.
It is worth stressing that since the stochastic source ends up dominating the evolution of the perturbations, the choice of initial conditions is actually not too relevant.This will be made clearer in Section 3.3, but for now let us illustrate it with a numerical example, redefining Q i to model the deviation from the Bunch-Davies initial conditions by multiplying the original Q i of (3.14) by some  (real) number ε BD .The effect of varying ε BD with respect to the case ε BD = 1 is shown in Fig. 3.Even for very large values of this parameter, ε BD ≃ 10 6 , we find that within roughly 1 e-fold (and several e-folds before horizon crossing) the solutions converge to the same value.
As we mentioned earlier, we find that in some cases the system of differential equations is numerically more stable if we do not impose the constraint of Eq. (2.19).This gives rise to an additional equation of motion (for the variable δq r ).The matrices for this 5 × 5 system are presented in Appendix E. We have checked that the numerical results using either set of equations are in agreement.The power spectrum for the benchmark point of the previous section obtained by solving either system is shown as a solid line in Fig. 4. The evolution of the perturbations for the mode k peak at which the power spectrum peaks is shown in Fig. 3.

Stochastic equations
In principle, to determine the probability distribution for the stochastic variable Φ, one should solve the Fokker-Planck equation (3.8).This is a rather difficult task.An alternative consists in estimating numerically the probability distribution by using a frequentist approach, i.e. by solving the system of Langevin equations (3.4) over many different realizations, where dW r N ≡ √ 2Re(ξ N )dN and dW i N ≡ √ 2Im(ξ N )dN are realvalued, independent Wiener increments. 8This is the approach that we adopt in this section.
We impose the following initial conditions, in accordance with Eqs.(2.20) and (3.14) The limits of integration are discussed below Eq. (3.14).We solve the Langevin system with a fixed time-step Runge Kutta method. 9The convergence of the solution was checked by successively decreasing the time-step.We found that decreasing the time step below ∆N = 10 −4 produces results for the averaged primordial spectrum that are indistinguishable at the percent level.The curvature perturbation and the corresponding power spectrum are determined by substituting the solution of the Langevin system into Eq.(3.1).
Fig. 4 shows (as light blue dots) a collection of 2160 stochastic realizations of the power spectrum for twenty different values of the wavenumber k.The dark blue dots represent the arithmetic average of all the realizations for each k.The continuous black curve, in turn, corresponds to the numerical solution of the matrix equation (3.10).The bottom panel of this same figure shows the relative difference between the frequentist approach and the matrix formalism solution.The result is a stochastic average which agrees with the matrix formalism results at the percent level.
The Langevin method provides for us not only the means to determine the first moment of the probability distribution of the power spectrum, but with enough sampling we can also determine the full distribution for P R (k) at each value of k.The left panel of Fig. 5 shows the normalized histogram for the 2160 realizations for log 10 P R (k) at k = 0.8 k peak for illustration.In this same panel we show as a vertical dashed blue line the corresponding expectation value over realizations, and as the vertical red dashed line the mean computed using the matrix formalism (presented in Section 3.1).The continuous black curve corresponds to a skew-normal fit to the (logarithmic) data.As a reminder, a random variable x is skew-normal distributed if its probability distribution function is given by where erfc(x) denotes the complementary error function and {µ, σ, α} are free parameters.Therefore, we find that the PDF of P R can be modelled as a skew-log-normal distribution. 10Defining for each we find that its probability distribution is very well approximated by a k-independent skew-normal distribution.The right panel of Fig. 5 shows the frequentist histogram for the full set of realizations for (3.18).Together with it we show the corresponding universal skew-normal fit (shown in solid black), with parameters {µ, σ, α} = {0.42,0.87, −4.15}.
A similar histogram can be created separately for each k, and we find that the standard deviations of the parameters {µ, σ, α} for each one of these histograms with respect to the corresponding values for the universal fit shown above are of order {3%, 2%, 9%}.
Note that the variance of the probability distribution function for the power spectrum is quite large.From Fig. 5 it is clear that, for a specific realization in a particular Hubble patch the spectrum can reach a value roughly one order of magnitude away from the 10 −2 value required to obtain f PBH ≃ 1, leading to either overproduction or underproduction of PBHs (according to the Gaussian estimate of the abundance).This effect can always be countered by adjusting any of the parameters in Γ that control the overall size of the average of the power spectrum, as discussed in Section 2.1, as well as the threshold for the collapse (on which the abundance depends exponentially within the Gaussian estimate).

Analytical approximation
To get a better understanding of the evolution of the perturbations and the shape of the primordial spectrum, it is useful to simplify the equations of motion in such a way that they can be solved analytically.Let us begin by noting that at late times, the only quantity that contributes to the curvature perturbation is δϕ, The second observation we make is that, in the equation of motion for δϕ, (2.16), we can neglect several terms and still reproduce the most important features of the spectrum, This approximation is obtained by discarding terms involving the potential (which are slow-roll suppressed), the metric perturbation (which can be checked numerically to be a good approximation), and assuming that Γ ∝ T 3 and the background remains in the attractor at all times, so that the attractor parameters defined in Eq. (2.5) indeed satisfy ϵ ϕ = ϵ ρ = 1.This last approximation is justified by Fig. 1, where it can be seen that the background quantities only leave the attractor for very brief periods.In addition, we have found numerically that the stochasticity of the system can be encoded via δρ r in (3.21), and therefore the original noise term on the right hand side of (2.16) can be dropped.
Let us explain more precisely this last approximation.If we have an explicit expression for δρ r as a function of time, then we can think of the δρ r term in Eq. (2.16) as a source term for δϕ, on stochastically, parametrically excited spectator fields during inflation [21,22].
the same footing as the ξ N term.Numerically, we find that the δρ r term dominates over the ξ N term, in the sense that one can set the latter to zero and still correctly reproduce the key features of the spectrum: the location of the peak, and the size of the spectrum; the latter within an order of magnitude of the full numerical result.Note that this does not mean that the noise ξ N is irrelevant.In fact, it is precisely the noise in Eq. (2.17) which determines δρ r , and thus in turn δϕ.In other words, the power spectrum of the comoving curvature perturbation is enhanced thanks to the source δρ r (whose value is set by the thermal noise) in the equation of motion for δϕ.
The strategy we will follow now is to propose a phenomenological parameterization for δρ r as a function of time and use it to solve Eq. (3.21).We will also assume that all background quantities can be approximated as piecewise-constant functions.The benchmark values we take for each quantity are shown in (the table of) Fig. 6, where we introduce the quantity These parameters have been chosen in order to obtain a primordial spectrum closely resembling in its main features the one derived for the dissipation coefficient (2.21) with the parameters (2.22).We assume the evolution proceeds in four different phases, which we label from 0 to 3. In phases 0 and 3 we have Γ ϕ = 0 and Γ ≪ H, so that we are in the weak dissipative regime.In phases 1 and 2 we have Γ ≫ H.During phase 1 we have Γ ϕ > 0, and during phase 2 we have Γ ϕ < 0. This parameterization is compared to the benchmark model (2.22) in Fig. 6.In addition, we parameterize the time evolution of the root mean square of δρ r with the following phenomenological expression, where the time N f at which the transition occurs is located a couple of e-folds before the horizon crossing time, where ∆N f is an O(1), k-independent constant.We take ∆N f = 2.1 for definiteness.Despite its simplicity (which of course cannot capture all the details of the full numerical solution), this parameterization is enough to understand the most important features of the spectrum.Since δρ r is a stochastic variable, it is not sufficient to parameterize its root mean square value, but we also need to know its correlation function.To make progress, we will assume that δρ r behaves like a Wiener increment, where the correlation function for ξ δρ N is The homogeneous solution to Eq. (3.21) can be written as ).The corresponding approximations as piecewise-constant functions from the table on the right are shown with dashed lines.Right: Benchmark parameters for the analytical calculation of the power spectrum.We take phases 1 and 2 to end at N 1 = 2.5 and N 2 = 5, respectively (we measure the number of e-folds from the end of phase 0, so that N 0 = 0, and we normalize a(N 0 ) = 1, see the figure on the left, as well as H = 7 × 10 −6 M p .
where δϕ ± are constants fixed by the initial conditions, J µ is the Bessel function of the first kind and This solution and its derivative can also be written in matrix form as δϕ + δϕ − .
(3.29)The constants µ and ν take different values in each one of the four phases.We denote their values in the j-th phase by µ j and ν j .The constants δϕ ± can be found by imposing continuity of the solution and its derivative at the end of each phase.We denote these constants by δϕ ±j in the j-th phase.We use N j to refer to the time at which each phase ends.In particular, phase 0 begins at −∞ and ends at N 0 = 0, and phase 3 ends at N 3 = ∞.To be as general as possible we keep our calculations generic for n + 1 phases, but we will eventually set n = 3.
Following the above procedure we can find the constants in the last phase In this equation, terms with smaller j should be placed at the end of the product. 12The total solution for δϕ, including both the homogeneous and inhomogeneous solutions is where G is the Green's function, which we will find below, and where Θ is the Heaviside step function.The constants in the homogeneous solution are obtained by imposing Bunch-Davies boundary conditions in the 0-th region,13 The expectation value of the power spectrum at late times is14 where Γ E (1 − µ n ) denotes Euler's Gamma function evaluated at (1 − µ n ) and Σ was defined in Eq. (3.22).The Green's function G appearing in the integrand of (3.35) is G(N, N ) = δϕ (1) (N )δϕ (2) ( N ) − δϕ (1) ( N )δϕ (2) (N ) N δϕ (1) ( N )δϕ (2) ( N ) − d d N δϕ (2) ( N )δϕ (1) ( N ) where N < N and δϕ (1,2) are two linearly independent solutions to the homogeneous equation.The calculation of the Green's function is simpler if instead of writing the homogeneous solutions as linear combinations of J µ and J −µ , we use J µ and Y µ (the Bessel function of the second kind, which is itself a linear combination of J µ and J −µ ).We therefore write which is completely equivalent to Eq. (3.27).Since the Green's function is independent of the boundary conditions chosen for the two linearly independent solutions, we can follow a slightly different procedure from before and arbitrarily choose some linearly independent set of constants in the final region instead of the first.The constants in the previous regions can then be found by matching the solutions and their derivatives at each boundary.We choose δ φ(1) +n , δ φ(1) −n = (0, 1) and δ φ(2) +n , δ φ(2) −n = (1, 0) for the two solutions.Figure 7: Left: the stochastic average of the power spectrum using the analytical approach is shown as a black solid line.The grey dashed line shows the result obtained with the matrix formalism of Section 3.1, see also Fig. 4. Right: homogeneous (solid) and inhomogeneous (dashed) solutions -the two terms in Eq. (3.31)-for modes leaving the horizon at the start and the end of the strongly dissipative phase in which Γ ≫ H.The inhomogeneous solution, which is independent of initial conditions, always dominates at late times, indicating the presence of an attractor in the equation of motion for the perturbations.The parameters chosen for both panels are shown in (the table of ) Fig. 6.
The reason for using Y µ instead of J −µ and choosing the constants in the final region instead of the first is that we obtain the following simple limits for the two independent solutions at late times, δϕ (1) If N is in the i-th region, the denominator of the Green's function becomes It is easy to show that this combination of constants is Since in the final region we have δ φ(1) Putting everything together, we find that the Green's function at late times is The integral in Eq. (3.35) then reads where we have also switched variables to v = e −N .These integrals can be found analytically in terms of hypergeometric functions.Now that we have all the necessary ingredients we can compute the power spectrum analytically by using Eqs.(3.35) and (3.44) and fixing the parameters as in (the table of) Fig. 6.To go from δϕ to R we use Eq.(3.20) in the late time limit, where the ratio between the two is approximately constant, see the central panel of Fig. 1.The resulting power spectrum is shown in Fig. 7.The overall size of the peak of the spectrum and the oscillations seen in Fig. 4 are present in the analytical solution.The oscillatory pattern observed for scales k > k peak in the power spectrum is attributed to Bessel functions appearing in the matching conditions imposed on the homogenous solution of the δϕ equation in the four phases.We find that, as with the numerical solution, the peak in the spectrum occurs for modes that leave the horizon around the end of the strongly dissipative phase.This is a consequence of the enhancement being an integrated effect, due to Eq. (3.44), as opposed to a local one.
Having an analytical solution allows us to understand why the initial conditions for the perturbations are irrelevant.All of the information about initial conditions is contained in the homogeneous solution (3.27) inside the integration constants δϕ ± .However, as shown in the right panel of Fig. 7, this solution is completely negligible at late times.The spectrum is completely dominated by the integral in Eq. (3.44), which is independent of initial conditions.This indicates the presence of an attractor in the equation of motion for the perturbations, as anticipated earlier.
The analytical approximation developed in this section is not enough to reproduce with accuracy the full averaged primordial spectrum.For instance, the actual slope of the log of the spectrum for k < k peak is about twice the value that the analytical approximation gives.However, this approximation is going to be useful in the next section to obtain a good estimate of the peak value of the gravitational wave signal induced at second order in perturbation theory.

Induced gravitational waves
In this section we compute the gravitational wave signal induced by scalar perturbations at second order.These gravitational waves are induced both during inflation and during the subsequent radiation era [24,25].The calculation is organized as follows.In Section 4.1 we write the equation of motion for tensor modes sourced by second order scalar perturbations (obtained by perturbing Einstein's equations) and derive explicit expressions for the source term in the two cases of interest, namely, when gravitational waves are induced in the inflationary epoch, and during radiation domination.We then solve this equation via the Green's function method.In Section 4.2 we compute the tensor power spectrum and in Section 4.3 relate it to the observable quantity of interest, which is the energy density of gravitational waves.In Section 4.5 we provide a numerical estimate of the latter quantity.

Second order scalar source
The equation of motion for tensor modes at second order in Fourier space and in terms of conformal time η is [24] h where primes denote derivatives with respect to conformal time (′ = d/dη), H = a ′ /a denotes the conformal Hubble factor, and the superscript s = (+, ×) stands for the polarization of the tensor mode.The source S s k is, in the Newtonian gauge and in the absence of anisotropic stress [26], The quantity e s ij is a symmetric transverse, traceless projector that satisfies e s ij (k)k i = 0 [26].The second term in the equation above can be alternatively written in terms of the total momentum perturbation15 δq using Eq.(2.11), In particular, during inflation, the total momentum perturbation is the sum of the radiation and inflaton components, The source term S s k of Eq. (4.2) before (pre) and after (post) inflation ends is respectively given by: where e s (k, p) ≡ e s ij (k)p i p j .To write the first expression we have neglected the first term in the brackets of Eq. (4.2), as well as the δq r term from Eq. (4.4).This approximation will be justified in Section 4.5.For the second expression we have assumed that the Universe enters a radiationdominated era after inflation ends, so that p = ρ/3.
As is customary, let us fix the time at which inflation ends as η = 0.The value of ψ at this time, which will be the initial condition for the post-inflationary source, is, on superhorizon scales and assuming the Universe enters a radiation-dominated era after inflation ends, where we have used Eq.(3.20).The post-inflationary source can then be rewritten as where the function Q is and T ψ k is the transfer function for ψ k , which is defined by and is obtained by solving16 The Green's functions for Eq.(4.1) during inflation and during the subsequent radiation-dominated era can be found as in Eq. (3.36).The results are, respectively, ) where we have used H = −1/η during inflation.The solution is therefore where T h k is the (linear) transfer function of h s k in the radiation era, and The lower integration limit η h is some early time at which we assume h s k (η h ) = 0.In other words, we assume that no second order gravitational waves have been induced at sufficiently early times.

The gravitational wave spectrum induced at second order
We now have almost all 17 the necessary ingredients to calculate the energy density of gravitational waves, which is related to the tensor power spectrum and is the observable quantity of interest.The expectation value of the tensor power spectrum late in the radiation era contains three terms, 18 (4.17) These three terms will in turn lead to three different contributions to the gravitational wave energy density.The first term corresponds to the gravitational waves induced during the inflationary epoch, and the third term corresponds to the gravitational waves induced during the subsequent radiationdominated era.The middle term mixes both contributions and its value typically lies between the other two.
To perform the rest of the calculation we will use the analytical results of Section 3.3.The reason for this is that to calculate the tensor power spectrum we need to take the quantum expectation value in addition to the stochastic one, in order to make the tensor power spectrum a deterministic quantity, as we did with P R .As discussed in Appendix D, when Eq. (3.31) is used to split δϕ into the homogeneous and inhomogeneous solutions to its equation of motion, finding the double expectation value is straightforward, since only the homogeneous solution is quantized, and only the inhomogeneous piece is stochastic.This splitting can be done only because the equation of motion for δϕ has been decoupled from the rest (up to the source term δρ r ), since the full system of differential equations cannot be solved by Green's function methods.Our simplifying approach should give a reasonable estimate of the actual result.
The first term in Eq. (4.17) is defined by 19 The quantum expectation value is already included inside P pre on the right-hand side -see Eq. (D.4) for the analogous scalar definition-so we only write the stochastic average.The two-point function on the left-hand side can be computed quantizing the inflaton field perturbation.The result is, using Eqs.(4.5) and (4.16), The four-point function of δϕ appearing in this equation can be computed using Eq.(3.31), together with Wick's theorem.This calculation is done in Appendix D. The resulting expression can be 17 The connection between the power spectrum of tensor modes and the energy density of gravitational waves will be completed in Section 4.3. 18In Section 4.4 we present for the first time a compact expression which includes the mixing term in standard (cold) inflation. 19As discussed in Appendix D, the brackets without subscripts denote a double expectation value, quantum and stochastic.
substituted back into the above equation and, using one of the Dirac deltas to perform the integral over l, we find where we have defined20 This quantity has the following properties It is useful to make the integrand in (4.20) manifestly real and symmetric under η ′ ↔ η ′′ .To do so, we take Eq.(4.19), rename the dummy variables η ′ ↔ η ′′ and sum the result with Eq. (4.19) itself.After using the second identity in Eq. (4.22), we find The final step consists in switching to spherical coordinates in the q integral, perform one of the angular integrals, and then make the change of variables This procedure amounts to the replacement [27, 28] Thus, we finally obtain for the first term in Eq. (4.17), To compute the term 2⟨P mix (k, η)⟩ S , we use where g k (η) is defined in Eq. (4.14), and similarly for the P post term.Following a completely analogous procedure to the one we just applied to compute ⟨P pre (k, η)⟩ S , we obtain and similarly, where we have used Eqs.(4.7) and (4.22) to relate Q δϕ to P R .This completes the calculation of the tensor power spectrum late into the radiation era.

The energy density of induced gravitational waves
The stochastic average of the energy density of gravitational waves is [29,30] where the bar denotes a time average over many wavelengths.The value of this quantity at some late time in the radiation era (when the temperature is T ) can be related to its value today (at temperature T 0 ) assuming entropy conservation [28], where g ⋆s is the effective number of entropic degrees of freedom.The time average can be obtained using The H 2 factor in the denominator of ⟨Ω GW ⟩ S will cancel out the 1/η 2 in these averages, yielding a finite result in the limit η → ∞. 21The energy density of gravitational waves today is therefore 2 K(ky, kz), (4.35) 21 The upper integration limit should be set to today, but since the scalar source decays quickly during the radiation era, the difference in the result using either limit is negligible [28].This choice simplifies the integrals in Eqs.(4.40, 4.41).
We have also neglected the evolution of the tensor modes during the late matter-dominated era, since the corresponding corrections are highly suppressed for the frequencies we are interested in, see e.g.[25].

Induced gravitational waves in the noiseless limit
For completeness, we present here the expression for the tensor power spectrum valid in the standard cold inflation case; that is, in the absence of the second term in the parentheses of Eq. (4.21).In this case the spectrum can be written as the square of a sum, where

.44)
The first term inside the square gives the gravitational waves induced during inflation [25] and the second one those produced during radiation domination.The mixing term obtained by developing the square has not been presented before in the literature.The gravitational wave energy density can then be found by using Eq.(4.31) and computing the time average with Eqs.(4.32), (4.33, (4.34).This result is thus useful for computing the full induced (inflation plus radiation) second order gravitational waves in standard (non-dissipative) inflation.

Approximating the time integrals
In this section we estimate the inflationary and post-inflationary contributions to the energy density of the induced gravitational waves in our scenario with a transient dissipative phase.Let us focus first on the K pre term of Eq. (4.36).This contribution depends on the lower integration limit η h (and as noted in [25] is formally divergent in the limit η h = −∞).We deal with this problem by integrating from a finite value of η h that we identify with the time at which the strongly dissipative phase begins.The assumption here is that the contribution from the source prior to this time is negligible.This assumption is reasonable since up to that time inflation proceeds as in the standard slow roll scenario (up to the presence of a weak dissipative term that does not alter the dynamics significantly), and we do not expect the corresponding gravitational wave signal to be peaked at any particular scale or exhibit any special features, in contrast to the piece arising due to the strongly dissipative phase.
In addition, we notice that the inflationary contribution to the energy density of gravitational waves diverges in the y = ∞ (infinite comoving wave number) limit.In principle, this divergence should be renormalized away by properly computing the induced gravitational wave signal using the in-in formalism.However, we just impose a cut-off which renders the result finite.
We have verified that our results do not depend on the cutoffs in conformal time and wave number unless unreasonably large values are chosen.We remark that only the K pre and K mix kernels can be affected by these ambiguities, unlike the post-inflationary contribution (which is finite).We also stress that the results of this section should be regarded as an accurate order of magnitude estimate of the overall size of the signal.
We can choose the time cutoff around the time at which the dissipative coefficient Γ begins to increase, which corresponds to the departure from standard cold inflation.In terms of the analytical calculation of Section 3.3, this corresponds to the beginning of phase 1.For the momentum cutoff we can choose k cutoff ∼ O(10−100)×k peak .The four-dimensional integral in ⟨Ω pre ⟩ S is quite difficult to perform.Fortunately, K pre is heavily peaked around a specific time, so the strategy we adopt is to approximate the time integrals by evaluating the integrand at this time and multiplying it by an appropriately chosen integration area.
To determine the point in parameter space at which the integrand is peaked, we use the fact that the integrand is symmetric under η ′ ↔ η ′′ and y ↔ z, so the set of maxima of the function must be symmetric under this transformation.If the function has a unique global maximum in some region (we do not prove that this is the case, but we have checked it numerically), then it follows that this maximum must be located along the surface with η ′ = η ′′ and y = z.On this surface where η max is the value of η at which the local maximum occurs, and ∆η max is the integration area, which must be appropriately chosen as a small square around η max requiring, for instance, that the integrand does not decrease by more than an order of magnitude or so, in such a way that the approximation holds.The integration area is, in terms of the number of e-folds, The function we need to maximize is therefore22 kF pre (0, η max ) 2 P δϕ (ky, η max )P δϕ (kz, η max ).(4.48) The quantity in Eq. (4.48) is shown in the left panel of Fig. 8 for the parameters in (the table of) Fig. 6.The discontinuity around N max = 5 e-folds is due to the fact that, as mentioned in Section 3.3, we take the background parameters as piecewise-constant functions for this calculation.Specifically, we take ϕ ′2 /a 2 (ρ+p) = 1 in phases 0 and 3, and ϕ ′2 /a 2 (ρ+p) = 0.02 in phases 1 and 2. In this figure we also take ∆N max = 2, which is clearly enough to account for the region in which the integrand is large.Changing this number by a factor of O(1) does not change our results.
To obtain the induced gravitational wave signal, we find the time η at which K pre is peaked for each k and perform the momentum integrals numerically.The time-dependent power spectrum P δϕ (k, η) Fig. 8.However, since the largest contribution to the integral comes from the region around ky = k peak , to make the calculation numerically less demanding we can simply take Nmax as the value at which the integrand, evaluated at ky = k peak , is peaked, and use the same value Nmax for all y.We have explicitly checked that the peak of the signal remains unchanged if the y-dependence of Nmax is taken into account.is calculated using the analytical formalism of Section 3.3, as we anticipated earlier.Specifically, it can be found by keeping the full expression for the Green's function in Eq. (3.36) instead of taking the N → ∞ limit.The resulting signal is shown in the right panel of Fig. 8.We find that the energy density of gravitational waves induced during inflation is much smaller than that of the gravitational waves induced during the radiation era.We do not show in this figure the mixed term from Eq. (4.36), but we find that it is well approximated by ⟨Ω mix ⟩ S ∼ ⟨Ω post ⟩ S ⟨Ω pre ⟩ S , and is therefore also suppressed with respect to the post-inflationary contribution.
We stress that approximating the integrand by its peak value is what allows us to neglect the subdominant terms involving ψ and δq r in Eq. (4.2).Let us show that this is a good approximation by estimating the relative contribution of each term in this equation during the inflationary epoch.Since the integrand in Eq. (4.2) is symmetric under p ↔ k − p, for the purpose of finding out which terms contribute the most at the point at which this integrand reaches its largest value, we can simply evaluate it at |p| = |k − p| as we did for K pre .We can then define where we have ignored the mixed terms in Eq. (4.2), since they are always subdominant.This quantity is shown in Fig. 9 for two different Fourier modes, both of which become super-Hubble near the end of the strongly dissipative phase.The figure illustrates that the time integral of the scalar source of noise is dominated by the δϕ term.Since the other contributions are at most of the same order and will therefore not change the size of the peak in the gravitational wave signal, we can neglect them.
Let us briefly summarize the results of this section.We have made three different approximations in this calculation.The first is that we approximated the time integrals as their peak value times an appropriately chosen area.The second is that we have neglected the contribution of the radiation and metric terms to the gravitational wave source during the inflationary epoch, and we have checked explicitly that the contribution of these terms is indeed subdominant.These two approximations are very good and should not change the order of magnitude of the result.If the time integrals are performed numerically and the metric and radiation perturbations are included in the source term, we expect that the size of the signal will change at most by a factor of O(1).The third and final approximation is related to the divergences in the far past and for large momenta due to the behaviour of scalar modes in the Bunch-Davies vacuum.We have assumed that this effect can be taken into account correctly by imposing reasonable cutoffs.We have checked that our results are independent of the cutoffs unless an unreasonably large choice is made.The uncertainty due to this approximation is difficult to quantify, but our results should be correct as an order of magnitude estimate.

Conclusions
We have shown that a period of dissipation lasting a few e-folds during inflation can lead to a large and peaked primordial spectrum of curvature perturbations P R at specific scales.This may provide the seeds of an abundant PBH population capable of accounting for all the dark matter of the Universe.The large power spectrum that is required in the most naive estimates (P R ∼ 10 −2 ) leads to a stochastic background of gravitational waves that we have computed using second order cosmological perturbation theory.As it is well known, given the current astrophysical bounds on the PBH abundance, this background of gravitational waves falls within the sensitivity reach of the future LISA interferometer.
The growth of P R is due to the stochastic nature of the dissipation.Thermal fluctuations of the radiation background that originates from the inflaton transferring its kinetic energy to lighter degrees of freedom act as a stochastic source for the curvature perturbation.This makes the mechanism very different from those operating in other inflationary setups that can produce abundant PBH dark matter (in particular, ultra slow-roll from an inflection point in single-field inflation).
We have described a method for computing the curvature power spectrum that allows to overcome the numerical difficulties associated to the stochastic character of the system of equations governing the time evolution of the perturbations.The method consists in formulating the problem as a system of deterministic differential equations for the two-point correlation functions in Fourier space.We have verified the validity and accuracy of the method by solving the stochastic equations for many realizations and computing statistical averages.The method we have presented can be useful in more general contexts in cosmological perturbation theory where stochastic sources are present, beyond the specific scenario we have discussed in this work.We have also shown that a single, simplified, stochastic differential equation with an analytical solution can be used to understand qualitatively the enhancement of the primordial spectrum.
Our analysis of the stochastic background of gravitational waves includes two aspects that make it richer than that of previous works.First, there is the unavoidable stochastic origin of the primordial curvature spectrum, which propagates into the calculation of tensor modes at second order and must be appropriately taken into account.In addition, we have computed not only the gravitational waves from modes entering the horizon during radiation domination, but also those generated at second order during inflation, and we have written an explicit expression for the term that describes the mixing between the two contributions.As a particular case, we also write for the first time an expression including the mixing term in standard (non-dissipative) inflation.These formulas could be useful in contexts beyond the scenario presented here.Applying the analytical approximation that we derive for the primordial curvature spectrum, we estimate the gravitational wave signal and find that it is dominated by the gravitational waves induced during the radiation-dominated era that follows inflation.
We have used a phenomenological parameterization to model the transfer of energy from the inflaton to the radiation background.Although there have been serious efforts on inflationary model-building aimed at obtaining significant dissipation throughout inflation (specifically, within the framework of warm inflation), we are not aware of any previous works proposing a localized dissipative period akin to the one we have explored.We have discussed, in Appendix A, a potential route to describe the microphysical origin of a peaked dissipative coefficient such as the one we propose by starting from a particular Lagrangian, but further work is needed to find concrete, well-motivated models.
The stochastic nature of the peak in the primordial spectrum is a general property of our scenario.We find that the spectrum presents a series of oscillations after the peak, a feature that is not commonly found in inflection-point models of PBH production.Similarly, the dip that is commonly seen before the peak of the spectrum in this class of models is absent in our scenario.Other details of the spectral shape and, consequently, the gravitational wave signal, are model dependent.It would be interesting to investigate whether there may be any features that could help distinguish our mechanism from other models should PBHs and a stochastic background of gravitational waves be detected in the adequate masses and frequencies of interest for dark matter.

A Microphysics of the dissipative coefficient
In this appendix we discuss a particular microphysical realization of a localized dissipative coefficient Γ during inflation.The purpose of the present discussion, however, is not to propose a definitive model but simply to show that obtaining a peaked dissipative coefficient that satisfies all the necessary constraints is, in principle, possible.Specifically, we introduce a Lagrangian which reproduces the form of the coefficient (2.21), assuming that the fields that couple to the inflaton are part of a thermalized bath.In order to keep the field content to a minimum we will consider a scenario where, besides the inflaton ϕ, only two additional degrees of freedom participate in the dynamics.The first one, denoted by σ and assumed to be a scalar, corresponds to a light radiation field in equilibrium.The second one, also a scalar and denoted by χ, corresponds to a heavy catalyst field.The large effective mass of χ arises via its coupling to the slowly rolling inflaton field with a non-zero vev.Through the coupling of the (unstable) χ with σ, the inflaton energy density can be efficiently dissipated into radiation.Indirect decay scenarios like this one are among the preferred mechanisms for realistic warm inflation models [32][33][34].One of the advantages of introducing a heavy catalyst field is to prevent the inflaton potential from receiving strong temperature corrections [34,35]. 23onsider the Lagrangian Here φ denotes a non-canonically normalized inflaton, due to the presence of the function K(φ) in its kinetic term.In this basis, the inflaton and the mediator χ interact via a four-legged contact term with coupling g.On the other hand, the term with coupling g connects the three fields.For a non-vanishing inflaton vev, this term induces the decay of χ into σ.The ellipsis corresponds to the interactions which are necessarily present in order to thermalize the light sector σ.If this sector is indeed in equilibrium, the presence of the ϕ → χ → σ channel modifies the equation of motion of the inflaton.In terms of the canonically normalized inflaton ϕ, related to φ via dϕ dφ = 1 the effective equation of motion for ϕ can be written as where the potential V (ϕ, T ) includes the thermal correction due to the presence of a bath of σ particles (fluctuation), while Γ encodes the production of σ quanta (dissipation).The appearance of a local dissipative term relies on the assumption that the microphysical processes which determine Γ operate at time-scales much smaller than those characteristic of the macroscopic slow roll of the inflaton and the expansion of the Universe (the adiabatic-Markovian approximation [32,36]).Additionally, we assume that the typical interaction time-scale between the constituents of the thermal bath is much shorter than the time-scales associated to the variation of the background quantities.In this approximation, the dissipative coefficient can be evaluated as [37][38][39] Γ(ϕ, where ω denotes the 0-th component of the 4-momentum p, n χ (ω) = (e ω/T −1) −1 is the Bose-Einstein distribution, and is the spectral density of the catalyst field χ.In this expression ω 2 p = |p| 2 + m 2 χ is the on-shell frequency and Γ χ is the decay width of χ.Note that our expression for Γ originates from the ϕ − χ interaction term in (A.1).In principle, the term proportional to g2 would also contribute to the dissipation rate, although this contribution is loop-suppressed [39].
In order to evaluate the dissipation coefficient, we need the decay width Γ χ .In general, this width must be computed in non-zero temperature QFT, and a general expression can be found in e.g.[39].For simplicity, we will restrict our discussion to the so-called low-temperature limit in what follows, which corresponds to assuming that T ≪ m χ .We remark however that in principle none of our assumptions forbid a peaked dissipative coefficient at higher temperatures.In the low-temperature limit the decay rate for the process χ → σσ may be written as in a frame boosted with respect to the rest frame of χ.In terms of this rate, the strong dissipative condition, necessary to maintain thermal equilibrium between the light σ and the heavy χ, and the adiabaticity requirement correspond to Moreover, in this low-temperature regime, the thermal correction to the inflaton potential can be neglected, V (ϕ, T ) ≃ Ṽ (φ(ϕ)) [34,35].Introducing the dimensionless quantities where Γ χ,0 = Γ χ (m χ , 0), and m 2 χ = g 2 φ 2 , c.f. (A.1), the dissipation rate (A.4) can be written as At low temperatures τ ≪ 1, the first term of the integrand can be simplified disregarding the u and γ χ -dependent terms in the denominator (equivalent to approximating the spectral function as ).The dissipation rate then simplifies to where we have used the fact that the p-integral can be well-approximated by 5τ 3 /2 for τ ≪ 1.
In addition, λ = 4.6 × 10 −14 , g ⋆ = 10 in (2.1), (2.23).The left panel shows the time evolution of the rate Γ, together with the expansion rate H and the temperature T in terms of the number of e-folds N .The right panel shows our consistency checks for the conditions in Eq. (A.7), as well as the small temperature condition T /m χ ≪ 1.The self-consistency of our set-up is therefore guaranteed.As a final remark, we note that we assume that the potential V (ϕ) = Ṽ [φ(ϕ)] supports slow-roll inflation, and in our particular example it is a quartic polynomial (2.23).Owing to the relation between the canonical and non-canonical forms of the inflaton field, given by (A.12), Ṽ (φ) would then be a relatively complicated function of φ.In any case, as we discuss in Section 2.1, the precise shape of the inflaton potential does not affect our conclusions.This is true as long as the potential does not possess any features that interfere with the effect of the dissipation rate Γ, which is the main quantity that determines the growth of fluctuations.We expect that the discussion presented in this appendix encourages further efforts to search for inflation scenarios which lead to peaked dissipative coefficients.

B Temperature-independent dissipative coefficient
In this appendix, we investigate the possibility that the dissipation coefficient is temperatureindependent and compute the power spectrum using the three methods presented previously.In order to assess the effects of the temperature dependence of the dissipative coefficient on the power spectrum, we consider a coefficient with the same dependence on the inflaton field as in Eq. (2.21): As a benchmark example, we choose the following set of parameters: (B.2) Given the fact that the dissipative coefficient is temperature-independent, in order to compute the curvature power spectrum one can directly apply the matrix formalism approach of Sec.3.1 or solve the system of Langevin equations as in Sec.3.2 by setting terms proportional to Γ T to zero.The power spectrum obtained using these two approaches in represented in Fig. 11 showing a good agreement between averages over stochastic realizations and the matrix formalism approach, typically at the percent level around the peak and of the order of 10 − 20% away from the peak.
In the analytical approach of Sec.3.3, we argued that the fluctuations δϕ were driven by the stochastic noise via the δρ r term.However, this term is proportional to Γ T and vanishes for the temperature-independent dissipative coefficient.Therefore, the driving stochastic term in the equation for the fluctuation δϕ is the explicit stochastic noise term, appearing on the right-hand-side of Eq. (2.16).By setting Γ T = 0, the approximate equation for δϕ can be written as after neglecting metric perturbations and slow-roll-suppressed terms.The correlation function for The solid black line represents the average of the power spectrum obtained via the deterministic matrix differential equation and the grey dashed line the results from the analytical approach.Bottom: Absolute value of the relative difference between the stochastic average and the matrix average of the power spectrum.The agreement for each k is at the percent level around the peak and of order 10 − 20% away from the peak.
which makes Eq. (B.3) formally identical to Eq. (3.21) upon substituting ξ N → ξ δρ N and rescaling the prefactor of the right-hand-side of Eq. (B.3) by the appropriate background quantities.Following the approach of Sec.3.3, we consider four time intervals.In the case of Eq. (B.3), the parameters s and S introduced in Eq. (3.32) and Eq.(3.33) become simply s = 3 and S = 1 throughout the four phases.The dissipative coefficient and relevant background quantities used for the analytical computation are listed in Table 1 for the four different phases.The resulting power spectrum is depicted in dashed-grey in Fig. 11.The analytical spectrum reproduces well the results from the matrix formalism approach at small k/k peak ≪ 1 and around the peak.The main differences can be seen in the oscillation pattern at larger k/k peak ≳ 1.The analytical approximation matches the numerical result significantly better than in the model with temperature-dependent Γ.The reason is that the calculation in Sec.3.3 involves more approximations, mainly neglecting ξ N with respect to δρ r and parameterizing δρ r by Eq. (3.23).
We conclude that the analytical approach developed in Sec.3.3 is consistent with both the matrix formalism and the stochastic approach presented in this paper.This approach allows to describe Table 1: Benchmark parameters for the analytical calculation of the power spectrum.We take phases 1 and 2 to end at N 1 = 2 and N 2 = 4, respectively (we measure the number of e-folds from the end of phase 0, so that N 0 = 0, and we normalize a(N 0 ) = 1 as well as H = 7 × 10 −6 M p .
and characterize the primordial power spectrum also for a temperature-independent dissipation rate.
We have shown that a temperature-dependent dissipative rate is actually not necessary in order to achieve an enhanced power spectrum.

C Stochastic differential equations
In this appendix we review some basic elements of stochastic differential equations, including the definition of a Wiener process and the derivation of the Fokker-Planck equation.Part of the discussion follows the presentation in [40].

C.1 Stochastic calculus
Let us consider a discrete-time equation ∆y = f (y, t)∆t + g(y, t)∆W t , (C.1) where ∆t are the time increments between consecutive time nodes.In other words, where T is the time range over which we consider the equation, and N is the number of time nodes in which said time range is divided.The increment ∆W t (known as Wiener increment) is randomly drawn from a Gaussian distribution at each time step, The ∆W t increments at every time step are independent from each other.The different solutions to Eq. (C.1) obtained via some finite sequence of random increments ∆W t are referred to as realizations.
Formally taking the continuum limit N → ∞ turns (C.1) into a stochastic differential equation dy = f (y, t)dt + g(y, t)dW t . (C.4) We will discuss this limit later on.For the moment, let us continue with non-infinitesimal time increments.
An important fact is that the only way for the distribution in Eq. (C.3) to be well-defined in the context of a problem like (C.1) is if σ 2 = ∆t.Let us illustrate this with a heuristic argument in a simplified case.If we consider the equation ∆y = ∆W t , a particular realization can be obtained by summing each random increment, y = N i ∆W t(i) .The variance of y, denoted by V (y), is simply the sum of the variance of each increment ∆W t(i) , since they are assumed to be independent.Let us impose σ 2 is equal to some power of ∆t, Then, the variance of y is This implies a consistency condition in the continuum limit.Indeed, taking N → ∞, V (y) vanishes for n > 1, and it diverges for n < 1.Thus, the only acceptable value is n = 1.Let us return to Eq. (C.1) and compute the expectation value of ∆W 2 t , where the expectation values ⟨• • • ⟩ are taken with respect to the distribution in Eq. (C.3), and we have used the fact that ⟨∆W t ⟩ = 0.By using Eq.(C.3), we can also show that the variance of the sum of the ∆W 2 t(i) is Since this quantity vanishes as N → ∞, we see that the sum of all the ∆W 2 t(i) in the continuum limit is not actually random, but deterministic, and must therefore be equal to its mean.Let us write this explicitly.In the continuum limit, sums of increments are promoted to integrals of differentials.We thus have We can now take the continuum limit of this equation.The ratio ∆z/∆t becomes the derivative dz dt .In the third term of the right-hand side, ∆W This is a deterministic equation and therefore standard methods for partial differential equations may be applied to it.

C.2 Correlators in Fourier space
Let us assume that the stochastic variable y depends not only on the time t but also on the spatial coordinates x.As mentioned earlier, the increments defining Wiener processes are often written as dW t = ξ t dt.The quantity ξ t can be thought of as a distribution, with correlation function given by The noise correlators in Fourier space are then where we have defined the shorthand δ 3 ± = δ 3 (k ± k ′ ).Following a similar procedure for the other entries, the entire correlation matrix can be computed For the real and imaginary parts of the noise Re(ξ t ) ≡ ξ r t and Im(ξ t ) ≡ ξ i t , we find In computing these entries it is necessary to use the parity of the δ-function, δ 3 (k) = δ 3 (−k).We therefore find that the real and imaginary parts of the noise are uncorrelated, but ξ t and ξ ⋆ t are not.

C.3 Derivation of the matrix equation
Throughout the rest of this appendix we will focus on a specific multivariate version of Eq. (3.4) which is close to the type of equations we deal with in our scenario.Let us consider24 where the stochastic time-dependent variable Φ(t) is an n-dimensional column complex vector, A is an n × n real matrix, σ is an n × m complex matrix, and ξ t is an m-dimensional column complex noise vector whose real and imaginary parts have correlation functions given by the multivariate analogue of Eq. (C.26).We have absorbed the overall 1/2 factor from Eq. (C.26) into the definition of ξ t , leading to the 1/ √ 2 factor in Eq. (C.27).We want to find the probability distribution P (Φ, Φ ⋆ , t), where Φ ⋆ obeys the equation of motion The noise vectors ξ t and ξ ⋆ t are not independent as per Eq.(C.25).In order to use a generalized multivariate version of the Fokker-Planck equation (C.22) for the probability distribution P (Φ, Φ ⋆ , t), we need to rewrite Eq. (C.27) in terms of uncorrelated noise sources.To this end, it is convenient to define the vectors Ψ ≡ (Φ T , Φ † ) T and χ where where α and Σ do not depend on Ψ.This notation might make α and Σ look like they have the same shapes, but α is a 2n × 2n matrix, whereas Σ is a 2n × 2m matrix (the 0 matrices inside have different shapes).The final step is to write this equation with correlated noises χ c (hence the subscript c) in terms of the uncorrelated noises χ u ≡ (ξ rT t , ξ iT t ) T .We have dΨ Since the noise sources in Eq. (C.31) are uncorrelated, Eq. (C.22) can be generalized directly, The sums in Eq. (C.33) can be expanded in terms of Φ and Φ ⋆ instead of Ψ.We find where the sum over each of k and l runs half the range than in (C.33).Finally, the equation of motion for the two-point statistical moments Q ≡ ⟨ΦΦ † ⟩ S , defined via

D Stochastic and quantum expectation values
As shown in Section 3.3, the solution to the simplified equation of motion for δϕ (3.21) is of the form where ξ N defines a white noise process satisfying (C.25) and S k is a function of time, the form of which is not relevant in what follows.We can quantize the field (δϕ k → δ φk ) by writing the homogeneous solution in terms of creation and annihilation operators Since we do not quantize the noise, we make the second term in Eq. (D.1) proportional to the identity operator.The quantum expectation value is defined via ⟨ , where |0⟩ defines the vacuum of the system, so we have The power spectrum is defined via We can show from Eq. (D.2) that ⟨δ φ(h) k δ φ(h) q ⟩ Q = (2π) 3 |δϕ (h) k | 2 δ 3 (k + q), (D.5) where 25

E Numerically stable system of equations
As discussed in Section 3.1, we find that ignoring the constraint in Eq. (2.19) and keeping Eq. (2.18) for δq r = 4ρ r δv r /3 is often numerically more stable.Here we present the matrices for the corresponding 5 × 5 system in the language of Section 3.1.We have checked that this system of equations yields the same results as the 4 × 4 system.Let Φ = (δρ r , δq r , ψ, δϕ, dδϕ/dN ) T .One should be careful with the order of the variables when comparing with the results of Section 3.1.The matrix A is dϕ dN , (E.5) 25 We have ignored some contact interaction terms that are not proportional to δ 3 (k + p).The upper integration limit arises because, if a < c < b, then Finally, we assume δq r = 0 initially, so the initial conditions matrix . (E.7)

F Fluctuation-dissipation theorem
A rigorous derivation of the stochastic term in Eq. (2.15) should be worked out in the context of the Schwinger-Keldysh formalism of non-equilibrium QFT.However, it is possible to gain some intuition on it from a classical toy model [41,42] in which a massive particle P moving along a coordinate q with momentum p is coupled to a system of N ≫ 1 harmonic oscillators with coordinates x i , momenta p i and natural frequencies ω i through coupling constants g i .The Hamiltonian of the system is particle Hamiltonian The last term of the Hamiltonian has been introduced ad hoc to simplify forthcoming calculations.It can be physically interpreted as a correction to the potential V (q) which ensures that the Hamiltonian can be written in a translationally invariant form when V = 0 (although these details are not relevant for us).The corresponding equations of motion are26 Should we know the initial conditions for both q, p and x i , p i , we could, in principle, integrate the full system as a mechanical problem by solving (F.2) and (F.3).However, since N ≫ 1, a more practical approach is a statistical one, in which we know q(0), p(0) and only some general properties of x i (0), p i (0).This suggests an interpretation of the harmonic oscillators as an environment with which the particle of interest P interacts.With this approach in mind, let us look for an effective description of the motion of P .We start by formally solving (F.3).The solution of the homogeneous equation is spanned by the functions sin(ω i x) and cos(ω i x).The solution to the full equation can be obtained computing a Green's function.Doing so, we get x i = x i (0) cos(ω i t) + p i (0) m i ω i sin(ω i t) + g i m i ω i t 0 ds sin[ω i (t − s)]q(s) .(F.4) Substituting this solution into Eq.(F.2) and integrating by parts we obtain M q + V ′ (q) + t 0 ds γ(t − s) q(s) = ξ(t), (F. i q(0) cos(ω i t) + p i (0) m i ω i sin(ω i t) . (F.7) The function γ(t) comes from the primitive of the sine inside the integral when integrating by parts.
The term with q(0) inside ξ(t) comes from the boundary term of the integration by parts evaluated in the lower integration limit.Finally, the boundary term evaluated in the upper limit exactly cancels the third term in (F.2).Notice that γ only depends on the intrinsic properties of the environment (i.e. the natural frequencies of the oscillators, ω i ), and the interaction between the particle P and the environment (i.e. the coupling constants g i ).We identify this function with the dissipative coefficient, since it appears multiplying the velocity of P in the equation of motion (even if it is under the integral sign).On the other hand, ξ depends on the initial conditions of the environment, whose only available description is, by fiat, statistical, so that x i (0) and p i (0) can be treated as random variables.Hence, we can identify ξ with a stochastic force.This simple example illustrates how deterministic dissipation (γ) and random fluctuations (ξ) arise together when coupling a system to a complicated environment.Furthermore, this model allows to extract information about the power spectrum of the noise ξ.Let us assume that the random variables x i (0) and p i (0) are homogeneously distributed along the classical orbits in phase space, and that the initial conditions of different oscillators are independent.Under these conditions, the classical virial theorem ensures that m i ω 2 i ⟨x i (0)x j (0)⟩ = 1 m i ⟨p i (0)p j (0)⟩ = δ ij ⟨E i ⟩, ⟨x i (0)p j (0)⟩ = 0, (F.8) where ⟨E i ⟩ is the expectation value of the energy of the i-th harmonic oscillator at t = 0. Furthermore, assuming the environment of oscillators is at thermal equilibrium at t = 0, the equipartition theorem implies that ⟨E i ⟩ = T .(F.9) Let us now compute the correlation of ξ(t) with ξ(s).From (F.7), we have ⟨ξ(t)ξ(s)⟩ = i,j g i g j ⟨x i (0)x j (0)⟩ cos(ω i t) cos(ω j s) + g i m i ω i g j m j ω j ⟨p i (0)p j (0)⟩ sin(ω i t) sin(ω j s).(F.10) The rest of the terms are zero either because they involve the average of one single random variable or because they involve the correlation of position and momentum.Substituting the correlators of initial conditions by their values, summing over one index and using (F.9), we get ⟨ξ(t)ξ(s)⟩ = T This illustrates the main idea of the fluctuation-dissipation theorem: the power of random fluctuations ⟨ξ(t)ξ(t ′ )⟩ associated to a dissipative process increases with the strength of the deterministic dissipative coefficient (and the temperature of the system).In the N → ∞ limit, it is common to write γ as where J(ω) is the so-called spectral density of modes of the environment.Let us consider the lowest order in a Taylor expansion of the spectral density, J(ω) = Γ 0 ω, with Γ 0 a constant.In this case γ(t) = 2Γ 0 δ(t).For this form of the dissipative coefficient, (F.5) reduces to M q + V ′ (q) + Γ 0 q(t) = ξ(t), (F.14) in which Γ 0 appears explicitly as a friction term.Similarly, the correlation in (F.12) yields ⟨ξ(t)ξ(s)⟩ = Γ 0 T δ(t − s).(F.15)

Figure 2 :
Figure 2: Stochastic average of the power spectrum computed for the benchmark parameters of Eq. (2.22) by varying g ⋆ (left panel) and m in units of M p (right panel).The horizontal axis is normalized at the scale k = k peak at which the peak in the benchmark spectrum (with g ⋆ = 8 and m = 1.4 × 10 −4 M p ) occurs.

Figure 3 :
Figure 3: Left: Time evolution of the averaged power spectrum of each perturbation as a function of the number of e-folds.Right: Effect of varying the initial conditions; see the discussion below Eq. (3.14).N cross denotes the time at which the scale k indicated in each panel crosses the horizon (k = aH).The matrix formalism of Section 3.1 has been used for both panels.

Figure 4 :
Figure 4: Top: stochastic average of the power spectrum of the comoving curvature perturbation, P R , for 22 different values of the comoving wave number, k, (dark blue dots).The number of realizations for each value of k is 2160.Each realization is represented as light blue dot.The solid black line represents the average of the power spectrum obtained via the deterministic matrix differential equation derived in Section 3.1.Bottom: Absolute value of the relative difference between the stochastic average and the matrix average of P R .The agreement for each k is at the percent level.

Figure 5 :
Figure 5: Left: Histogram of log 10 P R (k) for 2160 realizations for k = k peak , together with a skewnormal fit for the probability distribution function.Right: Histogram of the k-independent variable ∆ log 10 P R defined in Eq. (3.18) of 20 × 2160 realizations, together with a universal skew-normal fit for the probability distribution function.

4 ×Figure 6 :
Figure 6: Left: The dissipative coefficient Γ, the quantity Σ defined in (3.22) and the function |Γ ϕ |H −1 |dϕ/dN | as functions of the number of e-folds (and in units of M p ) for the model given by Eqs.(2.21) and(2.22).The corresponding approximations as piecewise-constant functions from the table on the right are shown with dashed lines.Right: Benchmark parameters for the analytical calculation of the power spectrum.We take phases 1 and 2 to end at N 1 = 2.5 and N 2 = 5, respectively (we measure the number of e-folds from the end of phase 0, so that N 0 = 0, and we normalize a(N 0 ) = 1, see the figure on the left, as well as H = 7 × 10 −6 M p .

6 Figure 8 :
Figure 8: Left: Maximum value of the integration kernel in Eq. (4.48) along the surface y = z for the parameters in (the table of ) Fig.6.Right: Gravitational wave signals induced during and after inflation compared with the LISA sensitivity curve[31].

Figure 9 :
Figure 9: The different components entering into the quantity S GW defined in Eq. (4.49) are shown as a function of time for two different values of k.We use the numerical results of Section 3.1.

B. 4 )Figure 11 :
Figure11: Top: stochastic average of the power spectrum for 22 different values of k (dark blue dots).The number of realizations for each value of k is 900 (represented as the small, light blue dots).The solid black line represents the average of the power spectrum obtained via the deterministic matrix differential equation and the grey dashed line the results from the analytical approach.Bottom: Absolute value of the relative difference between the stochastic average and the matrix average of the power spectrum.The agreement for each k is at the percent level around the peak and of order 10 − 20% away from the peak.

2
η < c but c < η ′ < b, so that the regions in the second double integral after the first equality do not overlap, making it impossible to satisfy the constraint η = η ′ imposed by the δ function and forcing the integral to vanish.An analogous result holds for a < b < c. and the f , g and h functions are defined as in Section 3.1.The vectors B and C are (dϕ/dN ) 2 − 4ρ r −3H 2 (dϕ/dN ) 0

i g 2 i m i ω 2 i
cos[ω i (t − s)] , (F.11) which implies ⟨ξ(t)ξ(s)⟩ = T γ(t − s).(F.12) We use this result in Section 3.1 to change the time variable from cosmic time to e-folds.Itô's rule also allows us to redefine the dependent variable of an stochastic differential equation.Suppose, for instance, that we want to write Eq. (C.4) in terms of a certain function z(y, t).Let us start by considering the Taylor expansion of z(y, t) to second order in non-infinitesimal increments of the stochastic variable y and time t, ∆W 2 t and higher powers all vanish in the continuum limit.For ∆t, this is immediate (indeed, T /N → 0 as N → ∞).For ∆W3t /∆t, the fact that it vanishes follows from Itô's rule: ∆W3t /∆t = ∆W t × ∆W 2 t /∆t → 0 × 1 = 0. Invoking again Itô's rule, ∆W 2 t goes to ∆t in the continuum limit, and therefore it also vanishes.On the other hand, the ratio ∆W t /∆t can be formally associated with the time derivative of the Wiener increment in the continuum limit, i.e.This equation, valid for any z(y, t), is known as Itô's lemma, and is the stochastic equivalent of the chain rule.As we will now see, Itô's lemma has a very useful and direct application.The main quantity of interest when solving a stochastic differential equation is the probability distribution P (y, t) for the stochastic variable to take the value y at time t.This distribution can be found via the Fokker-Planck equation, which we now derive.Let us take the expectation value of both sides of Eq. (C.18) with respect to the probability distribution P (y, t), assuming that z is not explicitly time-dependent, S to denote the expectation value with respect to P (y, t).Here we have used the fact that ⟨dW t ⟩ S = 0 in the first step and integrated by parts in the last step assuming that the probability distribution vanishes at the boundaries.On the other hand, the time derivative of the mean of z is also given by 17)⟩ = δ(t − t ′ ) , (C.17)i.e. a standard white noise process.In the remaining of this appendix and everywhere else in this work we use the notation ξ t = ξ(t).Taking the above considerations into account, 4(dϕ/dN )f ρ 0 h ψ + 4(dϕ/dN )f ψ h ϕ + 4(dϕ/dN )f ϕ h dϕ + 4(dϕ/dN )f dϕ