WarmSPy: a numerical study of cosmological perturbations in warm inflation

We present WarmSPy, a numerical code in Python designed to solve for the perturbations' equations in warm inflation models and compute the corresponding scalar power spectrum at CMB horizon crossing. In models of warm inflation, a radiation bath of temperature $T$ during inflation induces a dissipation (friction) rate of strength $Q \propto T^c/\phi^m$ in the equation of motion for the inflaton field $\phi$. While for a temperature-independent dissipation rate ($c=0$) an analytic expression for the scalar power spectrum exists, in the case of a non-zero value for $c$ the set of equations can only be solved numerically. For $c>0$ ($c<0$), the coupling between the perturbations in the inflaton field and radiation induces a growing (decaying) mode in the scalar perturbations, generally parameterized by a multiplicative function $G(Q)$ which we refer to as the scalar dissipation function. Using WarmSPy, we provide an analytic fit for $G(Q)$ for the cases of $c=\{3,1,-1\}$, corresponding to three cases that have been realized in physical models. Compared to previous literature results, our fits are more robust and valid over a broader range of dissipation strengths $Q\in[10^{-7},10^{4}]$. Additionally, for the first time, we numerically assess the stability of the scalar dissipation function against various model parameters, inflationary histories as well as the effects of metric perturbations. As a whole, the results do not depend appreciably on most of the parameters in the analysis, except for the dissipation index $c$, providing evidence for the universal behaviour of the scalar dissipation function $G(Q)$.

The most compelling solution to the horizon, flatness, and monopole problems is inflation [1][2][3][4][5][6][7][8][9].By virtue of accelerated expansion, inflation ensures that when the process concludes, the Universe achieves a state of significant flatness, homogeneity, and isotropy at the largest observable scales.In addition to resolving these fundamental issues, inflation offers a mechanism for generating the density fluctuations that ultimately give rise to the observed large-scale structures within the cosmic web.The distinctive pattern in the fluctuation spectrum is expected to align with the angular power spectrum observed in the cosmic microwave background (CMB) by various experiments, including recent data from the Planck satellite and BICEP/Keck Array [10,11], as well as measurements of small-scale CMB fluctuations obtained by the Atacama Cosmology Telescope [12][13][14] and the South Pole Telescope [15,16].These fluctuations, originating from quantum-mechanical effects, are inherently adiabatic and are attributed to perturbations in a scalar field known as the inflaton field, which drives the inflationary expansion.
Conventional models of inflation involve a single scalar field slowly rolling down a nearly flat potential, inducing a quasi-de Sitter phase during which the Universe inflates superluminally.Once the field nears the minimum of its potential, the vacuum energy converts to radiation, thereby reheating the Universe.
A well-established alternative framework to conventional inflation is warm inflation (WI), in which the inflaton is thermally coupled to a bath of radiation [17,18], see refs.[19][20][21] for earlier attempts.The inflaton continually sources the production of radiation.The resulting radiation bath induces a dissipation (friction) term in the equation of motion for the inflaton field, and also alleviates the need for a separate reheating phase at the end of inflation.For reviews of WI see refs.[22][23][24].
The presence of a radiation bath and dissipation not only alters the background dynamics of the inflaton field but also its perturbations.Specifically, the addition of these thermal effects can have a significant impact on the primordial power spectrum.Fluctuations produced in WI models are primarily thermal in origin, with quantum fluctuations being subdominant in the limit of a large dissipation rate between the inflaton and radiation sectors.Quantitative models of warm inflationary perturbations must therefore take into account the effect of this dissipation.The established approach in WI (which we also follow in this work) is based on the stochastic inflation formalism, in which short-wavelength dynamics backreact on the longwavelength modes via quantum and thermal noise terms in their equations of motion [25].The perturbations' evolution is then determined by solving the fully coupled equations for all scalar, radiation energy density, and radiation momentum perturbations.
Although an accurate calculation of the production of perturbations and the resulting primordial power spectrum in WI must be done numerically, one can gain insight by parameterizing the power spectrum analytically as [25][26][27][28]: where φ is the inflaton's background velocity, n BE = [exp(H/T )−1] −1 is the Bose-Einstein distribution, T is the radiation temperature, Q ≡ Υ/(3H) is a dimensionless parameter defining the strength of the dissipation, Υ is the dissipation rate, and H is the Hubble parameter during inflation (setting the scale of inflation); all quantities are evaluated at horizon crossing.The parenthetical terms correspond to the analytic power spectrum of a temperature-independent dissipation rate, and the multiplicative factor G(Q) accounts for the growth or suppression of inflaton fluctuations due to its direct coupling with radiation for a temperature-dependent dissipation rate.Here, we refer to G(Q) as the scalar dissipation function.Determining G(Q) can only be done numerically by solving the full set of perturbation equations found in WI and is computationally demanding, particularly for large values of Q.
In this work, we present-to our knowledge-the first open-source and publicly available code to compute G(Q) and the full, numerical power spectrum for warm inflationary models. 1  We present all analytic and numerical steps for this computation in a fully reproducible manner, with a goal of providing a clear guide for readers interested in performing their own warm inflationary computations.The novelty in our work is the public code as well as the extensive study done on the relationship between the scalar dissipation function G(Q) in eq.(1.1) and other model parameters.The system of perturbation equations was previously well known; we specifically follow the approach of [25].However, their solutions had not been studied numerically in such a comprehensive way as we do in this paper.Further, our numerical analysis extends to higher values of Q that had previously been studied.
Several codes that solve for the evolution of the inflaton field are already present in the literature, in the context of conventional (cold) inflation models without dissipation.Publicly available numerical packages that yield the power spectrum for multi-field inflation models include Pyflation [29], Inflation.jl[30], and MultiModeCode [31].Various codes have been implemented to numerically evaluate inflationary correlation functions for single-field models, including CppTransport [32], PyTransport [33], see also refs.[34,35], BINGO [36], the work by Chen et al. [37], and the work by Horner and Contaldi [38].Our work differs from these previous installments because it is designed to solve the equation of motion for the inflaton field along with the radiation bath content predicted in WI.The code includes various models for the inflaton potential and the coupling of the inflaton field to the radiation bath.The output generates the fitting formulas for the scalar dissipation function G(Q).
We present our work as follows.Section 2 summarizes the background dynamics of WI and presents the notation used.In section 3, we introduce the equations for the perturbations including the stochastic inflation formalism; we provide a simple proof to obtain the analytic expression that approximates the power spectrum of WI (eq.(1.1)), and define the scalar dissipation function G(Q).The various modules that compose our code WarmSPy are described in detail in section 4. The results for the function G(Q) using different models and inflationary histories are compared in section 5, where we also propose a new fitting function that extends currently available fits in the literature to wider range of values of Q, of relevance in some WI constructions in the literature [39][40][41].

Background evolution in warm inflation
In WI, the inflaton continually converts into radiation during inflation itself, as opposed to only reheating the universe once inflation ends.This dissipation mechanism is parameterized by the introduction of a non-negligible dissipation rate Υ = Υ(T, φ) that couples the inflaton field to the radiation energy density ρ r in their equations of motion [17,42]: where a dot denotes the derivative with respect to cosmic time and V ,φ = ∂V /∂φ.Conservation of the total energy of the system imposes that the energy lost by the inflaton field must be gained by the radiation's energy density, with the term on the RHS of eq.(2.2) encapsulating this process.The dissipation rate Υ describes the rate at which the energy in the inflaton field converts into radiation, and the detailed form of this term is model-dependent [43,44].Here, we parametrize the dissipation rate as: where C Υ is a dimensionless constant that depends on the underlying microphysical model, M is a mass scale specified by the model, while c and m are integers that can be either positive or negative.Successful models of WI have been built by coupling the inflaton with the radiation bath indirectly through heavy mediators [22,45,46] and directly by protecting the flatness of the inflaton potential from thermal and quadratic corrections [40,41,47].A discussion on the derivation of the dissipation coefficient using field theory tools is given in [48].
In the following, we will work with c = {3, 1, 0, −1}, motivated by the explicit constructions of dissipative rates in the literature.For example, a decay process of the inflaton into light fermions mediated by a heavy boson field predicts Υ ∝ T3 [46].Further, c = 3 also results when considering an axion field coupled to a non-Abelian gauge group, with the friction term related to the sphaleron transition rate obtained from lattice gauge theory [49,50].In the presence of a small mass to the final fermionic states, chirality-violating scattering suppresses the friction from the sphaleron transitions leading to Υ ∝ T [51].Moreover, in the warm little inflaton model, where the inflaton field is still taken to be a pseudo Nambu-Goldstone boson, Yukawa coupling to a few light fields have been shown to produce a dissipation coefficient with linear temperature dependence, i.e.Υ ∝ T [47], and in the high temperature regime an inverse temperature dependence, i.e.Υ ∝ T −1 [41].
The strength of the dissipation is typically parametrized by the dimensionless ratio Q: For Q 1, a strongly dissipative regime is achieved, while Q < 1 represents the weak regime of WI.In all cases, WI requires T > H, which is roughly the criterion for which thermal fluctuations dominate over quantum fluctuations [17].Throughout this work, we also assume that the radiation thermalizes on a time scale much shorter than 1/Υ [17,42], so that the energy density of radiation can be taken to be: where g * (T ) is the number of relativistic degrees of freedom for a radiation bath at temperature T . 2he Hubble expansion rate H ≡ ȧ/a is determined by the Friedmann equation, which reads: where GeV is the reduced Planck mass.Inflation takes place when the potential V is approximately flat and the potential energy dominates, such that the Hubble expansion rate H is approximately constant.During this period, known as the slow-roll regime of the inflaton field, higher order derivatives in eqs.( This regime can be parametrized by a set of slow-roll parameters , η and β, defined by: 3 which in order to meet the conditions in eq.(2.7), must themselves satisfy: [53][54][55].In other words, inflation quickly comes to an end when one of the slow-roll parameters defined above violates the slow-roll condition, i.e. = 1+Q, or |η| = 1+Q or |β| = 1+Q.More rigorously, the acceleration expansion continues as long as: and terminates when H = 1.
To summarize, the background dynamics of WI obeys the system of three ordinary differential equations described by eqs.(2.1), (2.2) and (2.6)-together with two constraints in eqs.(2.3) and (2.5)-which are evolved until H = 1.We can rewrite this system of equations in terms of the number of e-folds: which quantifies the amount of expansion that the scale factor undergoes during the accelerated expansion.Denoting derivatives with respect to N e with the dagger symbol (not to be confused with Hermitian conjugation), the system of differential equations for the background evolution reads: ) with φ = Hφ † , ρr = Hρ † r , and φ H 2 φ † † where we neglect slow-roll corrections. 4 Cosmological perturbations in warm inflation

Basic framework
The presence of a radiation bath and a dissipation rate not only alters the background dynamics of the inflaton field but also its perturbations.Specifically, thermal fluctuations of the radiation fluid source density perturbations that are subsequently transferred to the inflaton field as adiabatic curvature perturbations.This contrasts with standard cold inflation, in which the seeds for large-scale structure formation are simply generated by quantum fluctuations.To investigate scalar perturbations in WI, throughout this paper we follow the approach of [25].In this section we review the system of perturbation equations and solutions from that paper, although our notation differs in eqs.(3.16) to (3.23) in that we choose an explicit parameterization with respect to τ = ln(k/aH) for convenience.Further, when we solve the system of equations, to simplify the numerical implementation in sections 4 and 5, we set the dissipation strength Q to a constant value during the inflationary expansion, instead of letting it evolve as in other work [25]. 5In short, the system of perturbations' equations we investigate here is already well established.However, a comprehensive numerical examination of their solutions had not been performed to date.One of the scopes of our work is to fill this gap.Specifically, the key contribution of our paper, in addition to the public code presented in section 4, is the extensive study on the relationship between the scalar dissipation function G(Q) in eq.(1.1) and other model parameters.Further, our numerical analysis extends to higher values of Q that what has been explored to date.All the aforementioned contributions are thoroughly presented in section 5.
Returning now to the basic setup, we consider the perturbed FLRW metric in Newtonian gauge: We expand the inflaton field, radiation energy density, and radiation pressure around their background values: p r (x, t) = pr (t) + δp r (x, t) . (3.4) The evolution equations for the field and radiation perturbations follow from conservation of the energy-momentum tensor.In momentum space, the equations of motion for the radiation perturbations with comoving wavenumber k read [28]: where Ψ r is the momentum perturbation, w r 1/37 is the equation of state for the radiation fluid, and J r = −Υ φδφ and Q r = Υ φ2 are the momentum and radiation source terms, respectively.The radiation source perturbation δQ r is given by [28]: where according to its parametrization defined in eq.(2.3).
In addition to eqs.(3.5) and (3.6), we also have the evolution equation for the inflaton field perturbations δφ.This is obtained by perturbing the inflaton field's equation of motion about its background value, and adding two stochastic Gaussian noise terms corresponding to the quantum and thermal fluctuations [25]: ξ T describes the thermal noise fluctuations in the local approximation for the effective equation of motion of the inflaton field.It is connected to the dissipation coefficient through a Markovian fluctuation-dissipation relation, which in an expanding space-time is given by [25]: On the other hand, ξ q describes the back-reaction of short-distance quantum modes on longdistance field modes as in the standard stochastic inflation formalism.ξ q can thus be interpreted as a Gaussian white noise term, with two-point correlation function [24]: where n * denotes the inflaton statistical distribution due to the presence of the radiation bath.This is generally assumed to be the equilibrium Bose-Einstein distribution, i.e. n * = [exp(H/T ) − 1] −1 .We emphasize that eq.(3.11) is an approximation of the quantum noise term's two-point correlation function.The quantum noise amplitude is chosen to reproduce the analytical estimate for the scalar power spectrum in WI derived in [25], under a rigorous treatment of the stochastic inflation approach to WI perturbations.We provide an explicit proof of this in section 3.2.
The complete system of first order perturbation equations for WI is obtained by combining eqs.(3.5)-(3.9).To solve this system of equations, we perform a convenient change of variables from cosmic time t to: τ ≡ ln z, where z ≡ k/(aH) . (3.12) It can be easily shown that for any function f (t) it follows: where ≡ d/dτ , H has been defined in eq.(2.9), and η H ≡ H −1 d ln H /dt. Since the perturbations' equations are evolved until the CMB scale k * exits the horizon, approximately 60 e-folds before the end of inflation, both H and η H are 1 and can therefore be neglected.8Combining eqs.2.10 and 3.12, one notes that τ is related to the number of e-folds by a simple linear transformation: where τ k represents an arbitrary additive constant that can be chosen to ensure that τ denotes the number of e-folds as measured before horizon crossing of the k−mode.
The scalar power spectrum is determined from the definition of the comoving curvature perturbation R, rescaled by k 3/2 : where . . .stands for the ensemble average over different realizations of the noise terms and:

R =
i=φ,r ) with φ2 + V and pr ≡ p r = 1 3 ρ r .Finally, it is noteworthy to emphasize that by opting for the time variable τ and rescaling all the dynamical variables by k 3/2 , we effectively eliminate the k-dependence both in the perturbations' equations (eqs.3.20 to 3.23) and in the equation for the resulting scalar power spectrum (eq.(3.24)).Consequently, this approach enables the evaluation of a universal, kindependent solution to the equations, while also providing a convenient means to reintroduce the k-dependence by setting τ = 0, i.e. k = aH, at the time a specific mode k of interest crosses the horizon.In our case for the CMB mode k * , τ (k * ) = 0 corresponds to roughly 60 e-folds before the end of inflation.
Recently, ref. [56] presented an alternative formulation of the perturbations' equations (eqs.3.2 to 3.7 in their manuscript) that makes use of a Fokker-Planck approach to obtain a system of deterministic linear differential equations for the two-point correlation function of the perturbations.This approach-previously introduced in ref. [57]-differs from the methods employed in this work and consists of an interesting and novel pathway to solve the system of stochastic differential equations for the perturbations numerically.By employing this method, the need to repeatedly solve the complete system of stochastic differential equations for the perturbations is effectively circumvented.
In practice, a fair method for the assessment of the different approaches would consist e.g. in comparing the computational resources to achieve the desired accuracy.However, a discussion along this line for the Fokker-Planck approach has not been provided in ref. [56].While we have detailed the specifics of our parallelized code in section 4.3, we leave the study of the effectiveness of this formalism and its potential implementation in our publicly available code to future work.
In addition to the method employed, the formulation in ref. [56] departs from the one presented here in two ways: (1) the quantum noise term ξ q in the inflaton field perturbations is neglected, and (2) the thermal noise term ξ T has been incorporated in the radiation perturbations.Apart from these discrepancies, we have checked that the results from the two approaches agree as soon as the strong dissipative regime Q > 1 is considered, which is the relevant regime for the purposes of this work.Nevertheless, we briefly comment on these aspects.
With respect to point (1), in our framework, the quantum noise term is necessary to accurately consider the impact of sub-horizon quantum fluctuations on the behavior of long wavelength modes in the inflaton field.By doing so, we achieve a comprehensive understanding of the combined effects arising from both inflaton quantum and thermal fluctuations, succinctly captured by a single stochastic Langevin-like equation of motion.Notably, in this formulation, the cold inflation limit (Q = 0) emerges naturally, as eq.(3.20) aligns with the conventional Starobinsky stochastic inflation approach [58] in this regime.In contrast, in ref. [56] the cold limit is attained from the homogeneous solution of the Langevin equation for the inflaton perturbations sourced by thermal noise alone.This is achieved under a specific choice of inflaton perturbations quantization, which allows for the separation of contributions to the power spectrum into an homogeneous and inhomogeneous component.The two approaches exhibit inconsistency within the weak dissipative regime of warm inflation for Q 1. Specifically, in ref. [56], the scalar power spectrum generally appears smaller compared to our framework, seemingly underestimating the non-negligible contribution from the thermal fluctuations in the radiation bath and the thermal excitation of the inflaton perturbations for T > H.
Regarding point (2), the addition of a thermal noise term in the radiation perturbations had already been thoroughly studied in ref. [59] proving that the thermal noise has only a minor influence on the resulting scalar power spectrum, which becomes negligible for a nonzero temperature power of the dissipation rate, i.e. c = 0.9 Bastero-Gil et al. also emphasize in ref. [59] the ambiguity of incorporating this term, which comes from the fact that Einstein equations do not fully fix whether the dissipative noise source enters in the energy or in the momentum flux.In conclusion, the decision to include or neglect the thermal noise term in the radiation perturbations does not undermine either approach, rather presents a relatively minor yet unresolved matter within the context of warm inflation that does not modify the outcome of this work.

Analytic estimate of the power spectrum
To obtain the scalar power spectrum, one needs to solve the system of differential equations in eqs.(3.20)- (3.23).In general this can only be done numerically, since the radiation, inflaton and metric perturbations' equations are all mutually coupled.Specifically, from eqs. (3.20) and (3.21), we see that the inflaton fluctuations couple directly to the radiation fluctuations when c = 0.As first shown in [27], this coupling results in a growing (decaying) mode in the curvature power spectrum for c > 0 (< 0).
In this subsection, we will restrict our discussion to the one specific case in which an explicit analytic expression for the scalar power spectrum can in fact be obtained: the case where the temperature power of the dissipation rate is set to zero (c = 0) and the metric perturbations are set to zero by imposing ϕ = 0 in the Newtonian gauge. 10This second approximation is certainly valid for small to moderate values of the dissipation strength Q, as one can easily show that the coefficients of the ϕ and ϕ terms in the perturbations' equations are first order in the slow roll parameters and η.With these simplifications, the equations for the inflaton and the radiation perturbations decouple and one can obtain the curvature power spectrum using Green's function techniques.This was first done in [25], in which the authors also presented a rigorous treatment of the quantum noise term using the stochastic inflation formalism.Below, we provide a simplified proof that reproduces the analytical estimate obtained in [25], using the approximate Gaussian two-point correlation function for the quantum noise from eq. (3.11), thus justifying its form.
Following the notation in ref. [25], the equation satisfied by the inflaton perturbations δφ in eq.(3.20) can be written in terms of the variable z ≡ k/(aH) as: which has general solution: Here, the Green's function Φ(z, z ) is given in terms of Bessel functions: which is solved within z > z and where: Using Eq. (3.28) and the definition of the inflaton power spectrum ∆ 2 δφ in terms of the two-point correlation function for the scalar field perturbations, we obtain: In the last line, we use the expressions for the two-point correlation functions of the thermal and quantum noises in eqs.(3.10) and (3.11), and we assume that Υ, T , H and Q are approximately constant with respect to z.Since z 1 for the relevant observational scales for the inflaton perturbations, we can approximate the integral appearing in eq.(3.33) as: Here, Γ(x) is the Gamma function.To obtain the final equality, we (1) use the asymptotic form for the Bessel function Y α (z) when z 1, and (2) neglect the terms proportional to the slow-roll coefficients, such that In a spatially flat gauge, the comoving curvature perturbation takes the simple form R = Hδφ/ φ [60].From this it follows: In summary, the above analytical estimate for the scalar power spectrum in WI is only valid if the temperature power of the dissipation rate is to zero, i.e. c = 0, and if the metric perturbations are neglected, i.e. setting ϕ = 0 in the Newtonian gauge.In addition to these primary assumptions, to obtain the analytical expression in Eq. (3.37), we had to additionally assume that: (1) the temperature T and Hubble parameter H are constant with respect to z in order to bring them outside of the integral in eq.(3.33); (2) the term proportional to the slow-roll coefficient η could be neglected as such that we could set α ≈ ν.As we will show numerically in section 5.4, these two secondary assumptions, as well as neglecting the metric perturbations, do not significantly affect the resulting scalar power spectrum.In contrast, deviating from a temperature power of the dissipation rate equal to zero, i.e. c = 0, notably impacts the scalar power spectrum, particularly in the regime characterized by strong dissipation.This observation is what originally motivated the introduction of a multiplicative function of the dissipation strength G(Q), in conjunction with the existing analytical expression.
For general cases with c = 0, the analytic expression in eq.(3.37) must therefore be multiplied by a factor G(Q) which can only be determined numerically.The subsequent section delves into the definition and observational relevance of G(Q), which we also refer to as the scalar dissipation function.

The scalar dissipation function G(Q)
The analytic form of the power spectrum given by eq.(3.37) is the result used in most of the recent literature in WI.As emphasized above, once we include the coupling between the inflaton and radiation perturbations, i.e. set c = 0, and restore the metric perturbations, these equations can only be solved numerically.This results in a correction to eq. (3.37) which is usually parametrized by a multiplicative function of the dissipation strength G(Q), determined by fitting the numerically computed curvature power spectrum for a given form of the dissipation rate.In formulas, this reads: Once the scalar dissipation function G(Q) is specified, we can compute the spectral tilt n s and the tensor-to-scalar ratio r, and compare the obtained results with the observational constraints from the CMB at the pivot scale k * : Here, ∆ 2 T is the tensor power spectrum which in WI remains effectively unaltered relative to the standard cold inflationary scenario, since the gravitational interactions are weak. 11hus, the curvature power spectrum is the quantity primarily affected by the dissipation in WI, which in turn significantly alters the resulting spectral tilt n s and tensor-to-scalar ratio r, relative to standard cold inflation.In order to test WI models against current data, it is therefore crucial to obtain precise numerical fits of G(Q) for different forms of the dissipation rate and check their robustness across different inflaton potentials, model parameters and inflationary histories.This constitutes the main subject of this paper and of the publicly available code that we describe in detail in the following section.

Code implementation
The bulk of the code is contained in the Python script WI_Solver_utils.py,which is structured in four modules: (1) Inflaton_Model, (2) Background, (3) Perturbations and (4) Scalar_Dissipation_function, which respectively (1) define the main parameters for a given inflaton potential, solve for its (2) background and (3) perturbations' evolution, and finally use this information to (4) fit for the scalar dissipation function G(Q).An important simplification we make throughout this code is to take Q as a constant during the inflationary expansion.In section 5.1 we show explicitly that this assumption has no impact on the obtained fits for G(Q), while greatly simplifying the numerical implementation.e

Inflaton_Model
In the Inflaton_Model module, we define the inflaton potential and consequently all of the parameters that depend on it (such as the Hubble parameter H, its derivative with respect to cosmic time Ḣ, etc.), that will be used for all the subsequent computations for the background and perturbations' evolution.Several potentials have been proposed to embed the warm inflationary dynamics into a larger framework.Here, we review the ones currently available in the code, which generally represent some of the most studied cases in the literature.Monomial potential.The simplest type of potential we consider is the power-law form: where the index n describes the steepness of the potential and λ is dimensionful except when n = 4, for which it describes the self-coupling of the inflaton field.The form in eq.(4.1) has been adopted as the inflationary potential since the early conceptions of the theory, although it is now strongly disfavored by CMB data in combination with BAO and BICEP2/Keck Array data in the context of the standard cold inflation [11].In contrast, in WI, due to the different dynamics resulted from the dissipative effects, monomial potentials, specifically with even power law n, are consistent with CMB observations for a large range of the parameter values [66].
Hilltop-like potential.One of the simplest alternatives to monomial potentials is of the form: with n ≥ 1 and with the inflaton rolling off the top (plateau) of a potential, i.e. starting at φ φ f [67].We assume φ f to be sufficiently large such that inflation ends before the inflection point of the potential is reached.For a potential of the form in eq.(4.2), the field excursion during inflation is much smaller than the Planck scale, which is generally preferred from an effective field theory perspective.Hilltop-like potentials were found to be consistent with CMB data in both cold and warm inflation for wide regions of parameter space [66,68].
Natural inflation.The model of Natural Inflation [69] uses axion-like particles as natural candidates for the inflaton, in that their potentials possess an underlying shift symmetry that protects the required flatness of the potential to achieve both sufficient inflation and the correct amplitude of density perturbations for structure formation [70,71].The original variant of Natural Inflation used a cosine potential motivated by the case of the QCD axion: where Λ is related to the magnitude of the explicit breaking scale and N is the number of vacua.The model avoids the η-problem intrinsic to cold inflation and introduces a nearly flat potential suitable for slow-roll inflation.These symmetry properties also suppress radiative and thermal corrections which allow this model and more generally axion-like potentials to be trivially implemented in the context of WI [40,[72][73][74].This potential in the WI scenario was found to be consistent with observations for super-Planckian and marginal sub-Planckian values for the decay constant f [52,75].
β-exponential potential.A generalization of power-law (runaway) inflation models is the β-exponential potential of the type [76]: which arises from theories of gravity with higher derivative terms or compactified Kaluza-Klein models.Here, λ is a dimensionless coupling and the constant β parametrizes the deviation from the pure exponential model (for β → 0).In cold inflation scenarios, CMB data favors a non-minimal gravitational coupling of the inflaton in the case of β-exponential potentials [77], while these models are consistent with WI with a cubic dissipation rate [78].

Background
The Background module solves for the background evolution as a function of the number of e-folds N e according to the system of equations in eqs.(2.11)-(2.13),for a given choice of the inflaton potential and its parameters.The solver evolves the system of equations from N e = 0 to N e = N end e ≡ N pert e + N infl e , with the initial conditions chosen such that H = 1, i.e. inflation ends, when N e = N end e .N infl e refers to the number of e-folds between horizon crossing of the CMB mode k * and the end of inflation.We generally take this to be 60 to match CMB observations.N pert e refers to how long before horizon crossing is an individual mode initialized when evolving the perturbations.This sets the duration of numerical evolution of a particular mode-we take this value to be approximately 5 to 10 e-folds.Note that this choice does not necessarily initialize the perturbations at the start of the background solution, as doing so can lead to numerical stability issues.We elaborate on this in Section 4.3 below. 12 To find the value of φ 0 that guarantees N end e e-folds before inflation ends, we define the function Bg_solver_test which is an iterative solver of the background evolution that, starting from an initial guess φ guess 0 , scans the parameter space in φ via a standard gradient descent algorithm until the condition H (N end e ) = 1 is satisfied; we refer to Appendix B for more details.
For given values of Q and φ 0 , the additional required initial conditions on ρ r,0 and φ † 0 are determined from eqs. (2.11)-(2.13)assuming slow-roll, i.e. for φ † † 0 = 0 and ρ † 0 = 0.In formulas: Once the initial conditions are determined, we use the function Bg_solver_tau which additionally outputs all the relevant background quantities for the perturbations' equations as a function of τ ≡ ln(k/(aH)), such that they can be input directly in the perturbations solver. 13  12 To be more precise, the perturbation equations are only evolved for a specified interval in τ ∈ [τmax, τmin] that contains τ = 0, i.e. the time corresponding to CMB horizon crossing.The number of e-folds N pert e is then safely chosen such that it contains said τ range up to CMB horizon crossing, i.e.N pert e ⊇ [τmax, 0].See section 4.3 for more details. 13Note that solving the background with respect to the number of e-folds Ne is a convenient choice since as shown in eq.(3.15), Ne is linearly related to the time variable τ .This implies that the background evolution with respect to Ne and τ is the same up to a linear transformation of the independent variable, i.e. the background solution can be directly inputted in the perturbations' equations.
We emphasize here that the simplifying assumption of a constant dissipation strength Q during the inflationary expansion implies that the background evolution is independent of both the temperature and inflaton field's power-law exponents in the dissipation rate Υ, respectively c and −m according to eq. (2.3).This feature of the code is crucial as it allows for an efficient evaluation of the scalar dissipation function G(Q), as we will show in greater detail in sec. 5.However, it also implies that currently WarmSPy should not be directly employed to study the background dynamics in WI, where the explicit evolution of Q is essential.Note, that the effectiveness of our code resides in the modularity and the capability to solve for the inflaton perturbations given various backgrounds for the evolution of the inflaton field.In any case, in future updates of WarmSPy, we plan to enhance the code's capabilities by incorporating the explicit evolution of Q and enable a precise study of the background dynamics in WI.

Perturbations
The Perturbations module takes as input (1) the background evolution computed in the Background module for a given choice of the inflaton potential, and (2) the values for c and m, respectively the temperature and negative inflaton field's power-law exponents in the dissipation rate Υ.The module numerically solves the perturbations' stochastic equations of motion, (3.20)-(3.23),for a specified interval in τ ∈ [τ max , τ min ] divided into N τ steps.This interval is chosen to include τ = 0, the time corresponding to horizon crossing of the mode being evolved.Since in this study we are interested in the CMB mode k * , we specifically set τ = 0 at the time corresponding to N inf e e-folds before the end of inflation, i.e. when the CMB mode k * crosses the horizon, in formulas this corresponds to: 14 τ (N pert e ) = 0. (4.6) The system of equations is then solved over N ξ realizations of the stochastic noises and from that the scalar power spectrum at τ = 0 is computed according to its definition in eq.(3.24), i.e. as the ensemble average of the comoving curvature perturbation R. The default τ interval is given by {τ max , τ min } = {6.5, −1}.Additionally, we choose N τ = 5 × 10 5 time steps and N ξ = 1024 realizations of the noise terms.The value of τ min is set so that all modes are evolved for 1 e-fold after horizon crossing.On the other hand, the choices on τ max , N τ and N ξ were made in order to ensure numerical stability in the solution, without compromising the overall speed of the code.
In this regard, the selection of τ max plays a crucial role in preserving the analytic accuracy and numerical stability of our solution, as it is closely tied to the initialization of the k-mode evolution.As noted in Section 4.2, we do not necessarily initialize the mode at the start of inflation.Rather, we start the mode's evolution at some fixed time between 5 and 10 efolds before horizon crossing.This choice becomes particularly important for scenarios where inflation lasts a considerable duration, such as O(1000) e-folds.In such cases, the term e 2τ in the equation of motion (3.20) becomes extremely large at early times.Recalling that τ denotes the number of e-folds as measured before horizon crossing (see eq. (3.15)), this yields an O(e 2000 ) term in eq.(3.20) when the mode is initialized.The level of precision required to resolve both these early-time large terms and the other smaller terms that are relevant at late times would dramatically slow down any numerical evolution of the perturbations' equations.To overcome this issue, we exploit the fact that most of the subhorizon evolution is trivial, since the mode will merely behave as a Minkowski harmonic oscillator in the deep subhorizon [79].If we initialize the mode sufficiently early that it still falls in the Minkowski oscillator limit, our choice of initial conditions will still be accurate and applicable.This typically holds around 5 or more e-folds before horizon crossing, which eliminates the need to begin the mode's numerical evolution at the start of inflation.In the above example of inflation lasting O(1000) e-folds, instead of the aforementioned O(e 2000 ) blow-up in the equations of motion, this simplification only has a term of O(e 10 ), which is much more manageable numerically.With our default choice for τ max , the evolution of the perturbations' equations is initialized 6.5 e-folds before horizon exit, which recovers the subhorizon oscillator limit just discussed.In section 5.3, we investigate the robustness of the obtained G(Q) fits for different choices of τ max in greater detail.
As already discussed in section 3.1, to solve for the evolution of the perturbations we also need to specify a set of initial conditions, which for simplicity are taken to be: where T,0 is determined via eq.(3.19), for a given value of the dissipation strength Q and the corresponding background quantities H(τ max ) and T (τ max ) computed in the Background module.We emphasize here that the choice of initial conditions is actually not relevant.In fact, the stochastic sources dominate the evolution of the perturbations, washing out any information regarding the initial state of the perturbations [28,57].
This module also provides the option to turn on or off (1) the metric perturbations, i.e. setting ϕ = 0 or = 0, and (2) the corrections of first and second order in H and η H from the chain rule when changing variables from cosmic time t to τ ; see eqs.(3.13) and (3.14).By default, we neglect both effects, since, as we will show explicitly in section 5.3, including them merely slows down the execution of the code without impacting the obtained fits for G(Q).
Finally, the evaluation of the stochastic differential equations is completely parallelized via the function Pool_Solver, over the number of available cores N cores in a given machine.This significantly reduces the total execution time of the code roughly by a factor of N cores .For a machine with N cores = 128 and given the default τ interval with N τ = 5 × 10 5 , it takes roughly 20 seconds to solve the perturbations' equations N ξ = 1024 times.
To summarize, the Perturbations module outputs the scalar power spectrum at CMB horizon crossing and its standard deviation as a function of the dissipation strength Q, which in this work we take to be in the interval Q ∈ [10 −7 , 10 4 ].

Scalar_Dissipation_function
The Scalar_Dissipation_function module takes the background and perturbations' evolution as inputs, computed respectively in the Background and Perturbations modules for a given choice of the inflaton potential and values of c and m.The scalar dissipation function G(Q) is then computed according to its definition in eq.(3.38).This module then determines the best fitting function G(Q) via the method of least squares.Specifically, to allow for a direct comparison with previous work, we employ the polynomial fitting functions most prominently used in the literature due to its simplicity [47,52,66,80]: In the case of an inverse temperature dependence in the dissipation rate (c = −1), G(Q) pol fits our data well.On the other hand, for c > 0, this fitting function does not robustly capture the behaviour of the scalar dissipation function over the entire range of values of Q investigated, i.e. for Q ∈ [10 −7 , 10 4 ].To correct for this, when c > 0 we propose a logarithmic fitting function of the form: As we show in detail in section 5, G log (Q) fits the data well over the entire range of values of Q investigated, for the same number of free parameters as the polynomial fit, totaling four parameters.

Results
In this section, we present the numerical results, obtained with WarmSPy, for the scalar dissipation function G(Q).We provide precise fits of G(Q) for different values of the power-law temperature dependence in the dissipation rate, namely c ∈ {3, 1, −1}, and investigate the relationship of the dissipation function G(Q) to additional model parameters, such as the inflaton potential and the inflaton field power-law dependence of the dissipation rate.Unless otherwise stated, all results presented below were obtained with all the main parameters of the code set to their default values, specifically {τ max , τ min } = {6.5, −1}, N τ = 5 × 10 5 , N ξ = 1024, N pert e = 10 and ϕ = 0.15

Method validation and comparison with previous literature
As mentioned in section 4, the evolution of the dissipation strength Q is neglected during the inflationary expansion.Such an approximation is one of the main differences in the methodology between our work and the previous literature.Our choice greatly simplifies the evaluation of the initial conditions to obtain N tot e e-folds of inflation while keeping a high level of accuracy in the numerical solutions as discussed in the following sections.In fact, if we let Q vary according to eq. ( 2.3), the background evolution becomes dependent on the exponents c and m.This increases the total number of computations required to obtain the scalar power spectrum for different combinations of c and m, as one needs to redetermine the initial conditions each time.On the other hand, with our method, for a given inflaton potential and value of N tot e , we only need to find the initial conditions once, regardless of the values of c and m.
Physically, taking Q to be constant is equivalent to enforcing a specific evolution of the dimensionless constant C Υ in the dissipation rate Υ, see eq. ( 2.3), in such a way as to keep the same dissipation strength Q throughout the inflationary expansion.However, we would expect that once the numerical scalar spectrum is divided by the analytical result, any potential alteration to the shape of the scalar power spectrum would get washed out, effectively leading to the same G(Q).We confirm this assumption numerically by comparing the solution found with a constant Q against what is obtained by allowing Q to vary according to eq. (2.3).
In figure 1, we show the comparison between the scalar dissipation function G(Q) obtained with our code (M.A.V.F. in the legend) both for a constant and varying dissipation strength Q, and for different combinations of c and m (see the caption of figure 1 for more  [47,66,80] (labeled "Fit 1" in the legend) and in ref. [39] (labelled "Fit 2" in the legend), see text for additional details.In summary, this figure demonstrates the consistency of our code in the resulting G(Q) when either keeping Q constant or allowing it to evolve, as well as its agreement with the well-established G(Q) fits found in the existing literature.details).Results are obtained assuming that the inflaton rolls down the quartic potential in eq. ( 4.1) with n = 4 and the self-interacting constant λ = 10 −14 , and setting N infl e = 60.In addition, figure 1 shows the fitting functions G(Q) obtained in refs.[47,66,80] (labelled "Fit 1" in the legend) where the cases c = {3, 1, −1} are explored respectively, along with the analytic form given in eq.(4.8).These fitting functions were also obtained for a quartic inflaton potential in refs.[47,66,80], but for a smaller range of values Q 100.As expected, the data produced with our code match these fitting functions extremely well for Q < 100 and only diverge significantly for larger values of the dissipation strength and for c = 3.Finally, we also plot a more recent fitting function computed in ref. [39] (labelled "Fit 2" in the legend) for the case of a runaway potential, a cubic temperature dependence c = 3 in the dissipation rate, and a wider range of values Q 10 3 .This fitting function reproduces our data more accurately than the basic polynomial form, and does so up to Q ∼ 10 3 .This was not necessarily expected given that our data and the fitting function from ref. [39] were obtained for different inflaton potentials, and it hints at the universality of the scalar dissipation function G(Q) for a given value of c, an aspect we explore more in depth in the following section.Finally, it is also important to emphasize that the better fit achieved by the scalar dissipation function G(Q) proposed in ref. [39] is also due to its cumbersome form with 11 free parameters compared to the 4 from the simpler polynomial fit in eq.(4.8).
Overall, the main results from the comparison shown in figure 1 can be summarized as follows: (1) there is no difference in the resulting scalar dissipation function G(Q) when Q is kept constant or left evolving during the inflationary expansion; (2) the scalar dissipation function G(Q) evaluated with our code reproduces accurately the well-known G(Q) fits previously obtained in the literature.
This discussion also raises two pertinent issues for further investigation: (1) the feasibility of generating precise fits of the scalar dissipation function G(Q) for a large range of Q values, i.e. for Q ∈ [10 −7 , 10 4 ], with a minimal number of free parameters; (2) the relationship (if any) of the scalar dissipation function G(Q) for a set value of c to additional model parameters, such as the inflaton potential or the value of m.These aspects will be thoroughly examined in the next section.

Numerical fits of the scalar dissipation function G(Q)
In figure 2, we plot the numerical solutions obtained with WarmSPy, along with the corresponding fits of the scalar dissipation function G(Q) for different values of the power-law temperature dependence in the dissipation rate (see eq. (2.3)), namely c ∈ {3, 1, −1}.The parameters obtained from the fitting functions in eqs.(4.8) and (4.9) to the numerical results for G(Q) for shown in figure 2.
These results were obtained under the following assumptions: (1) there is no inflaton field dependence in the dissipation rate, i.e. m = 0; (2) the inflaton rolls down the monomial potential in eq.(4.1) with n = 4 and the self-interacting constant λ = 10 −14 ; (3) the number of e-folds after CMB horizon crossing is set to N infl e = 60; (4) the metric perturbations are neglected, i.e. ϕ = 0.While these assumptions might appear excessively specific and restrictive, the detailed analysis presented in section 5.3 will demonstrate that the fits of the scalar dissipation function G(Q) shown here possess broader validity, extending beyond the fulfillment of all aforementioned assumptions.In short, we find that the behaviour of the scalar dissipation function G(Q) is susceptible only to the temperature dependence of the dissipation rate c.For the case of c = {3, 1} (c = −1), this behaviour is accurately described by the fitting function from eq. (4.9) (eq.(4.8)) respectively for the parameters shown in table 1, throughout the entire range of values of Q considered.Specifically for c = {3, 1}, the fractional residuals plotted in the lower panel of figure 2 are always O(1), meanwhile for c = −1 the fit is even more robust with fractional residuals O(0.2).
For a positive temperature dependence in the dissipation rate, we also test the polynomial fitting function from eq. (4.8) which is generally unable to capture the behavior of the scalar dissipation function G(Q) over the entire range of Q values considered.While this is not evident from the fractional residuals shown in the bottom panel of figure 2 which are roughly the same order as those for the logarithmic fit; one can note that for both c = {3, 1}: (1) the best-fit for the parameter B is driven to an extreme small value ∼ {10 −14 , 10 −17 } respectively; (2) the fit largely underestimates the numerical solution in the range Q ∈ (1, 10 2 ).This seems to be at least partially related to the small dip in the magnitude of G(Q) for Q ∼ O(0.1 − 1), which cannot be captured by the simple polynomial fitting function.In addition, we also observed that the best fit parameters and the accuracy of the polynomial fit depends strongly on the range of values of Q to which it is applied.This explains the significant difference in the best-fit parameters found here against those obtained in previous work, despite the consistency of the numerical solutions in the overlapping range of values of Q, see figure 1.Given these considerations and acknowledging the value of a simple fitting function for G(Q) in calculating cosmological observables, we suggest restricting the evaluation of the polynomial fit to a narrow range in Q that is relevant to the specific inflaton model under examination.

Universality of the scalar dissipation function G(Q)
In this section, we assess the relationship of the scalar dissipation function G(Q) for a set value of the temperature power in the dissipation rate c ∈ {3, 1, −1}, to additional model parameters.Specifically, in the order in which they will be presented, we investigate the dependence of G(Q) from: (1) the negative inflaton field power of the dissipation rate m; (2) the inflaton potential; (3) the number of e-folds before and after CMB horizon crossing, i.e. by respectively changing τ max and N infl e .Finally we also analyze (4) the impact of turning on the metric perturbations, i.e. setting ϕ = 0, and/or the corrections of first and second order in H and η H (see differential chain rule in eqs.(3.13) and (3.14)).
The negative inflaton field power of the dissipation rate m: We discuss the impact of considering different values for the index m (see eq. ( 2.3) for its definition) for a set temperature dependence in the dissipation rate.Specifically, for this study we assume that the inflaton rolls down the monomial potential in eq.(4.1) with n = 4 and the self-interacting constant λ = 10 −14 , and setting N infl e = 60.In the upper panel of figure 3 we then plot the resulting scalar dissipation function G(Q) for a cubic (c = 3) and inverse (c = −1) temperature dependence of the dissipation rate and different values of m, respectively m = {2, 0} and m = {−2, 0}. 16For a given value of c, the scalar dissipation function G(Q) is independent of m.Specifically, the lower panel of figure 3 shows the residuals between the two values of m assessed for each set value of c = {3, −1}, which all lie below O(0.25) for the entire range of dissipation strength Q considered.Additionally, one can note that, while still small, the residuals for the case with m = 0 = c − 1 against the m = 0 case are mostly negative as Q 10 2 .This implies that the case with m = 0 slightly underestimates the growth (suppression) of the scalar power spectrum in the strong dissipative regime when c = 3 (c = −1).
These observations can be understood analytically by the fact that the m dependence in the inflaton field perturbation equation is slow-roll suppressed in eq.(3.20), i.e. m φQ/(φH) βQ/(1 + Q) which is negligible for Q 1.In the strong dissipative regime (Q 1), this term cannot be ignored and contributes to the growth (suppression) of the inflaton field perturbations for c = {3, −1} respectively.Interestingly, even in this regime, this term has only a minor influence in the resulting scalar power spectrum as the coefficient proportional to c multiplying δ ρr dominates the overall evolution of the inflaton perturbations.
The inflaton potential: We discuss the impact of considering different models for the inflaton potential, following the framework discussed in section 4.1.For this study, we assume m = 0 in eq.For each panel, we show the numerical results of G(Q) for the various potentials considered, the fractional residuals with respect to the best fits for G(Q) from table 1 and the fractional residuals with respect to the numerical solution for the monomial potential.In summary, this figure demonstrates that G(Q) is independent on the inflaton potential chosen.
The results obtained are presented in figure 4 which is divided in 3 self-contained subfigures respectively for the cases c = {3, 1, −1}.Each sub-figure is in turn composed of an upper panel with the numerical solutions of G(Q) for the different inflaton potentials considered, and two lower panels showing the fractional residuals respectively in reference to the best fits for G(Q) from table 1 and to the numerical solution for the monomial potential.
For a given value of c, the scalar dissipation function G(Q) is independent of the inflaton potential.Specifically, the two lower panels of each sub-figure in figure 4 show that the fractional residuals in reference to the best fits for G(Q) and to the numerical solution for a monomial potential are respectively always O(1) and O(0.3), over the entire range of dissipation strength Q considered.In addition, one can note that the shape of the residuals with respect to the best fits for G(Q) reported in table 1 (and obtained using a monomial potential) is the same regardless of the inflaton potential considered, further emphasizing that such fits are universally valid and independent of the inflaton model implemented.To the best of our knowledge, this observation is shown here explicitly for the first time.The implication is very powerful: one can compute with our code the dissipation function G(Q) for a convenient choice of the inflaton potential, such as the monomial case, and apply it to any, perhaps more complex, inflaton model of interest.
In any case, one may still wonder if the results just presented depend on the specific values chosen for the parameters of the different inflaton potentials.The choice was conveniently made in order to guarantee that inflation can last for a sufficient number of e-folds and also roughly reproduce the amplitude of the scalar power spectrum measured by CMB data, i.e. ∆ 2 R (k * ) 2.1 × 10 −9 with k * = 0.05 Mpc −1 [83].Nevertheless, for a given inflaton potential the resulting scalar dissipation function is independent of the inflaton parameters chosen.This is expected as the main effect of a different choice of parameters is in the overall scale of the scalar power spectrum, which gets integrated out in order to obtain G(Q).For the interested reader, we refer to figure 8 in Appendix A which shows this explicitly for the monomial potential.
Number of e-folds before and after CMB horizon crossing: We assess the changes due to a variation in the number of e-folds after CMB horizon crossing N infl e and the quantity τ max that parametrizes the duration of the inflation stage considered before CMB horizon crossing.For this study, we assume no field dependence on the dissipation rate, i.e. m = 0, and we consider the monomial potential in eq.With respect to both scenarios analyzed, a change in the number of e-folds before or after CMB horizon crossing effectively results in a modification of the evolution of the back- ground quantities and/or the initial conditions in the perturbations' equations.As previously mentioned in section 4.3, changing the initial conditions has no appreciable impact on the solution to the perturbations' equations because the stochastic sources end up dominating the evolution, washing out any information regarding the initial state.On the other hand, changing the evolution of the background quantities in the perturbations' equations affects the resulting power spectrum, as one would expect.This evidently gets erased when the power spectrum resulting from the numerical computation is divided by the analytical estimate of the power spectrum to yield the scalar dissipation function G(Q).In other words, for a given value of Q, the different evolution of the background quantities only alters the scale of the scalar power spectrum which is then removed when evaluating G(Q).
Finally, it is useful to clarify here the relationship between the length of the τ interval N τ , the chosen value of τ max and the dissipation strength Q, in terms of how they impact the numerical stability of the solution of the perturbations' equations.Generally, as one would expect, the larger the value of τ max , the larger the value for N τ is required.The relationship between τ max and N τ is highly non-linear such that choosing τ max 8 requires setting N τ 10 7 which significantly slows down the code execution time.This makes such choice for τ max undesirable.In addition, we note two more relevant aspects: (1) for a given τ max , the smaller the value of Q, the larger N τ needs to be to fully capture the numerical solution of the perturbations; (2) a larger value of Q requires a larger value of τ max to allow the numerical solution to stabilize to its classical ensemble-average result.These considerations emphasize that varying N τ and τ max in terms of the dissipation strength Q would be the optimal choice from a computational standpoint, particularly if we want to extend the evaluation of the perturbations to values of the dissipation strength Q above ∼ 10 5 .For simplicity, in this work we choose a fixed value for τ max and N τ over the entire range of dissipation strength Q considered, which we have checked numerically to guarantee a stable solution.
Metric perturbations and corrections in H and η H : We now assess the impact of the metric perturbations (MPs) ϕ in eq.(3.1) and the corrections from the slow-roll parameters H and η H on the numerical results.In fact, as discussed in section 4.3, both MPs and the slow-roll corrections to the dynamics of the perturbations have been neglected in the previous analyses, however they can be reinserted through a switch that turns them back on in the numerical solution.For this study, we assume no field dependence on the dissipation rate, i.e. m = 0, and we consider the dynamics obtained from the monomial potential in eq.(4.1) with the index n = 4 and the self-interacting coupling constant λ = 10 −14 .We additionally set N infl e = 60.In the upper panel of figure 6 we compare the resulting scalar dissipation function G(Q) for a given temperature dependence in the dissipation rate c ∈ {3, 1, −1} without the metric perturbations and the corrections in H and η H (triangle), with the corrections from metric perturbations included in the computation (square) or with both corrections from MPs and the slow-roll parameters included (diamond).
Neither the metric perturbations nor the corrections in the slow-roll parameters leads to a significant impact on the scalar dissipation function G(Q) for the values of the index c considered.More specifically, the lower panel of figure 6 shows the fractional residuals between the case with ϕ = 0 and no H and η H corrections against the two other cases where ϕ = 0 and the H and η H corrections are included or not, see the caption of figure 6 for additional details.For all the values of c assessed and in the entire range of dissipation strength Q considered, the fractional residuals are O(0.25) and centered around zero.
In regards to the corrections from the slow-roll parameters, this result was trivially expected since the equations for the perturbations are evolved up to the time of CMB horizon crossing, that is N infl e = 60 e-folds before the end of inflation, when both H and η H are safely 1 and can therefore be safely neglected.On the other hand, concerning the metric perturbations, the coefficients of the ϕ and ϕ terms in the perturbations' equations are first order in the slow roll parameters and η and are thus negligible in the weak dissipative regime for Q < 1.In the strong dissipative regime (Q 1), while these terms should not be neglected they evidently have no impact in the resulting scalar power spectrum as the coefficients proportional to c multiplying δ ρr , exclusively drive the enhancement or suppression of the scalar power spectrum respectively when c is positive or negative.This is a general result that we observed repeatedly in this section: G(Q) only depends on the dissipation strength Q and the power-law temperature dependence in the dissipation rate c which, in the strong dissipative regime, uniquely drives the enhancement or the suppression of the scalar power spectrum.We quantify the precision of the analytical estimate of the scalar power spectrum from eq. (3.37) against our numerical solutions for the case in which the dissipation rate Υ does not depend on temperature and the inflaton field, i.e. {c, m} = {0, 0}. Figure 7 shows the analytical estimate and the numerical solution of the scalar power spectrum for different choices of the inflaton potential, normalized at Q = 10 −7 to the amplitude measured at recombination, namely ∆ 2 R (k * ) 2.1 × 10 −9 with k * = 0.05 Mpc −1 [83].Overall, the analytical approximation works well for the entire range of dissipation strength Q considered.Numerically, this is presented in the lower panel of figure 7, which shows the fractional residuals between the two methods for each inflaton potential used, and they are all O(1).Additionally, one can note that the residuals are always positive for Q O(1) and negative otherwise.This indicates that the analytical approximation marginally under(over)estimates the numerical solution of the scalar power spectrum for Q O(1) (Q O(1)).This behaviour is a direct consequence of the two main assumptions that allowed us to obtain the analytical form of the scalar power spectrum in eq.(3.37), namely: (1) we assumed the temperature T and Hubble parameter H to be constant with respect to z in order to bring them outside of the integral in eq.(3.33); (2) we neglected the term proportional to the slow-roll coefficient η by taking α ≈ ν as defined in eqs.(3.30) and (3.31). 17With respect to assumption (1), T and H are actually both monotonically increasing functions of z and therefore by computing the scalar power spectrum with respect to their value at z = 0, we are generally always underestimating the true value of the scalar power spectrum.This is particularly relevant for the quantum noise term whose amplitude is quadratically proportional to H. Thus, in the weak dissipative regime, i.e. for Q 1, when the quantum noise term still dictates the overall amplitude of the scalar power spectrum, the analytical approximation underestimates the value of the numerical solution, in agreement with our findings.On the other hand, assumption (2) is robust in the weak dissipative regime but becomes obsolete in the strong dissipative regime, for Q 1, when the slow-roll parameter η can be of O(Q).In this case, one can easily show that the approximated solution in eq.(3.34) to the integral in eq.(3.33), obtained by setting α ≈ ν, overestimates the value of the integral when α = ν.Particularly, in the strong dissipative regime, this effect is greater than the underestimation of the temperature T in the thermal noise term arising from assumption (1).Thus, overall from the above argument, we expect the analytical approximation to slightly overestimate the value of the numerical solution for Q 1, which is again in agreement with our findings.

Summary & Conclusion
In this paper we presented WarmSPy which consists-to our knowledge-the first open-source and publicly available code to compute the full, numerical scalar power spectrum for a given warm inflationary (WI) scenario and the scalar dissipation function G(Q) defined in eq.(3.38).The code is available at github.com/GabrieleMonte/WarmSPy.git.The stochastic nature of the problem is faced via a Monte Carlo approach that accounts for the thermal and quantum noise terms.We use the formalism above to assess the dependence of G(Q) on various model parameters and inflationary histories, as well as the effects of metric perturbations.We find that G(Q) only depends on (1) the dissipation strength Q, and (2) the dissipation rate's power-law temperature dependence c, see eq. (2.3), which in the strong dissipative regime Q 1 uniquely drives the enhancement or suppression of the scalar power spectrum.In this sense, the scalar dissipation function G(Q) can be thought of as an universal function: for a given value of c, G(Q) effectively describes the modification to the power spectrum in eq.(1.1) regardless of the inflaton model implemented, its history, and possible additional inflaton field dependencies in the dissipation rate or addition of metric perturbations.
Throughout this work, we specifically set c ∈ {3, 1, −1} motivated by the explicit constructions of dissipative rates in the literature.For all these cases, we report in table 1 the best-fit parameters for the proposed fitting functions of G(Q), respectively eqs.(4.9) and (4.8) for c ∈ {3, 1} and c = −1.Compared to previous work in the literature, the fits we present are valid in a much broader range of the dissipation strength, namely for Q ∈ [10 −7 , 10 4 ], with a minimal number of free parameters.The robustness and broader validity of our numerical fits for the function G(Q) will contribute to phenomenological studies of WI, and specifically help in precisely testing WI models against current and future CMB bounds.

1
Method validation and comparison with previous literature 17 5.2 Numerical fits of the scalar dissipation function G(Q) 19 5.3 Universality of the scalar dissipation function G(Q) 21 5.4 Robustness of the analytical estimate of the scalar power spectrum 28 1 Introduction )

14
Recall that we evaluate the background evolution in the interval Ne ∈ [0, N end e ] with N end e ≡ N pert e +N infl e .This means that the time of CMB horizon crossing, i.e.N inf e e-folds before the end of inflation, corresponds to Ne = N end e − N infl e = N pert e .In other words, τ (k * ) = 0 when Ne = N pert e .

Figure 1 :
Figure 1: The scalar dissipation function G(Q) appearing in eq.(3.38) as a function of the dissipation strength Q ≡ Γ/(3H), for different values of the set of indices (c, m) describing the power-law dependence of the dissipation rate to T and φ respectively, see eq. (2.3).Triangle: (c, m) = (3, 2); square: (c, m) = (1, 0); diamond: (c, m) = (−1, −2).Results are shown both for the code employing a constant value of Q (blue markers) and for an evolving Q (red markers).Also shown are the fitting functions obtained from refs.[47,66,80] (labeled "Fit 1" in the legend) and in ref.[39] (labelled "Fit 2" in the legend), see text for additional details.In summary, this figure demonstrates the consistency of our code in the resulting G(Q) when either keeping Q constant or allowing it to evolve, as well as its agreement with the well-established G(Q) fits found in the existing literature.

Figure 2 :
Figure 2: Upper panel: The numerical solutions obtained with our code for the scalar dissipation function G(Q) as a function of Q, along with an assessment of the goodness of the fitting functions in eqs.(4.8) and (4.9) with the fitting parameters reported in the caption.We have considered a monomial potential, see eq. (4.1), with n = 4 and λ = 10 −14 .Results are shown for a different temperature dependence in the dissipation rate (see eq. (2.3)): c = 3 (red), c = 1 (blue), c = −1 (green).Lower panel: The fractional residuals between the numerical result and the two different choices of the fit for each value of c, over the range of the dissipation strength Q considered.

Figure 3 :
Figure 3: Upper panel: An assessment of the impact of the index m in Eq. (2.3) on the numerical results of the scalar dissipation function G(Q), for different temperature dependence in the dissipation rate, namely c = 3 (red), c = −1 (green).Lower panel: The fractional residuals between the different choices of the index m over the range of the dissipation strength Q considered.In summary, this figure demonstrates that G(Q) is independent on the value of m.

Figure 4 :
Figure 4: The results of our numerical code for the various inflaton potentials discussed in section 4.2 and for a different temperature dependence in the dissipation rate: c = 3 (upper panel), c = 1 (middle panel), c = −1 (lower panel).Details for the potentials used are in the main text.For each panel, we show the numerical results of G(Q) for the various potentials considered, the fractional residuals with respect to the best fits for G(Q) from table 1 and the fractional residuals with respect to the numerical solution for the monomial potential.In summary, this figure demonstrates that G(Q) is independent on the inflaton potential chosen.
(4.1) with n = 4 and the self-interacting coupling constant λ = 10 −14 .The resulting scalar dissipation function G(Q) is compared for a given temperature dependence in the dissipation rate c ∈ {3, 1, −1} and different choices of the number of e-folds before and after CMB horizon crossing, respectively changing τ max and N infl e .Specifically, we present two main scenarios in the upper and lower plot of figure 5: (1) we fix τ max = 6.5 and vary the number of e-folds after CMB horizon crossing with N infl e = {40, 60, 80}; (2) we fix N infl e = 60 and vary the number of e-folds before CMB horizon crossing over which we solve the perturbations' equations by setting τ max = {6.5, 7, 7.5}.For a given value of c, the scalar dissipation function G(Q) is independent of the values for N infl e and τ max .Specifically, the two lower panels in figure 5 for each main plot show the fractional residuals between the different values of N infl e and τ max assessed for the range of indexes c = {3, 1, −1} considered, respectively.The residuals all lie below O(0.25) for the entire range of dissipation strength Q considered.

Figure 5 :
Figure 5: The scalar dissipation function G(Q) obtained for different choices of the number of e−folds before and after CMB horizon crossing, respectively changing τ max (lower panel) and N infl e (upper panel).Details on the values of τ max and N infl e can be found in the text.In each panel, we show the numerical results for G(Q) for different temperature dependence in the dissipation rate, c = 3 (red), c = 1 (blue), c = −1 (green), and the corresponding fractional residuals with respect to the fiducial case with N infl e = 60 and τ max = 6.5.In summary, this figure demonstrates that G(Q) is robust against changes in the inflationary history.

Figure 6 :
Figure 6: Upper panel: The scalar dissipation function G(Q) obtained from considering different sources of corrections to the inflaton perturbations: no corrections (triangle), corrections from metric perturbations (square), corrections from both metric perturbations and slow-roll parameters (diamond).Results are shown for the index c = 3 (red), c = 1 (blue), c = −1 (green).Lower panel: The fractional residuals between the results obtained with no corrections against i) setting only MPs different from zero or ii) both corrections from MPs and slow-roll parameters included in the computation.In summary, this figure demonstrates that G(Q) is independent on the metric perturbations and slow-roll corrections.

Figure 7 :
Figure 7: Upper panel: The scalar power spectrum for c = 0 as a function of the dissipation strength for different models for the inflaton potential: quartic (green), hilltop (yellow), natural inflation (blue), and β-exponential (red), see text for the choice of the parameters.Results are shown for the analytical expression in eq.(3.37) (solid lines) and for the computation with the numerical code (triangle marks).Lower panel: The fractional residuals from the two methods over the range of the dissipation strength Q considered.In summary, this figure demonstrates that the analytical estimate to the scalar power spectrum reproduces quite well the numerical results over the entire range of Q considered, with residuals always below O(1).